#include "ninox.h"
#include <ctype.h>
#include <string.h>

static unsigned int _SubSample8(unsigned char *ptr, int img_wid, int bpp, int x_size, int y_size);
static unsigned int _SubSample16(unsigned char *ptr, int img_wid, int bpp, int x_size, int y_size);

/*
 * Generate and return an array of subsampled images, count is how many images to make, 
 * and is also the number of elements of xsizes[] and ysizes[] which give the sizes
 * of the subsampled images.
 *
 * rdepth is the requested depth, the list must contain sub-images of this depth
 * 
 * Region(r) defines the area in the original image to use. If NULL then the cutout region in
 * the source image is used
 */
struct Image **
GenerateSubSampledImages(struct Image *img, struct Region *r, int *xsizes, int *ysizes, int rdepth, int count)
	{
	struct Image **list;
	unsigned char *buffer = img->data;
	int width = img->width;
	int height = img->height;
	int depth = img->depth;
	int w,h,x,y,x1,y1,x2,y2;
	int subsample_x,subsample_y;
	int i,bpp,n,x_inc;
	int x_last,y_last;
	int pixels_per_row;
	unsigned char *buf;
	unsigned short *ubuf;
	char fname[1024];

	if (depth != 16 && depth != 8 && depth != 24) {
	   Print("GenerateSubSampledImages: Depth must be 8,16 or 24, not %d\n",depth);
	   return NULL;
	   }

	if (count>32) {
	   Print("GenerateSubSampledImages: maximum of 32 images allowed\n");
	   return NULL;
	   }

	list = (struct Image **)malloc(sizeof(struct Image **) * count);

	if (r != NULL) {
	   x1=r->x1; y1=r->y1;
	   x2=r->x2; y2=r->y2;
	   }
	else {
	   x1 = img->cutout.x1; x2 = img->cutout.x2;
	   y1 = img->cutout.y1; y2 = img->cutout.y2;
	   }

	Clip(width,height,&x1,&y1,&x2,&y2);

	pixels_per_row = x2-x1+1;
	bpp = depth/8;

//WriteImage(img);

	for(i=0; i<count; ++i) {

	   /* Desired width and height */
	   w = xsizes[i]; h=ysizes[i];

	   /* size of subsampled pixel */
	   subsample_x = width/w;
	   subsample_y = height/h;

	   x_last = x2 - subsample_x + 1;
	   y_last = y2 - subsample_y + 1;
	   x_inc = subsample_x * bpp;

	   // Allocate buffer for result
	   buf = (unsigned char *)malloc(w * h * rdepth/8);
	   ubuf = (unsigned short *)buf;

	   if (buf == NULL) {
	     Print("GenerateSubSampledImages: Out of memory");
	     exit(1);
	     }

	   switch(rdepth) {
		case 8:
	   		for(y=y1,n=0; y <= y_last; y += subsample_y) {
	      		   unsigned char *ptr = buffer + (y*width+x1) * bpp;
	      		   for(x=x1; x <= x_last; x+=subsample_x, ptr += x_inc) {
		  		buf[n++] = _SubSample8(ptr,width,bpp,subsample_x,subsample_y);
		  		}
	     		   }
			break;
		case 16:
	   		for(y=y1,n=0; y <= y_last; y += subsample_y) {
	      		   unsigned char *ptr = buffer + (y*width+x1) * bpp;
                           int rowpixels=0;
                           for(rowpixels=0,x=x1; x <= x_last; x+=subsample_x, ptr += x_inc, ++rowpixels) {
                                ubuf[n++] = _SubSample16(ptr,width,bpp,subsample_x,subsample_y);
                                }
                           while(rowpixels++ < w) ubuf[n++]=0;  // Pad with black
                           }
			break;
		case 24:
			Print("GenerateSubSampledImages: Depth %d not supported\n",rdepth);
			exit(1);
			break;
		}

	   /* Generate image */
	   sprintf(fname,"%s_%d.fit",img->src_fname,i);
	   list[i] = CreateFitsImage(fname, w, h, rdepth, buf, -1, -1.0);

//ShowImage(list[i]);
//Usleep(5000*1000);
	   if (! list[i]) {
	  	Print("GenerateSubSampledImages: error generating file %s\n",fname);
		exit(1);
		}
	   }

	return list;
	}

/*
 * Subsample a region starting at *ptr of size X size pixels.
 * Always return a 16bpp result
 */
static unsigned int
_SubSample16(unsigned char *ptr, int img_wid, int bpp, int x_size, int y_size)
	{
	int x,y,val=0;
	unsigned short v,*uptr;
	double d;
	
	switch(bpp) {
	   case 1:
		for(y=0; y<y_size; ++y) {
	   	   for(x=0; x<x_size; ++x) val += ((unsigned short)ptr[x]) << 8;
		   ptr += img_wid;
	   	   }
		break;

	   case 2:
		uptr = (unsigned short *)ptr;
		for(y=0; y<y_size; ++y) {
	   	   for(x=0; x<x_size; ++x) val += uptr[x];
		   uptr += img_wid;
	   	   }
		break;

	   case 3:
		// 24bpp data, 8 bit BGR, convert to monochrome
		for(y=0; y<y_size; ++y) {
	   	   for(x=0; x<x_size*3; x+=3) {
			d = ((double)ptr[x])*0.114;
			d += ((double)ptr[x+1])*0.587;
			d += ((double)ptr[x+2])*0.299;
			val += ((unsigned short)(d))<<8;
			}
		   ptr += img_wid*3;
	   	   }
		break;

	   default:
		Print("SubSample: Unsupported depth %d\n",bpp);
		exit(1);
	   }

	return val / (x_size * y_size);
	}

/*
 * Subsample a region starting at *ptr of size X size pixels.
 * Always return a 8bpp pixel
 */
static unsigned int
_SubSample8(unsigned char *ptr, int img_wid, int bpp, int x_size, int y_size)
	{
	int x,y,val=0;
	unsigned short v,*uptr;
	double d;
	
	switch(bpp) {
	   case 1:
		for(y=0; y<y_size; ++y) {
	   	   for(x=0; x<x_size; ++x) val += ptr[x];
		   ptr += img_wid;
	   	   }
		break;

	   case 2:
		uptr = (unsigned short *)ptr;
		for(y=0; y<y_size; ++y) {
	   	   for(x=0; x<x_size; ++x) val += uptr[x]>>8;
		   uptr += img_wid;
	   	   }
		break;

	   case 3:
		// 24bpp data, 8 bit BGR, convert to monochrome
		for(y=0; y<y_size; ++y) {
	   	   for(x=0; x<x_size*3; x+=3) {
			d = ((double)ptr[x])*0.114;
			d += ((double)ptr[x+1])*0.587;
			d += ((double)ptr[x+2])*0.299;
			val += d;
			}
		   ptr += img_wid*3;
	   	   }
		break;

	   default:
		Print("SubSample: Unsupported depth %d\n",bpp);
		exit(1);
	   }

	return val / (x_size * y_size);
	}

// Test mode: subsample input image n times with a random offset between samples
// Write out samples as images named "sample_n.fit"

int
Test_Subsample(struct Image *img)
	{
	// calculate the (width,height) of the grid in source pixels
	struct Image **sub;
	struct fits_info *fi;
	struct Region r;
	int offset_x, offset_y;
	int offset_delta = Subsample_max_offset - Subsample_min_offset;
	int i,j;
	int output_width,output_height;
	unsigned short *uptr;
	int v;
	int *num;	// Array to remap file names to randomise the sequence
	FILE *z;

	Print("Test Mode: Subsample\n");
	stack_init();

	// initialise random number generator
	Srandom();

	// Standard type is 16bpp monochrome
	img = ConvertImage(img,IMG_FIT,16);

	img->dst_fname = strdup("reference.fit");
	WriteImage(img);

	z = fopen("offsets.txt","w");

	num = (int *)ZeroMalloc(Subsample_nframes * sizeof(int));
	for(i=0; i<Subsample_nframes; ++i) num[i] = i;
	for(i=0; i<Subsample_nframes; ++i) {
	   int j = Random(Subsample_nframes);
	   int tmp = num[i]; num[i] = num[j]; num[j] = tmp;
	   }

	for(i=0; i<Subsample_nframes; ++i) {
	   char fname[32];

	   sprintf(fname,"%s%d.fit",Subsample_prefix,num[i]);
	   if (! Subsample_scan) {
		// random offsets
	      	offset_x = Subsample_min_offset + offset_delta * ((double)Random(10000)/(double)10000.0);
	      	offset_y = Subsample_min_offset + offset_delta * ((double)Random(10000)/(double)10000.0);
	 	}
	   else {
		// non-random, scan through all values
		int mod = Subsample_max_offset + 1;
		offset_x = Subsample_min_offset + (i % mod);
		offset_y = Subsample_min_offset + ((i/mod) % mod);
		}

	   r.x1 = offset_x;
	   r.y1 = offset_y;
	   r.x2 = img->width-1;
	   r.y2 = img->height-1;

	   // Calculate the size of the subsampled image
	   output_width = img->width / Subsample_pixel_width;
	   output_height = img->height / Subsample_pixel_height;

	   Print("Generating image %d: %s (%d,%d) with offsets (%d,%d)\n",
		i,fname,
		output_width, output_height,
		offset_x,offset_y);

	   // Generate 1 subsampled image
	   sub = GenerateSubSampledImages(img, &r, &output_width, &output_height, 16, 1);
	   sub[0]->dst_fname = strdup(fname);

	   if (Subsample_noise > 0) {
	      // Add small amount of noise to image
	      uptr = (unsigned short *)(sub[0]->data);
	      for(j=0; j<output_width * output_height; ++j) {
	   	v = *uptr + ((double)Random(RAND_MAX)/(double)RAND_MAX) * (*uptr) * (double)Subsample_noise/100.0;
		if (v>65535) v=65535;
		*(uptr++) = v;
		}
	      }

	   if (Subsample_upscale>1) {
	      upscale_image(sub[0],Subsample_upscale);
	      }

	   if (Subsample_smooth > 1) {
      	      smooth_image(sub[0],Subsample_smooth);
	      }

	   sub[0] = ConvertToType(sub[0],IMG_BMP,8);
	   WriteImage(sub[0]);
	   stack_frame(sub[0]);

	   DestroyImage(sub[0]);

	   fprintf(z,"%s: (%d,%d)\n",fname,offset_x, offset_y); fflush(z);
	   }

	fclose(z);

	write_stack_file("stack.fit");
	return 1;
	}