add temperature compensation + new colormap for better false color generation

image.c

@@ -23,18 +23,18 @@
 #include <stdio.h>
 #include <string.h>
 #include <sndfile.h>
-
+#include <math.h>
 
 typedef struct {
 double slope;
 double offset;
 } rgparam;
 
-static void rgcomp(double *x, rgparam *rgpr)
+static void rgcomp(double x[16], rgparam *rgpr)
 {
-/*const double y[9] = { 0.106,0.215,0.324,0.434,0.542,0.652,0.78,0.87 ,0.0}; */
-const double y[9] = { 0.11872,0.2408,0.36288,0.48608,0.60704,0.73024,0.8736,0.9744 , 0.0 }; 
-const double yavg=0.48819556;
+/* 0.106,0.215,0.324,0.434,0.542,0.652,0.78,0.87 ,0.0 */
+const double y[9] = { 31.1,63.0,95.0,127.2,158.9,191.1,228.6,255.0, 0.0 }; 
+const double yavg=(y[0]+y[1]+y[2]+y[4]+y[5]+y[6]+y[7]+y[8])/9.0;
 double xavg;
 double sxx,sxy;
 int i;
@@ -59,19 +59,22 @@ static double rgcal(float x,rgparam *rgpr)
 return(rgpr->slope*x+rgpr->offset);
 }
 
+
+static double tele[16];
+static nbtele;
+
 int Calibrate(float **prow,int nrow,int offset)
 {
 
-double tele[3000];
-double frm[9];
+double teleline[3000];
+double wedge[16];
 rgparam regr[30];
 int n;
-int k,mk;
-int shift;
-int mshift,channel;
+int k;
+int mtelestart,telestart;
+int channel=-1;
 float max;
 
-
 printf("Calibration ");
 fflush(stdout);
 
@@ -79,65 +82,63 @@ fflush(stdout);
 for(n=0;n<nrow;n++) {
 	int i;
 
-	tele[n]=0.0;
+	teleline[n]=0.0;
 	for(i=3;i<43;i++)
-		tele[n]+=prow[n][i+offset+909];
-	tele[n]/=40;
+		teleline[n]+=prow[n][i+offset+909];
+	teleline[n]/=40;
 }
 
-if(nrow<128) {
+if(nrow<192) {
 	fprintf(stderr," impossible, not enought row\n");
 	return (0);
 }
 
 /* find telemetry start */
-max=-1000.0;mshift=0;
-for(n=60;n<nrow-3;n++) {
+max=0.0;mtelestart=0;
+for(n=nrow/3-64;n<2*nrow/3-64;n++) {
 	float df;
 
-	df=(tele[n-4]+tele[n-3]+tele[n-2]+tele[n-1])-(tele[n]+tele[n+1]+tele[n+2]+tele[n+3]);
-	if(df> max){ mshift=n; max=df; }
+	df=(teleline[n-4]+teleline[n-3]+teleline[n-2]+teleline[n-1])/(teleline[n]+teleline[n+1]+teleline[n+2]+teleline[n+3]);
+	if(df> max){ mtelestart=n; max=df; }
 }
 
-mshift-=64;
-shift=mshift%128;
-if(nrow<shift+8*9) {
+mtelestart-=64;
+telestart=mtelestart%128;
+
+if(mtelestart <0 || nrow<telestart+128) {
 	fprintf(stderr," impossible, not enought row\n");
 	return (0);
 }
 
-/* compute regression params */
-for(n=shift,k=0;n<nrow-8*9;n+=128,k++) {
+/* compute wedges and regression */
+for(n=telestart,k=0;n<nrow-128;n+=128,k++) {
 	int j;
 
-	for(j=0;j<9;j++) {
+	for(j=0;j<16;j++) {
 		int i;
 
-		frm[j]=0.0;
+		wedge[j]=0.0;
 		for(i=1;i<7;i++)
-			frm[j]+=tele[n+j*8+i];
-		frm[j]/=6;
+			wedge[j]+=teleline[n+j*8+i];
+		wedge[j]/=6;
 	}
-	rgcomp(frm,&(regr[k]));
 
+	rgcomp(wedge,&(regr[k]));
 
-	if(n==mshift) {
+	if (k==1) {
 		/* channel ID */
-		int i;
-		float cv,df;
-
-		for(i=1,cv=0;i<7;i++)
-			cv+=tele[n+15*8+i];
-		cv/=6;
-		
-		for(i=0,max=10000.0,channel=-1;i<6;i++) {
-			df=cv-frm[i]; df=df*df;
-			if (df<max) { channel=i; max=df; }
+		for(j=0,max=10000.0,channel=-1;j<6;j++) {
+			float df;
+			df=wedge[15]-wedge[j]; df=df*df;
+			if (df<max) { channel=j; max=df; }
 		}
+		for(j=0;j<16;j++)
+			tele[j]=rgcal(wedge[j],&(regr[k]));
 	}
 }
-mk=k;
+nbtele=k;
 
+/* calibrate */
 for(n=0;n<nrow;n++) {
 	float *pixelv;
 	int	i,c;
@@ -148,27 +149,134 @@ for(n=0;n<nrow;n++) {
 		int k,kof;
 
 		pv=pixelv[i+offset];
-		k=(n-shift)/128;
-		kof=(n-shift)%128;
+
+		k=(n-telestart)/128;
+		kof=(n-telestart)%128;
 		if(kof<64) {
-			if(k<1)
-		 		pv=rgcal(pv,&(regr[k]));
-			else
-		 		pv=rgcal(pv,&(regr[k]))*(64+kof)/128.0+
-		 		   rgcal(pv,&(regr[k-1]))*(64-kof)/128.0;
+			if(k<1) {
+				pv=rgcal(pv,&(regr[k]));
+			} else {
+	 			pv=rgcal(pv,&(regr[k]))*(64+kof)/128.0+
+	 		   	rgcal(pv,&(regr[k-1]))*(64-kof)/128.0;
+			}
 		} else {
-			if((k+1)>=mk)
-		 		pv=rgcal(pv,&(regr[k]));
-			else
-		 		pv=rgcal(pv,&(regr[k]))*(192-kof)/128.0+
-		 		   rgcal(pv,&(regr[k+1]))*(kof-64)/128.0;
+			if((k+1)>=nbtele) {
+				pv=rgcal(pv,&(regr[k]));
+			} else {
+	 			pv=rgcal(pv,&(regr[k]))*(192-kof)/128.0+
+	 		   	rgcal(pv,&(regr[k+1]))*(kof-64)/128.0;
+			}
 		}
 
-		if(pv>1.0) pv=1.0;
+		if(pv>255.0) pv=255.0;
+		if(pv<0.0) pv=0.0;
+		pixelv[i+offset]=pv;
+	}
+}
+printf("Done\n");
+return(channel);
+}
+
+/* ------------------------------temperature calibration -----------------------*/
+extern int satnum;
+const struct {
+        float d[4][3];
+        struct {
+                float vc,A,B;
+        } rad[3];
+        struct {
+                float Ns;
+                float b[3];
+        } cor[3];
+} satcal[2]=
+#include "satcal.h"
+
+typedef struct {
+double  Nbb;
+double  Cs;
+double  Cb;
+int ch;
+} tempparam;
+
+/* temperature compensation for IR channel */
+static void tempcomp(double t[16], int ch,tempparam *tpr)
+{
+double Tbb,T[4];
+double C;
+double r;
+int n;
+
+	tpr->ch=ch-3;
+
+	/* compute equivalent T black body  */
+	for (n=0;n<4;n++) {
+		float d0,d1,d2;
+
+		C=t[9+n]*4;
+		d0=satcal[satnum].d[n][0];
+		d1=satcal[satnum].d[n][1];
+		d2=satcal[satnum].d[n][2];
+		T[n]=d0;
+		T[n]+=d1*C; C=C*C; 
+		T[n]+=d2*C;
+	}
+	Tbb=(T[0]+T[1]+T[2]+T[3])/4.0;
+	Tbb=satcal[satnum].rad[tpr->ch].A+satcal[satnum].rad[tpr->ch].B*Tbb;
+
+	/* compute radiance Black body */
+	C=satcal[satnum].rad[tpr->ch].vc;
+	tpr->Nbb=c1*C*C*C/(exp(c2*C/Tbb)-1.0);
+	/* store Count Blackbody and space */
+	tpr->Cs=1023; /* don't know how to get it in the APT telemetry any idea ? */
+	tpr->Cb=t[14]*4.0;
+}
+
+static double tempcal(float Ce,tempparam *rgpr)
+{
+double Nl,Nc,Ns,Ne;
+double T,vc;
+
+	Ns=satcal[satnum].cor[rgpr->ch].Ns;
+	Nl=Ns+(rgpr->Nbb-Ns)*(rgpr->Cs-Ce*4.0)/(rgpr->Cs-rgpr->Cb);
+	Nc=satcal[satnum].cor[rgpr->ch].b[0]+
+	   satcal[satnum].cor[rgpr->ch].b[1]*Nl+
+	   satcal[satnum].cor[rgpr->ch].b[2]*Nl*Nl;
+
+	Ne=Nl+Nc;
+
+	vc=satcal[satnum].rad[rgpr->ch].vc;
+	T=c2*vc/log(c1*vc*vc*vc/Ne+1.0);
+	T=(T-satcal[satnum].rad[rgpr->ch].A)/satcal[satnum].rad[rgpr->ch].B;
+
+	/* rescale to range 0-255 for -60-+40 °C */
+	T=(T-273.15+60)/100*255.0;
+
+	return(T);
+}
+int Temperature(float **prow,int nrow,int channel,int offset)
+{
+tempparam temp;
+int n;
+
+printf("Temperature ");
+
+	tempcomp(tele,channel,&temp);
+
+
+for(n=0;n<nrow;n++) {
+	float *pixelv;
+	int	i,c;
+
+	pixelv=prow[n];
+	for(i=0;i<954;i++) {
+		float pv;
+
+		pv=tempcal(pixelv[i+offset],&temp);
+
+		if(pv>255.0) pv=255.0;
 		if(pv<0.0) pv=0.0;
 		pixelv[i+offset]=pv;
 	}
 }
 printf("Done\n");
-return (channel);
 }

satcal.h

@@ -0,0 +1,42 @@
+{/* calibration coeff from NOAA KLM POES satellite user guide */
+{/* NOAA-15 */
+{ /* PRT coeff d0,d1,d2 */
+{276.60157 , 0.051045 , 1.36328E-06},
+{276.62531 , 0.050909 , 1.47266E-06},
+{276.67413 , 0.050907 , 1.47656E-06},
+{276.59258 , 0.050966 , 1.47656E-06}
+},
+{ /* channel radiance coeff vc,A,B*/
+{925.4075 , 0.337810 , 0.998719}, /* channel 4 */
+{839.8979 , 0.304558 , 0.999024}, /* channel 5 */
+{2695.9743 , 1.621256 , 0.998015} /* channel 3B */
+},
+{ /*  nonlinear radiance correction Ns,b0,b1,b2 */
+{-4.50 , 4.76 , -0.0932 , 0.0004524}, /* channel 4 */
+{-3.61 , 3.83 , -0.0659 , 0.0002811},  /* channel 5*/
+{0.0,0.0,0.0,0.0}  /* channel 3B*/
+}
+},
+{/* NOAA 17 */
+{ /* PRT coeff d0,d1,d2 */
+{276.628 , 0.05098 , 1.371e-06},
+{276.538 , 0.05098 , 1.371e-06},
+{276.761 , 0.05097 , 1.369e-06},
+{276.660 , 0.05100 , 1.348e-06}
+},
+{ /* channel radiance coeff vc,A,B*/
+{926.2947 , 0.271683 , 0.998794}, /* channel 4 */
+{839.8246 , 0.309180 , 0.999012}, /* channel 5 */
+{2669.3554 , 1.702380 , 0.997378} /* channel 3B */
+},
+{ /*  nonlinear radiance correction Ns,b0,b1,b2 */
+{-8.55 , 8.22 , -0.15795 , 0.00075579}, /* channel 4 */
+{-3.97 , 4.31 , -0.07318 , 0.00030976}, /* channel 5 */
+}
+}
+};
+
+
+const float c1=1.1910427e-5;
+const float c2=1.4387752;
+