#include "ninox.h"

static double * Comp;
static struct Image *Ref;

// Support for detecting and masking out dead pixels
#define MAX_DEAD_PIXELS 512
static int DeadPixels = 0;
static int DeadPixelList[MAX_DEAD_PIXELS];

int
LoadDarkFrameRef(char *fname)
	{
	int npix,i;
	unsigned short *uptr;
	double total=0,avg;
	struct Image *tmp;
	int DeadThreshHold = 25 * 256;

	tmp = LoadImage(fname,"darkreference");
	if (! tmp) {
	   Print("Error loading DarkFrame image '%s'\n",fname);
	   exit(1);
	   }

	// Make sure it's 16bpp monochrome
	Ref = ConvertImage(tmp,IMG_FIT,16);

	Print("[DarkFrame image: %dx%dx%d ",Ref->width,Ref->height,Ref->depth);
	npix = Ref->width * Ref->height;

	if (! npix) {
	   Print("Error loading DarkFrame image '%s': no pixels?\n",fname);
	   exit(1);
	   }

	uptr = (unsigned short *) Ref->data;

	if (Comp) free(Comp);
	Comp = Malloc(npix * sizeof(double));

	DeadPixels = 0;
	for(i=total=0; i<npix; ++i) {
	   if (uptr[i] >= DeadThreshHold && DeadPixels < MAX_DEAD_PIXELS) DeadPixelList[DeadPixels++] = i;
	   else total += uptr[i];
	   }

	avg = (double)total / (npix - DeadPixels);
	Print("Average: %f, %d dead pixels found\n",avg,DeadPixels);

	for(i=0; i<npix; ++i) Comp[i] = (double)uptr[i];

	Print("LoadDarkRef: Loaded %d entries from '%s'\n",npix,fname);

	return(npix);
	}

int
SubtractDarkFrame(struct Image *img)
	{
	int dc,i;
	unsigned char *ptr = (unsigned char *)img->data;
	unsigned short *uptr = (unsigned short *)img->data;
	double val;
	int w = img->width;
	int h = img->height;
	int npix = w * h;

	if (Ref->width == 0 || Ref->height == 0) {
	   Print("SubtractDarkFrame: Error: no reference loaded\n");
	   return 0;
	   }

	if (img->width != Ref->width || img->height != Ref->height) {
	   Print("SubtractDarkFrame: Image has different dimensions to reference\n");
	   return(0);
	   }

Print("Subtracting...\n");

	// Dark Frame subtraction
	switch(img->depth) {
	   case 8:
		for(i=0; i<npix; ++i) if (ptr[i]) {
			double d = Comp[i];
			double val = ptr[i];

			if (Ref->depth==16) d /= 256;
			val -= d;
			if (val<0) val=0; if (val>255) val=255;
			ptr[i] = val;
			}
		break;

	   case 16:
		for(i=0; i<npix; ++i) if (uptr[i]) {
			double d = Comp[i];
			double val = uptr[i];

			val -= d;
			if (val<0) val=0; if (val>65535) val=65535;
			uptr[i] = val;
			}

		break;
	   default:
		Print("SubtractDarkFrame: Unsupported depth %d\n",img->depth);
		return 0;
	   }

	// Dead pixel mapping
	switch(img->depth) {
	   case 8:
		for(i=0; i< DeadPixels; ++i) {
		   int n = DeadPixelList[i];
		   int y = n / w;
		   int x = n % w;
		   double t;
		   if (x>0 && y>0 && x < w-1 && y < h-1) {
		   	t = ptr[n-w-1]; t += ptr[n-w]; t += ptr[n-w+1];
		   	t += ptr[n-1]; t += ptr[n+1]; t += ptr[n+w-1];
		   	t += ptr[n+w]; t += ptr[n+w+1];
		   	t /= 8; if (t<0) t=0; if (t>255) t=255;
		   	ptr[n] = (int) t;
		   	}
		   }
		break;

	   case 16:
		for(i=0; i< DeadPixels; ++i) {
		   int n = DeadPixelList[i];
		   int y = n / w;
		   int x = n % w;
		   double t;
		   if (x>0 && y>0 && x < w-1 && y < h-1) {
		   	t = uptr[n-w-1]; t += uptr[n-w]; t += uptr[n-w+1];
		   	t += uptr[n-1]; t += uptr[n+1]; t += uptr[n+w-1];
		   	t += uptr[n+w]; t += uptr[n+w+1];
		   	t /= 8; if (t<0) t=0; if (t>65535) t=65535;
		   	uptr[n] = (unsigned short)(t/8);
		   	}
		   }
		break;
		}

	return(1);
	}
