From c4c771ad4c3eac7c08c6f4d24135ead677b87506 Mon Sep 17 00:00:00 2001
From: Thierry Leconte <f4dwv@users.sourceforge.net>
Date: Sat, 11 Oct 2003 15:09:36 +0000
Subject: [PATCH] change calibration : use linear regression

---
 image.c | 150 ++++++++++++++++++++++++++++++++++++++++----------------
 main.c  |  76 +++++++++++++++++++++++-----
 2 files changed, 170 insertions(+), 56 deletions(-)

diff --git a/image.c b/image.c
index 5a9c95a..ca050fa 100644
--- a/image.c
+++ b/image.c
@@ -24,84 +24,148 @@
 #include <string.h>
 #include <sndfile.h>
 
+
+typedef struct {
+double slope;
+double offset;
+} rgparam;
+
+static void rgcomp(double *x, 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;
+double xavg;
+double sxx,sxy;
+int i;
+
+for(i=0,xavg=0.0;i<9;i++)
+	xavg+=x[i];
+xavg/=9;
+for(i=0,sxx=0.0;i<9;i++) {
+	float t=x[i]-xavg;
+	sxx+=t*t;
+}
+for(i=0,sxy=0.0;i<9;i++) {
+	sxy+=(x[i]-xavg)*(y[i]-yavg);
+}
+rgpr->slope=sxy/sxx;
+rgpr->offset=yavg-rgpr->slope*xavg;
+
+}
+
+static double rgcal(float x,rgparam *rgpr)
+{
+return(rgpr->slope*x+rgpr->offset);
+}
+
 int Calibrate(float **prow,int nrow,int offset)
 {
-/* const double cal[8] = { 0.106,0.215,0.324,0.434,0.542,0.652,0.78,0.87 }; */
-const double cal[8] = { 0.1155 ,0.234,0.353,0.473,0.59,0.71,0.85,0.95 };
-
-float tele[2000];
-float frm[16];
-float slv[8];
-int n,i;
-int shift,channel;
+
+double tele[3000];
+double frm[9];
+rgparam regr[30];
+int n;
+int k,mk;
+int shift;
+int mshift,channel;
 float max;
 
 
 printf("Calibration ");
 fflush(stdout);
 
-if(nrow<128) {
-	fprintf(stderr," impossible, not enought row\n");
-	return (0);
-}
-
-/* build telemetry values */
+/* build telemetry values lines */
 for(n=0;n<nrow;n++) {
+	int i;
+
 	tele[n]=0.0;
 	for(i=3;i<43;i++)
 		tele[n]+=prow[n][i+offset+909];
 	tele[n]/=40;
 }
 
-max=0.0;shift=0;
-for(n=60;n<nrow-64;n++) {
+if(nrow<128) {
+	fprintf(stderr," impossible, not enought row\n");
+	return (0);
+}
+
+/* find telemetry start */
+max=-1000.0;mshift=0;
+for(n=60;n<nrow-3;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){ shift=n; max=df; }
+	if(df> max){ mshift=n; max=df; }
 }
-for(n=0;n<16;n++) {
-	frm[n]=0.0;
-	for(i=1;i<7;i++)
-		frm[n]+=tele[shift-64+n*8+i];
-	frm[n]/=6;
+
+mshift-=64;
+shift=mshift%128;
+if(nrow<shift+8*9) {
+	fprintf(stderr," impossible, not enought row\n");
+	return (0);
 }
 
-/* channel I.D */
-max=100000.0;
-for(n=0;n<6;n++) {
-	float df;
-	df=fabs(frm[n]-frm[15]);
-	if(df<max) { 
-		max=df;
-		channel=n;
+/* compute regression params */
+for(n=shift,k=0;n<nrow-8*9;n+=128,k++) {
+	int j;
+
+	for(j=0;j<9;j++) {
+		int i;
+
+		frm[j]=0.0;
+		for(i=1;i<7;i++)
+			frm[j]+=tele[n+j*8+i];
+		frm[j]/=6;
 	}
-}
-fflush(stdout);
+	rgcomp(frm,&(regr[k]));
+
 
-/* compute slopes */
-slv[0]=cal[0]/(frm[0]-frm[8]);
-for(n=1;n<8;n++) {
-	slv[n]=(cal[n]-cal[n-1])/(frm[n]-frm[n-1]);
+	if(n==mshift) {
+		/* 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; }
+		}
+	}
 }
+mk=k;
 
 for(n=0;n<nrow;n++) {
 	float *pixelv;
 	int	i,c;
 
 	pixelv=prow[n];
-	for(i=0;i<909;i++) {
+	for(i=0;i<954;i++) {
 		float pv;
+		int k,kof;
 
 		pv=pixelv[i+offset];
+		k=(n-shift)/128;
+		kof=(n-shift)%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;
+		} 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(pv<frm[8]) {pv=0.0; goto conv; }
-		if(pv<frm[0]) {pv=(pv-frm[8])*slv[0]; goto conv; }
-		for(c=1;c<8;c++) 
-		    if(pv<frm[c]) {pv=(pv-frm[c-1])*slv[c]+cal[c-1]; goto conv; }
- 	        pv=(pv-frm[7])*slv[7]+cal[7];
 		if(pv>1.0) pv=1.0;
-conv:
+		if(pv<0.0) pv=0.0;
 		pixelv[i+offset]=pv;
 	}
 }
diff --git a/main.c b/main.c
index 055fb8f..6c5efe1 100644
--- a/main.c
+++ b/main.c
@@ -163,9 +163,9 @@ for(n=0;n<nrow;n++) {
 	for(i=0;i<909;i++) {
 		int x,y;
 
-		x=(int)(pixelv[i+85]*255.0);
-		y=(int)(pixelv[i+1125]*255.0);
-		pixel[i]=cmap[x][y];
+		y=(int)(pixelv[i+85]*255.0);
+		x=(int)(pixelv[i+1125]*255.0);
+		pixel[i]=cmap[y][x];
 	}
 	png_write_row(png_ptr,(png_bytep)pixel);
 }
@@ -232,6 +232,47 @@ png_destroy_read_struct(&png_ptr,&info_ptr,NULL);
 return(0);
 }
 
+#ifdef DEBUG
+unsigned int distrib[256][256];
+int Distrib(char *filename,float **prow,int nrow)
+{
+int n;
+int x,y;
+int max=0;
+FILE *df;
+
+for(y=0;y<256;y++)
+	for(x=0;x<256;x++)
+                distrib[y][x]=0;
+
+for(n=0;n<nrow;n++) {
+        float *pixelv;
+        png_color pixel[909];
+        int     i;
+
+        pixelv=prow[n];
+        for(i=0;i<909;i++) {
+
+                y=(int)(pixelv[i+85]*255.0);
+                x=(int)(pixelv[i+1125]*255.0);
+		if(x>255) x=255;
+		if(x<0) x=0;
+		if(y>255) y=255;
+		if(y<0) y=0;
+                distrib[y][x]+=1;
+		if(distrib[y][x]> max) max=distrib[y][x];
+        }
+}
+df=fopen(filename,"w");
+fprintf(df,"P2\n#max %d\n",max);
+fprintf(df,"256 256\n255\n");
+for(y=0;y<256;y++)
+	for(x=0;x<256;x++)
+		fprintf(df,"%d\n",(int)((255.0*(double)(distrib[y][x]))/(double)max));
+fclose(df);
+}
+#endif
+
 extern int Calibrate(float **prow,int nrow,int offset);
 extern int optind,opterr;
 extern char *optarg;
@@ -299,25 +340,34 @@ ImageOut(pngfilename,prow,nrow,2080,0);
 }
 
 /* Channel A */
-ch=Calibrate(prow,nrow,85);
-if(ch>0) {
-	a=1;
-	if(((strchr(imgopt,(int)'a')!=NULL) || (strchr(imgopt,ch)!=NULL))) {
+if(((strchr(imgopt,(int)'a')!=NULL) || (strchr(imgopt,(int)'c')!=NULL))) {
+	ch=Calibrate(prow,nrow,85);
+	if(ch>0) {
+		a=1;
+		if(strchr(imgopt,(int)'a')!=NULL) {
 		sprintf(pngfilename,"%s/%s-%s.png",pngdirname,name,chid[ch]);
-	ImageOut(pngfilename,prow,nrow,909,85);
+		ImageOut(pngfilename,prow,nrow,954,85);
+		}
 	}
 }
 
 /* Channel B */
-ch=Calibrate(prow,nrow,1125);
-if(ch>0) {
-	b=1;
-	if(((strchr(imgopt,(int)'b')!=NULL) || (strchr(imgopt,ch)!=NULL))) {
+if(((strchr(imgopt,(int)'b')!=NULL) || (strchr(imgopt,(int)'c')!=NULL))) {
+	ch=Calibrate(prow,nrow,1125);
+	if(ch>0) {
+		b=1;
+		if(strchr(imgopt,(int)'b')!=NULL) {
 		sprintf(pngfilename,"%s/%s-%s.png",pngdirname,name,chid[ch]);
-		ImageOut(pngfilename,prow,nrow,909,1125);
+		ImageOut(pngfilename,prow,nrow,954,1125);
+		}
 	}
 }
 
+#ifdef DEBUG
+	sprintf(pngfilename,"%s/%s-d.pnm",pngdirname,name);
+	Distrib(pngfilename,prow,nrow);
+#endif
+
 /* color image */
 if(a && b && strchr(imgopt,(int)'c')!=NULL){
 	sprintf(pngfilename,"%s/%s-c.png",pngdirname,name);