You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 
 
 
 

301 regels
5.8 KiB

  1. /*
  2. * Atpdec
  3. * Copyright (c) 2003 by Thierry Leconte (F4DWV)
  4. *
  5. * $Id$
  6. *
  7. * This library is free software; you can redistribute it and/or modify
  8. * it under the terms of the GNU Library General Public License as
  9. * published by the Free Software Foundation; either version 2 of
  10. * the License, or (at your option) any later version.
  11. *
  12. * This program is distributed in the hope that it will be useful,
  13. * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14. * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  15. * GNU Library General Public License for more details.
  16. *
  17. * You should have received a copy of the GNU Library General Public
  18. * License along with this library; if not, write to the Free Software
  19. * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  20. *
  21. */
  22. #include <stdio.h>
  23. #include <string.h>
  24. #include <sndfile.h>
  25. #include <math.h>
  26. #define REGORDER 3
  27. typedef struct {
  28. double cf[REGORDER+1] ;
  29. } rgparam;
  30. static void rgcomp(double x[16], rgparam *rgpr)
  31. {
  32. /*{ 0.106,0.215,0.324,0.433,0.542,0.652,0.78,0.87 ,0.0 }; */
  33. const double y[9] = { 31.07,63.02,94.96,126.9,158.86,191.1,228.62,255.0, 0.0 };
  34. extern void polyreg(const int m,const int n,const double x[],const double y[],double c[]);
  35. polyreg(REGORDER,9,x,y,rgpr->cf);
  36. }
  37. static double rgcal(float x,rgparam *rgpr)
  38. {
  39. double y,p;
  40. int i;
  41. for(i=0,y=0.0,p=1.0;i<REGORDER+1;i++) {
  42. y+=rgpr->cf[i]*p;
  43. p=p*x;
  44. }
  45. return(y);
  46. }
  47. static double tele[16];
  48. static double Cs;
  49. static nbtele;
  50. int Calibrate(float **prow,int nrow,int offset)
  51. {
  52. double teleline[3000];
  53. double wedge[16];
  54. rgparam regr[30];
  55. int n;
  56. int k;
  57. int mtelestart,telestart;
  58. int channel=-1;
  59. float max;
  60. printf("Calibration ");
  61. fflush(stdout);
  62. /* build telemetry values lines */
  63. for(n=0;n<nrow;n++) {
  64. int i;
  65. teleline[n]=0.0;
  66. for(i=3;i<43;i++) {
  67. teleline[n]+=prow[n][i+offset+956];
  68. }
  69. teleline[n]/=40.0;
  70. }
  71. if(nrow<192) {
  72. fprintf(stderr," impossible, not enought row\n");
  73. return (0);
  74. }
  75. /* find telemetry start in the 2nd third */
  76. max=0.0;mtelestart=0;
  77. for(n=nrow/3-64;n<2*nrow/3-64;n++) {
  78. float df;
  79. df=(teleline[n-4]+teleline[n-3]+teleline[n-2]+teleline[n-1])/(teleline[n]+teleline[n+1]+teleline[n+2]+teleline[n+3]);
  80. if(df> max){ mtelestart=n; max=df; }
  81. }
  82. mtelestart-=64;
  83. telestart=mtelestart%128;
  84. if(mtelestart <0 || nrow<telestart+128) {
  85. fprintf(stderr," impossible, not enought row\n");
  86. return (0);
  87. }
  88. /* compute wedges and regression */
  89. for(n=telestart,k=0;n<nrow-128;n+=128,k++) {
  90. int j;
  91. for(j=0;j<16;j++) {
  92. int i;
  93. wedge[j]=0.0;
  94. for(i=1;i<7;i++)
  95. wedge[j]+=teleline[n+j*8+i];
  96. wedge[j]/=6;
  97. }
  98. rgcomp(wedge,&(regr[k]));
  99. if (k==nrow/256) {
  100. int i,l;
  101. /* telemetry calibration */
  102. for(j=0;j<16;j++) {
  103. tele[j]=rgcal(wedge[j],&(regr[k]));
  104. }
  105. /* channel ID */
  106. for(j=0,max=10000.0,channel=-1;j<6;j++) {
  107. float df;
  108. df=wedge[15]-wedge[j]; df=df*df;
  109. if (df<max) { channel=j; max=df; }
  110. }
  111. /* Cs computation */
  112. for(Cs=0.0,i=0,j=n;j<n+128;j++) {
  113. double csline;
  114. for(csline=0.0,l=3;l<43;l++)
  115. csline+=prow[n][l+offset];
  116. csline/=40.0;
  117. if(csline>50.0) {
  118. Cs+=csline;
  119. i++;
  120. }
  121. }
  122. Cs/=i;
  123. Cs=rgcal(Cs,&(regr[k]));
  124. }
  125. }
  126. nbtele=k;
  127. /* calibrate */
  128. for(n=0;n<nrow;n++) {
  129. float *pixelv;
  130. int i,c;
  131. pixelv=prow[n];
  132. for(i=0;i<1001;i++) {
  133. float pv;
  134. int k,kof;
  135. pv=pixelv[i+offset];
  136. k=(n-telestart)/128;
  137. if(k>=nbtele) k=nbtele-1;
  138. kof=(n-telestart)%128;
  139. if(kof<64) {
  140. if(k<1) {
  141. pv=rgcal(pv,&(regr[k]));
  142. } else {
  143. pv=rgcal(pv,&(regr[k]))*(64+kof)/128.0+
  144. rgcal(pv,&(regr[k-1]))*(64-kof)/128.0;
  145. }
  146. } else {
  147. if((k+1)>=nbtele) {
  148. pv=rgcal(pv,&(regr[k]));
  149. } else {
  150. pv=rgcal(pv,&(regr[k]))*(192-kof)/128.0+
  151. rgcal(pv,&(regr[k+1]))*(kof-64)/128.0;
  152. }
  153. }
  154. if(pv>255.0) pv=255.0;
  155. if(pv<0.0) pv=0.0;
  156. pixelv[i+offset]=pv;
  157. }
  158. }
  159. printf("Done\n");
  160. return(channel);
  161. }
  162. /* ------------------------------temperature calibration -----------------------*/
  163. extern int satnum;
  164. const struct {
  165. float d[4][3];
  166. struct {
  167. float vc,A,B;
  168. } rad[3];
  169. struct {
  170. float Ns;
  171. float b[3];
  172. } cor[3];
  173. } satcal[2]=
  174. #include "satcal.h"
  175. typedef struct {
  176. double Nbb;
  177. double Cs;
  178. double Cb;
  179. int ch;
  180. } tempparam;
  181. /* temperature compensation for IR channel */
  182. static void tempcomp(double t[16], int ch,tempparam *tpr)
  183. {
  184. double Tbb,T[4];
  185. double C;
  186. double r;
  187. int n;
  188. tpr->ch=ch-3;
  189. /* compute equivalent T black body */
  190. for (n=0;n<4;n++) {
  191. float d0,d1,d2;
  192. C=t[9+n]*4.0;
  193. d0=satcal[satnum].d[n][0];
  194. d1=satcal[satnum].d[n][1];
  195. d2=satcal[satnum].d[n][2];
  196. T[n]=d0;
  197. T[n]+=d1*C; C=C*C;
  198. T[n]+=d2*C;
  199. }
  200. Tbb=(T[0]+T[1]+T[2]+T[3])/4.0;
  201. Tbb=satcal[satnum].rad[tpr->ch].A+satcal[satnum].rad[tpr->ch].B*Tbb;
  202. /* compute radiance Black body */
  203. C=satcal[satnum].rad[tpr->ch].vc;
  204. tpr->Nbb=c1*C*C*C/(exp(c2*C/Tbb)-1.0);
  205. /* store Count Blackbody and space */
  206. tpr->Cs=Cs*4.0;
  207. tpr->Cb=t[14]*4.0;
  208. }
  209. static double tempcal(float Ce,tempparam *rgpr)
  210. {
  211. double Nl,Nc,Ns,Ne;
  212. double T,vc;
  213. Ns=satcal[satnum].cor[rgpr->ch].Ns;
  214. Nl=Ns+(rgpr->Nbb-Ns)*(rgpr->Cs-Ce*4.0)/(rgpr->Cs-rgpr->Cb);
  215. Nc=satcal[satnum].cor[rgpr->ch].b[0]+
  216. satcal[satnum].cor[rgpr->ch].b[1]*Nl+
  217. satcal[satnum].cor[rgpr->ch].b[2]*Nl*Nl;
  218. Ne=Nl+Nc;
  219. vc=satcal[satnum].rad[rgpr->ch].vc;
  220. T=c2*vc/log(c1*vc*vc*vc/Ne+1.0);
  221. T=(T-satcal[satnum].rad[rgpr->ch].A)/satcal[satnum].rad[rgpr->ch].B;
  222. /* rescale to range 0-255 for +40-60 °C */
  223. T=(-T+273.15+40)/100.0*255.0;
  224. return(T);
  225. }
  226. int Temperature(float **prow,int nrow,int channel,int offset)
  227. {
  228. tempparam temp;
  229. int n;
  230. printf("Temperature ");
  231. tempcomp(tele,channel,&temp);
  232. for(n=0;n<nrow;n++) {
  233. float *pixelv;
  234. int i,c;
  235. pixelv=prow[n];
  236. for(i=0;i<1001;i++) {
  237. float pv;
  238. pv=tempcal(pixelv[i+offset],&temp);
  239. if(pv>255.0) pv=255.0;
  240. if(pv<0.0) pv=0.0;
  241. pixelv[i+offset]=pv;
  242. }
  243. }
  244. printf("Done\n");
  245. }