// VexCOL2CIR.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "stdio.h"
#include "math.h"
#include "stdlib.h"
#include <memory.h>



int main(int argc, char* argv[])
{

	struct RGBkuva
{
  unsigned char R;
  unsigned char G;
  unsigned char B;
};

struct VexPix
{
   unsigned short R;
   unsigned short G;
   unsigned short B;
   unsigned short IR;
};

	int i, j, index;
	FILE *FilePointer,*InFilePointer2, *OutFilePointer, *Out2FilePointer,*Out3FilePointer,*Out4FilePointer, *Out5FilePointer; 
// COL.TIF, PAN.TIF, outputs: RGBIR(uint16), RGB(CIR) uint8, PAN uint8, PANsharpCIR uint8, PAN (uint16)	
	FilePointer = fopen (argv[1],"rb");
	if (FilePointer == NULL) 
	{
		 printf("%s %s\n",argv[1], " Did not Open");
		 return 0;
	}
	InFilePointer2 = fopen (argv[2],"rb");
	if (InFilePointer2 == NULL) 
	{
		 printf("%s %s\n",argv[2], " Did not Open");
		 return 0;
	}
	OutFilePointer = fopen (argv[3],"wb");
		if (OutFilePointer == NULL) 
	{
		 printf("%s %s\n",argv[3], " Did not Open");
		 return 0;
	}

	Out2FilePointer = fopen (argv[4],"wb");
		if (Out2FilePointer == NULL) 
	{
		 printf("%s %s\n",argv[4], " Did not Open");
		 return 0;
	}
	Out3FilePointer = fopen (argv[5],"wb");
	if (Out3FilePointer == NULL) 
	{
		 printf("%s %s\n",argv[5], " Did not Open");
		 return 0;
	}

	   Out4FilePointer = fopen (argv[6],"wb");
	if (Out4FilePointer == NULL) 
	{
		 printf("%s %s\n",argv[6], " Did not Open");
		 return 0;
	}

	Out5FilePointer = fopen (argv[7],"wb");
	if (Out5FilePointer == NULL) 
	{
		 printf("%s %s\n",argv[7], " Did not Open");
		 return 0;
	}

	// uint(16)*4 Vexcel imahe
    VexPix *ima;
 	ima =  (VexPix *) malloc(3680*2400*8);
    unsigned short *RowOfReds, *RowOfGreens, *RowOfBlues, *RowOfIRs;
    RowOfReds =  (unsigned short *) malloc(3680*sizeof(unsigned short));
    RowOfGreens =  (unsigned short *) malloc(3680*sizeof(unsigned short));
    RowOfBlues =  (unsigned short *) malloc(3680*sizeof(unsigned short));
    RowOfIRs =  (unsigned short *) malloc(3680*sizeof(unsigned short));

	// uint(16)*3 RGB image
	RGBkuva *CIR_ima;
 	CIR_ima =  (RGBkuva *) malloc(3680*2400*3) ;

	double gammaIR = 0.6;
	double gammaR = 0.2;
	double gammaG = 0.2;
    
	// Vexcel PAN Image-versions uint(16) and uint(8)
    unsigned short *rawima;
	unsigned short *uint16ima;
    unsigned char  *byteima;
 	uint16ima =  (unsigned short *) malloc(90439680*sizeof(unsigned short));
	rawima =  (unsigned short *) malloc(90439680*sizeof(unsigned short));
 	byteima =  (unsigned char *) malloc(90439680*sizeof(unsigned char));

	if (rawima == NULL) {printf ("Muistinvaraus ei onnistunut\n");}
    if (rawima != NULL) {printf ("Muistinvaraus onnistui\n");}
	if (byteima == NULL) {printf ("Muistinvaraus ei onnistunut\n");}
    if (byteima != NULL) {printf ("Muistinvaraus onnistui\n");}
	if (ima == NULL) {printf ("Muistinvaraus ei onnistunut\n");}
    if (ima != NULL) {printf ("Muistinvaraus onnistui\n");}
	if (RowOfReds == NULL) {printf ("Muistinvaraus ei onnistunut\n");}
    if (RowOfReds != NULL) {printf ("Muistinvaraus onnistui\n");}
	if (RowOfGreens == NULL) {printf ("Muistinvaraus ei onnistunut\n");}
    if (RowOfGreens != NULL) {printf ("Muistinvaraus onnistui\n");}
	if (RowOfBlues == NULL) {printf ("Muistinvaraus ei onnistunut\n");}
    if (RowOfBlues != NULL) {printf ("Muistinvaraus onnistui\n");}
	if (RowOfIRs == NULL) {printf ("Muistinvaraus ei onnistunut\n");}
    if (RowOfIRs != NULL) {printf ("Muistinvaraus onnistui\n");}
	if (CIR_ima == NULL) {printf ("Muistinvaraus ei onnistunut\n");}
    if (CIR_ima != NULL) {printf ("Muistinvaraus onnistui\n");}


 // Read PAN image
	fseek(InFilePointer2,long(8), SEEK_SET);
    fread(&rawima[0], sizeof(unsigned short), 90439680, InFilePointer2);
    unsigned short maxvalue, minvalue;
	unsigned char maxchar;
	maxvalue = 0;
	minvalue = 20000;
	for (i=0;i<=90439680-1;i++)
	{
		if (rawima[i]>=maxvalue) {maxvalue = rawima[i];}
        if (rawima[i] > 1 && rawima[i]<=minvalue) {minvalue = rawima[i];}
	}
	printf("%s %d\n", "PAN Maximum: ", maxvalue);
	printf("%s %d\n", "PAN Minimum: ", minvalue);

   fclose(InFilePointer2);
   int l, hh, tile, base, rivi, sarake, indeksi;
   int TS = 512;
   int TS2;
   TS2 = TS*TS;
   double gamma = 0.15;
   double Gcorr;

   Gcorr = double(maxvalue) / (pow(double(maxvalue), gamma));
   maxvalue = unsigned short(pow(double(maxvalue), gamma) * Gcorr);
   minvalue = unsigned short(pow(double(minvalue), gamma) * Gcorr);

   // Remove tiling
  printf ("Tiilityksen purku\n");
  printf("%s %s\n", "With simoultaneous output of ", argv[7]);
   // Output uint16 PAN image; ima[] before correction!
  fseek(Out5FilePointer,long(0),SEEK_SET);

  for (l=1; l<=15;l++)  // tilerows
    {
     for (hh=1; hh<=23; hh++) // tilecolumns
	 {
  	   tile = tile + 1;    
	   base = (tile-1) * 512 * 512; 
	   for (i=1;i<=512;i++)  
		 {
   	       for (j=1; j<=512;j++) 
			{

             rivi = (l - 1) * 512 + j; // row in the untiled image
	         sarake = (hh - 1) * 512 + i; // column in untiled image
			 indeksi = (rivi-1) * 11776 + sarake-1;
			 index = i+(j-1)*TS+(hh-1)*TS2+(l-1)*23*TS2-1;
			 // Gamma correct
             uint16ima[indeksi] = rawima[index];
			 rawima[index] = unsigned short(pow(double(rawima[index]), gamma) *Gcorr );
             byteima[indeksi] = unsigned char( floor( (double(rawima[index])-double(minvalue)) / double(double(maxvalue)-double(minvalue)) * double(255)));
			}
		 }
	 }
	 }

  printf("%s %s\n", "Done, now writing output of ", argv[5]);
  // Print all rows in image
  int ii, jj, kk,kkk;
  printf ("Kopioin PAN datan\n");
  free(rawima);
  unsigned short  *Pan_ima;
  Pan_ima =  (unsigned short *) malloc(86250000*2);

  for (ii=1;ii<=7500;ii++)
  {
	  kk=(ii-1)*(11500);
	  jj=(ii-1)*(11776);
      fwrite(&uint16ima[jj], long(2), 11500, Out5FilePointer);
	  fwrite(&byteima[jj],1,11500,Out3FilePointer);
	  for (kkk=0;kkk<=11499;kkk++)
	  {
	    Pan_ima[kk+kkk]=byteima[jj+kkk];
      }
  }

  free(uint16ima);
  free(byteima);
//  free(RowOfReds);
//  free(RowOfGreens);
//  free(RowOfBlues);
//  free(RowOfIRs);
  printf("%s\n", "Copying of PAN Done! Proceeding with CIR, reading it");
  fclose(Out3FilePointer);
  fclose(Out5FilePointer);
//  return 0;

  VexPix maxs;
	maxs.R = 0; maxs.G = 0; maxs.B = 0; maxs.IR = 0; 
    VexPix mins;
	mins.R = 10000; mins.G = 10000; mins.B = 10000; mins.IR = 10000; 
    double GcorrIR, GcorrR, GcorrG;

for (i=1;i<=2400;i++)
{
	fseek(FilePointer,(8 + 3680 * i * 8 - 3680 * 8),SEEK_SET);
    fread(&RowOfReds[0], 2, 3680, FilePointer);
//	fseek(FilePointer,long(8 + 3680 * i * 8 - 3680 * 6),SEEK_SET);
    fread(&RowOfGreens[0], 2, 3680, FilePointer);
//	fseek(FilePointer,long(8 + 3680 * i * 8 - 3680 * 4),SEEK_SET);
    fread(&RowOfBlues[0], 2, 3680, FilePointer);
//	fseek(FilePointer,long(8 + 3680 * i * 8 - 3680 * 2),SEEK_SET);
    fread(&RowOfIRs[0], 2, 3680, FilePointer);
    // Loop each row, assign and find minima and maxima
	for (j=1;j<=3680;j++)
	{
	index = (i-1)*(3680)+j-1;
    ima[index].R = RowOfReds[j-1];
	if (ima[index].R > maxs.R) {maxs.R = ima[index].R;}   
	if (ima[index].R < mins.R) {mins.R = ima[index].R;}   
    ima[index].G = RowOfGreens[j-1];
	if (ima[index].G > maxs.G) {maxs.G = ima[index].G;}   
	if (ima[index].G < mins.G) {mins.G = ima[index].G;}   
    ima[index].B = RowOfBlues[j-1];
	if (ima[index].B > maxs.B) {maxs.B = ima[index].B;}   
	if (ima[index].B < mins.B) {mins.B = ima[index].B;}   
    ima[index].IR = RowOfIRs[j-1];
	if (ima[index].IR > maxs.IR) {maxs.IR = ima[index].IR;}   
	if (ima[index].IR < mins.IR) {mins.IR = ima[index].IR;}   
	}
}

free(RowOfReds);
free(RowOfGreens);
free(RowOfBlues);
free(RowOfIRs);

// Output the RGBIR image before scaling
fseek(OutFilePointer,long(0),SEEK_SET);
fwrite(&ima[0].R, long(8), 3680*2400, OutFilePointer);

printf("%s\n", "RGBIR uint(16) Written! Proceeding with CIR uint8");

// GammaCorrect the CIR image

    GcorrIR = double(maxs.IR) / (pow(double(maxs.IR), gammaIR));
    GcorrR = double(maxs.R) / (pow(double(maxs.R), gammaR));
    GcorrG = double(maxs.G) / (pow(double(maxs.G), gammaG));
	maxs.IR = unsigned short(pow(double(maxs.IR), gammaIR) * GcorrIR);
	mins.IR = unsigned short(pow(double(mins.IR), gammaIR) * GcorrIR);
	maxs.R = unsigned short(pow(double(maxs.R), gammaR) * GcorrR);
	mins.R = unsigned short(pow(double(mins.R), gammaR) * GcorrR);
	maxs.G = unsigned short(pow(double(maxs.G), gammaG) * GcorrG);
	mins.G = unsigned short(pow(double(mins.G), gammaG) * GcorrG);

for (i=1;i<=2400;i++)
{
	for (j=1;j<=3680;j++)
	{
	index = (i-1)*(3680)+j-1;
// Gamma korjaus
  //  Gcorr = double(16384) / (pow(double(16384), gammaIR));
	ima[index].IR = unsigned short(pow(double(ima[index].IR), gammaIR) * GcorrIR);
  //  Gcorr = double(16384) / (pow(double(16384), gammaR));
	ima[index].R = unsigned short(pow(double(ima[index].R), gammaR) * GcorrR );
  //  Gcorr = double(16384) / (pow(double(16384), gammaG));
	ima[index].G = unsigned short(pow(double(ima[index].G), gammaG) * GcorrG );

	CIR_ima[index].R = unsigned char(floor( (double(ima[index].IR) - double(mins.IR))/(double(maxs.IR)-double(mins.IR))*double(255))); 
	CIR_ima[index].G = unsigned char(floor( (double(ima[index].R) - double(mins.R))/(double(maxs.R)-double(mins.R))*double(255))); 
	CIR_ima[index].B = unsigned char(floor( (double(ima[index].G) - double(mins.G))/(double(maxs.G)-double(mins.G))*double(255)));

	}
}

printf ("%d %d %d %d %d %d %d %d\n",mins.R,mins.G,mins.B,mins.IR,maxs.R,maxs.G,maxs.B,maxs.IR);

// Now ima[] has the R,G,B,IR -structure of unsigned shorts
// We output it and RGB image of the CIR image

// Print the RGB-image of uint8
fseek(Out2FilePointer,long(0),SEEK_SET);
fwrite(&CIR_ima[0].R, long(3), 3680*2400, Out2FilePointer);
// Print the RGBIR-image of uint16

printf("%s\n", "RGB CIR uint8 Written! Proceeding with PANsharpened CIR in RGB uint8");

fclose(FilePointer);
fclose(OutFilePointer);
fclose(Out2FilePointer);


free(ima);

// PAN sharpen for CIR
// We have Pan_ima[] with PAN    11776 x 2400
//         CIR_ima[] win CIR-RGB 3680 x 2400

  // Allocate space for the pansharpened CIR image
  RGBkuva  *PanCIR;
  PanCIR =  (RGBkuva *) malloc(86250000*3);
  double min, max, delta;
  int col, row, k;
  double R, G, B, Pi, h, s, v;
  double f, p, q, t;
  Pi=3.141592654;
  double mini, intensity;
  
  for (i=1;i<=11500;i++)
  {
	  for (j=1;j<=7500;j++)
	  {
         col = int(ceil(double(i)/3.125));
         row = int(ceil(double(j)/3.125));

//		 printf("%d %d\n", col, row);

         R = double(CIR_ima[(row-1)*3680+col-1].R)/double(255);
         G = double(CIR_ima[(row-1)*3680+col-1].G)/double(255);
         B = double(CIR_ima[(row-1)*3680+col-1].B)/double(255);
         
		 intensity = double(Pan_ima[(j-1)*11500+i-1])/double(255);
    
		 mini = R;
	     max = R;
	
	     if (G>max) {max=G;}
         if (B>max) {max=B;}
     	 if (G<mini) {mini=G;}
         if (B<mini) {mini=B;}

	     v = max;				// v
	     delta = max - mini;
	    if( max != 0 )
		  s = delta / max;		// s
	     else {
		 // r = g = b = 0		// s = 0, v is undefined
		 s = 0;
		 h = -1;
		 goto donehsv;
		 }
		 if( R == max )
			 h = ( G - B ) / delta;		// between yellow & magenta
		 	else if( G == max )
				h = 2 + ( B - R ) / delta;	// between cyan & yellow
				else
					h = 4 + ( R - G ) / delta;	// between magenta & cyan
				    h = h*60;				// degrees
					if( h < 0 )
						h += 360;

donehsv:;

// NOw compute back
 v = intensity;
 if( s == 0 ) {
		// achromatic (grey)
		R = G = B = intensity;
		goto done;
	}


	h /= 60;			// sector 0 to 5

	k = int(floor( h ));
	f = h - k;			// factorial part of h
	p = v * ( 1 - s );
	q = v * ( 1 - s * f );
	t = v * ( 1 - s * ( 1 - f ) );
	switch( k ) {
		case 0:
			R = v;
			G = t;
			B = p;
			break;
		case 1:
			R = q;
			G = v;
			B = p;
			break;
		case 2:
			R = p;
			G = v;
			B = t;
			break;
		case 3:
			R = p;
			G = q;
			B = v;
			break;
		case 4:
			R = t;
			G = p;
			B = v;
			break;
		default:		// case 5:
			R = v;
			G = p;
			B = q;
			break;
	}




done:;

         PanCIR[(j-1)*11500+i-1].R = unsigned char(R*double(255));
         PanCIR[(j-1)*11500+i-1].G = unsigned char(G*double(255));
         PanCIR[(j-1)*11500+i-1].B = unsigned char(B*double(255));
			  

	  }
  }

printf("%s\n", "Wrinting PANsharpened CIR");

fseek(Out4FilePointer,long(0),SEEK_SET);
fwrite(&PanCIR[0].R, long(3), 11500*7500, Out4FilePointer);

free(CIR_ima);
free(PanCIR);


	return 0;
}