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.
 
 
 
 
 

324 lines
6.7 KiB

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