選択できるのは25トピックまでです。 トピックは、先頭が英数字で、英数字とダッシュ('-')を使用した35文字以内のものにしてください。
 
 
 
 
 

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