您最多选择25个主题 主题必须以字母或数字开头,可以包含连字符 (-),并且长度不得超过35个字符
 
 
 
 
 

297 行
5.7 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.434,0.542,0.652,0.78,0.87 ,0.0 }; */
  33. const double y[9] = { 31.1,63.0,95.0,127.2,158.9,191.1,228.6,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 csline[3000];
  54. double wedge[16];
  55. rgparam regr[30];
  56. int n;
  57. int k;
  58. int mtelestart,telestart;
  59. int channel=-1;
  60. float max;
  61. printf("Calibration ");
  62. fflush(stdout);
  63. /* build telemetry values lines */
  64. for(n=0;n<nrow;n++) {
  65. int i;
  66. teleline[n]=0.0;
  67. csline[n]=0.0;
  68. for(i=3;i<43;i++) {
  69. teleline[n]+=prow[n][i+offset+909];
  70. csline[n]+=prow[n][i+offset-50];
  71. }
  72. teleline[n]/=40;
  73. csline[n]/=40;
  74. }
  75. if(nrow<192) {
  76. fprintf(stderr," impossible, not enought row\n");
  77. return (0);
  78. }
  79. /* find telemetry start */
  80. max=0.0;mtelestart=0;
  81. for(n=nrow/3-64;n<2*nrow/3-64;n++) {
  82. float df;
  83. df=(teleline[n-4]+teleline[n-3]+teleline[n-2]+teleline[n-1])/(teleline[n]+teleline[n+1]+teleline[n+2]+teleline[n+3]);
  84. if(df> max){ mtelestart=n; max=df; }
  85. }
  86. mtelestart-=64;
  87. telestart=mtelestart%128;
  88. if(mtelestart <0 || nrow<telestart+128) {
  89. fprintf(stderr," impossible, not enought row\n");
  90. return (0);
  91. }
  92. /* compute wedges and regression */
  93. for(n=telestart,k=0;n<nrow-128;n+=128,k++) {
  94. int j;
  95. for(j=0;j<16;j++) {
  96. int i;
  97. wedge[j]=0.0;
  98. for(i=1;i<7;i++)
  99. wedge[j]+=teleline[n+j*8+i];
  100. wedge[j]/=6;
  101. }
  102. rgcomp(wedge,&(regr[k]));
  103. if (k==1) {
  104. int i;
  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. /* telemetry computation */
  112. for(j=0;j<16;j++)
  113. tele[j]=rgcal(wedge[j],&(regr[k]));
  114. /* Cs computation */
  115. Cs=0.0;i=0;
  116. for(j=n;j<n+128;j++)
  117. if(csline[j]>50.0) {
  118. Cs+=csline[j];
  119. i++;
  120. }
  121. Cs/=128;
  122. Cs=rgcal(Cs,&(regr[k]));
  123. }
  124. }
  125. nbtele=k;
  126. /* calibrate */
  127. for(n=0;n<nrow;n++) {
  128. float *pixelv;
  129. int i,c;
  130. pixelv=prow[n];
  131. for(i=0;i<954;i++) {
  132. float pv;
  133. int k,kof;
  134. pv=pixelv[i+offset];
  135. k=(n-telestart)/128;
  136. if(k>=nbtele) k=nbtele-1;
  137. kof=(n-telestart)%128;
  138. if(kof<64) {
  139. if(k<1) {
  140. pv=rgcal(pv,&(regr[k]));
  141. } else {
  142. pv=rgcal(pv,&(regr[k]))*(64+kof)/128.0+
  143. rgcal(pv,&(regr[k-1]))*(64-kof)/128.0;
  144. }
  145. } else {
  146. if((k+1)>=nbtele) {
  147. pv=rgcal(pv,&(regr[k]));
  148. } else {
  149. pv=rgcal(pv,&(regr[k]))*(192-kof)/128.0+
  150. rgcal(pv,&(regr[k+1]))*(kof-64)/128.0;
  151. }
  152. }
  153. if(pv>255.0) pv=255.0;
  154. if(pv<0.0) pv=0.0;
  155. pixelv[i+offset]=pv;
  156. }
  157. }
  158. printf("Done\n");
  159. return(channel);
  160. }
  161. /* ------------------------------temperature calibration -----------------------*/
  162. extern int satnum;
  163. const struct {
  164. float d[4][3];
  165. struct {
  166. float vc,A,B;
  167. } rad[3];
  168. struct {
  169. float Ns;
  170. float b[3];
  171. } cor[3];
  172. } satcal[2]=
  173. #include "satcal.h"
  174. typedef struct {
  175. double Nbb;
  176. double Cs;
  177. double Cb;
  178. int ch;
  179. } tempparam;
  180. /* temperature compensation for IR channel */
  181. static void tempcomp(double t[16], int ch,tempparam *tpr)
  182. {
  183. double Tbb,T[4];
  184. double C;
  185. double r;
  186. int n;
  187. tpr->ch=ch-3;
  188. /* compute equivalent T black body */
  189. for (n=0;n<4;n++) {
  190. float d0,d1,d2;
  191. C=t[9+n]*4.0;
  192. d0=satcal[satnum].d[n][0];
  193. d1=satcal[satnum].d[n][1];
  194. d2=satcal[satnum].d[n][2];
  195. T[n]=d0;
  196. T[n]+=d1*C; C=C*C;
  197. T[n]+=d2*C;
  198. }
  199. Tbb=(T[0]+T[1]+T[2]+T[3])/4.0;
  200. Tbb=satcal[satnum].rad[tpr->ch].A+satcal[satnum].rad[tpr->ch].B*Tbb;
  201. /* compute radiance Black body */
  202. C=satcal[satnum].rad[tpr->ch].vc;
  203. tpr->Nbb=c1*C*C*C/(exp(c2*C/Tbb)-1.0);
  204. /* store Count Blackbody and space */
  205. tpr->Cs=Cs*4.0;
  206. tpr->Cb=t[14]*4.0;
  207. }
  208. static double tempcal(float Ce,tempparam *rgpr)
  209. {
  210. double Nl,Nc,Ns,Ne;
  211. double T,vc;
  212. Ns=satcal[satnum].cor[rgpr->ch].Ns;
  213. Nl=Ns+(rgpr->Nbb-Ns)*(rgpr->Cs-Ce*4.0)/(rgpr->Cs-rgpr->Cb);
  214. Nc=satcal[satnum].cor[rgpr->ch].b[0]+
  215. satcal[satnum].cor[rgpr->ch].b[1]*Nl+
  216. satcal[satnum].cor[rgpr->ch].b[2]*Nl*Nl;
  217. Ne=Nl+Nc;
  218. vc=satcal[satnum].rad[rgpr->ch].vc;
  219. T=c2*vc/log(c1*vc*vc*vc/Ne+1.0);
  220. T=(T-satcal[satnum].rad[rgpr->ch].A)/satcal[satnum].rad[rgpr->ch].B;
  221. /* rescale to range 0-255 for -60-+40 °C */
  222. T=(T-273.15+60)/100*255.0;
  223. return(T);
  224. }
  225. int Temperature(float **prow,int nrow,int channel,int offset)
  226. {
  227. tempparam temp;
  228. int n;
  229. printf("Temperature ");
  230. tempcomp(tele,channel,&temp);
  231. for(n=0;n<nrow;n++) {
  232. float *pixelv;
  233. int i,c;
  234. pixelv=prow[n];
  235. for(i=0;i<954;i++) {
  236. float pv;
  237. pv=tempcal(pixelv[i+offset],&temp);
  238. if(pv>255.0) pv=255.0;
  239. if(pv<0.0) pv=0.0;
  240. pixelv[i+offset]=pv;
  241. }
  242. }
  243. printf("Done\n");
  244. }