#define cimg_debug 2

#include "CImg.h"
#include <ctype.h>
#include <valarray>
#include <assert.h>
#include <math.h>
#define ACTG char actg[256]; memset(actg, 255, 256); \
	actg[(int)'A']=(char)0; actg[(int)'C']=(char)1; actg[(int)'T']=(char)2; actg[(int)'G']=(char)3; 

using namespace cimg_library;

class CGenome {
public:
        CGenome(char* filename) {
                ACTG; char ch; unsigned int pos, fileLen, read, i;
                FILE* file = fopen(filename, "rb"); if (!file) { printf("Can't open %s\n", (filename)?filename:"<null>"); }
                fseek(file, 0, SEEK_END); fileLen = ftell(file); assert(fileLen);
                m_data = CImg<unsigned char>(fileLen,1); fseek(file, 0, SEEK_SET);
                read = (unsigned int)fread(m_data.data, 1, fileLen, file); assert(read == fileLen);
                
		pos = 0; i=0; 
		if (m_data[0] == '>') while (m_data[i++] != '\n') ; else printf("Not fasta format - %s.\n",filename);
                while (i<fileLen) { 
                        ch = m_data[i++];
                        if(ch=='>') {printf("Two or more sequence found.\n"); break; }
                        if (actg[(unsigned char)ch]>=0) m_data[pos++] = ch;
                }
                m_data.resize(pos,1,1,1,0);
        }
        char *getData()  { return (char*)m_data.data; }
        unsigned int getLength() { return m_data.width; }
private:
        CImg<unsigned char> m_data;
};

char buf[1024];

void gen2asc(char* filename, int len) {
	char *buffer; unsigned int i, width = 1<<len;

	CGenome gen(filename);
	char* data = gen.getData();

	CImg<float> src = CImg<>(width,width).fill(0); unsigned int X = 0, Y = 0;
	for (i=0; i<gen.getLength(); i++) {
		buffer = (data+i);

		char ch = toupper(*buffer);
		int z = (ch=='A')?0:((ch=='C')?1:((ch=='T')?2:((ch=='G')?3:-1)));
		if (z>=0) { X*=2; X+=z%2; X%=width; Y*=2; Y+=z/2; Y%=width; src(X,Y)=src(X,Y)+1;	}
	}

	*buf = 0; assert(strlen(filename)<1000); 
	sprintf(buf, "%s.%d.asc", filename, len); src.save(buf);
	std::printf("%s\n", buf);
}

void usage(char* filename) {
	std::fprintf(stderr, "Usage: %s [-n] <fasta files>\n\tn=1..9\n\n", filename);
	exit(-1);
}

int main(int argc,char **argv) {
	int len = 9, from = 1;
	if (argc < 2)
		usage(argv[0]);

	if (argv[1][0] == '-') {
		len = atoi(argv[1]+1);
		if (len<1)
			usage(argv[0]);
		from = 2;
	}
	
	for (int i=from; i < argc; i++) {
		gen2asc(argv[i], len);
	}

	return 0;
}