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.
 
 
 
 
 

331 line
8.3 KiB

  1. /*
  2. * This file is part of Aptdec.
  3. * Copyright (c) 2004-2009 Thierry Leconte (F4DWV), Xerbo (xerbo@protonmail.com) 2019
  4. *
  5. * Aptdec is free software: you can redistribute it and/or modify
  6. * it under the terms of the GNU General Public License as published by
  7. * the Free Software Foundation, either version 2 of the License, or
  8. * (at your option) any later version.
  9. *
  10. * This program is distributed in the hope that it will be useful,
  11. * but WITHOUT ANY WARRANTY; without even the implied warranty of
  12. * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  13. * GNU General Public License for more details.
  14. *
  15. * You should have received a copy of the GNU General Public License
  16. * along with this program. If not, see <https://www.gnu.org/licenses/>.
  17. *
  18. */
  19. #include <stdio.h>
  20. #include <string.h>
  21. #include <sndfile.h>
  22. #include <math.h>
  23. #include "offsets.h"
  24. #include "messages.h"
  25. #define REGORDER 3
  26. typedef struct {
  27. double cf[REGORDER + 1];
  28. } rgparam;
  29. extern void polyreg(const int m, const int n, const double x[], const double y[], double c[]);
  30. // Compute regression
  31. static void rgcomp(double x[16], rgparam * rgpr) {
  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. polyreg(REGORDER, 9, x, y, rgpr -> cf);
  35. }
  36. // Convert a value to 0-255 based off the provided regression curve
  37. static double rgcal(float x, rgparam *rgpr) {
  38. double y, p;
  39. int i;
  40. for (i = 0, y = 0.0, p = 1.0; i < REGORDER + 1; i++) {
  41. y += rgpr->cf[i] * p;
  42. p = p * x;
  43. }
  44. return(y);
  45. }
  46. static double tele[16];
  47. static double Cs;
  48. static int nbtele;
  49. void histogramEqualise(float **prow, int nrow, int offset, int width){
  50. // Plot histogram
  51. int histogram[256] = { 0 };
  52. for(int y = 0; y < nrow; y++)
  53. for(int x = 0; x < width; x++)
  54. histogram[(int)floor(prow[y][x+offset])]++;
  55. // Find min/max points
  56. int min = -1, max = -1;
  57. for(int i = 5; i < 250; i++){
  58. if(histogram[i]/width/(nrow/255.0) > 1.0){
  59. if(min == -1) min = i;
  60. max = i;
  61. }
  62. }
  63. //printf("Min Value: %i, Max Value %i\n", min, max);
  64. // Spread values to avoid overshoot
  65. min -= 5; max += 5;
  66. // Stretch the brightness into the new range
  67. for(int y = 0; y < nrow; y++)
  68. for(int x = 0; x < width; x++)
  69. prow[y][x+offset] = (prow[y][x+offset]-min) / (max-min) * 255;
  70. }
  71. // Brightness calibrate, including telemetry
  72. void calibrateBrightness(float **prow, int nrow, int offset, int width, int telestart, rgparam regr[30]){
  73. offset -= SYNC_WIDTH+SPC_WIDTH;
  74. for (int n = 0; n < nrow; n++) {
  75. float *pixelv = prow[n];
  76. for (int i = 0; i < width+SYNC_WIDTH+SPC_WIDTH+TELE_WIDTH; i++) {
  77. float pv = pixelv[i + offset];
  78. // Blend between the calculated regression curves
  79. /* TODO: this can actually make the image look *worse*
  80. * if the signal has a constant input gain.
  81. */
  82. int k, kof;
  83. k = (n - telestart) / FRAME_LEN;
  84. if (k >= nbtele)
  85. k = nbtele - 1;
  86. kof = (n - telestart) % FRAME_LEN;
  87. if (kof < 64) {
  88. if (k < 1) {
  89. pv = rgcal(pv, &(regr[k]));
  90. } else {
  91. pv = rgcal(pv, &(regr[k])) * (64 + kof) / FRAME_LEN +
  92. rgcal(pv, &(regr[k - 1])) * (64 - kof) / FRAME_LEN;
  93. }
  94. } else {
  95. if ((k + 1) >= nbtele) {
  96. pv = rgcal(pv, &(regr[k]));
  97. } else {
  98. pv = rgcal(pv, &(regr[k])) * (192 - kof) / FRAME_LEN +
  99. rgcal(pv, &(regr[k + 1])) * (kof - 64) / FRAME_LEN;
  100. }
  101. }
  102. pv = CLIP(pv, 0, 255);
  103. pixelv[i + offset] = pv;
  104. }
  105. }
  106. }
  107. // Get telemetry data for thermal calibration/equalization
  108. int calibrate(float **prow, int nrow, int offset, int width) {
  109. double teleline[3000] = { 0.0 };
  110. double wedge[16];
  111. rgparam regr[30];
  112. int n, k;
  113. int mtelestart = 0, telestart;
  114. int channel = -1;
  115. // Calculate average of a row of telemetry
  116. for (n = 0; n < nrow; n++) {
  117. float *pixelv = prow[n];
  118. // Average the center 40px
  119. for (int i = 3; i < 43; i++) teleline[n] += pixelv[i + offset + width];
  120. teleline[n] /= 40.0;
  121. }
  122. // The minimum rows required to decode a full frame
  123. if (nrow < 192) {
  124. fprintf(stderr, ERR_TELE_ROW);
  125. return(0);
  126. }
  127. /* Wedge 7 is white and 8 is black, which will have the largest
  128. * difference in brightness, this will always be in the center of
  129. * the frame and can thus be used to find the start of the frame
  130. */
  131. double max = 0.0;
  132. for (n = nrow / 3 - 64; n < 2 * nrow / 3 - 64; n++) {
  133. float df;
  134. // (sum 4px below) / (sum 4px above)
  135. df = (teleline[n - 4] + teleline[n - 3] + teleline[n - 2] + teleline[n - 1]) /
  136. (teleline[n + 0] + teleline[n + 1] + teleline[n + 2] + teleline[n + 3]);
  137. // Find the maximum difference
  138. if (df > max) {
  139. mtelestart = n;
  140. max = df;
  141. }
  142. }
  143. // Find the start of the first frame
  144. telestart = (mtelestart - FRAME_LEN/2) % FRAME_LEN;
  145. // Make sure that theres at least one full frame in the image
  146. if (nrow < telestart + FRAME_LEN) {
  147. fprintf(stderr, ERR_TELE_ROW);
  148. return(0);
  149. }
  150. // For each frame
  151. for (n = telestart, k = 0; n < nrow - FRAME_LEN; n += FRAME_LEN, k++) {
  152. float *pixelv = prow[n];
  153. int j;
  154. // Turn each wedge into a value
  155. for (j = 0; j < 16; j++) {
  156. // Average the middle 6px
  157. wedge[j] = 0.0;
  158. for (int i = 1; i < 7; i++) wedge[j] += teleline[(j * 8) + n + i];
  159. wedge[j] /= 6;
  160. }
  161. // Compute regression on the wedges
  162. rgcomp(wedge, &(regr[k]));
  163. // Read the telemetry values from the middle of the image
  164. if (k == nrow / (2*FRAME_LEN)) {
  165. int l;
  166. // Equalise
  167. for (j = 0; j < 16; j++) tele[j] = rgcal(wedge[j], &(regr[k]));
  168. /* Compare the channel ID wedge to the reference
  169. * wedges, the wedge with the closest match will
  170. * be the channel ID
  171. */
  172. float min = 10000;
  173. for (j = 0; j < 6; j++) {
  174. float df;
  175. df = tele[15] - tele[j];
  176. df *= df;
  177. if (df < min) {
  178. channel = j;
  179. min = df;
  180. }
  181. }
  182. // Cs computation, still have no idea what this does
  183. int i;
  184. for (Cs = 0.0, i = 0, j = n; j < n + FRAME_LEN; j++) {
  185. double csline;
  186. for (csline = 0.0, l = 3; l < 43; l++)
  187. csline += pixelv[l + offset - SPC_WIDTH];
  188. csline /= 40.0;
  189. if (csline > 50.0) {
  190. Cs += csline;
  191. i++;
  192. }
  193. }
  194. Cs /= i;
  195. Cs = rgcal(Cs, &(regr[k]));
  196. }
  197. }
  198. nbtele = k;
  199. calibrateBrightness(prow, nrow, offset, width, telestart, regr);
  200. return(channel + 1);
  201. }
  202. // --- Temperature Calibration --- //
  203. extern int satnum;
  204. #include "satcal.h"
  205. typedef struct {
  206. double Nbb;
  207. double Cs;
  208. double Cb;
  209. int ch;
  210. } tempparam;
  211. // IR channel temperature compensation
  212. static void tempcomp(double t[16], int ch, tempparam *tpr) {
  213. double Tbb, T[4];
  214. double C;
  215. tpr -> ch = ch - 4;
  216. // Compute equivalent T blackbody temperature
  217. for (int n = 0; n < 4; n++) {
  218. float d0, d1, d2;
  219. C = t[9 + n] * 4.0;
  220. d0 = satcal[satnum].d[n][0];
  221. d1 = satcal[satnum].d[n][1];
  222. d2 = satcal[satnum].d[n][2];
  223. T[n] = d0;
  224. T[n] += d1 * C;
  225. C = C * C;
  226. T[n] += d2 * C;
  227. }
  228. Tbb = (T[0] + T[1] + T[2] + T[3]) / 4.0;
  229. Tbb = satcal[satnum].rad[tpr->ch].A + satcal[satnum].rad[tpr->ch].B * Tbb;
  230. // Compute radiance blackbody
  231. C = satcal[satnum].rad[tpr->ch].vc;
  232. tpr->Nbb = c1 * C * C * C / (expm1(c2 * C / Tbb));
  233. // Store count blackbody and space
  234. tpr->Cs = Cs * 4.0;
  235. tpr->Cb = t[14] * 4.0;
  236. }
  237. // IR channel temperature calibration
  238. static double tempcal(float Ce, tempparam * rgpr) {
  239. double Nl, Nc, Ns, Ne;
  240. double T, vc;
  241. Ns = satcal[satnum].cor[rgpr->ch].Ns;
  242. Nl = Ns + (rgpr->Nbb - Ns) * (rgpr->Cs - Ce * 4.0) / (rgpr->Cs - rgpr->Cb);
  243. Nc = satcal[satnum].cor[rgpr->ch].b[0] +
  244. satcal[satnum].cor[rgpr->ch].b[1] * Nl +
  245. satcal[satnum].cor[rgpr->ch].b[2] * Nl * Nl;
  246. Ne = Nl + Nc;
  247. vc = satcal[satnum].rad[rgpr->ch].vc;
  248. T = c2 * vc / log1p(c1 * vc * vc * vc / Ne);
  249. T = (T - satcal[satnum].rad[rgpr->ch].A) / satcal[satnum].rad[rgpr->ch].B;
  250. // Rescale to 0-255 for -60'C to +40'C
  251. T = (T - 273.15 + 60.0) / 100.0 * 256.0;
  252. return(T);
  253. }
  254. // Temperature calibration wrapper
  255. void temperature(float **prow, int nrow, int ch, int offset){
  256. tempparam temp;
  257. printf("Temperature... ");
  258. fflush(stdout);
  259. tempcomp(tele, ch, &temp);
  260. for (int n = 0; n < nrow; n++) {
  261. float *pixelv = prow[n];
  262. for (int i = 0; i < CH_WIDTH; i++) {
  263. float pv = tempcal(pixelv[i + offset], &temp);
  264. pv = CLIP(pv, 0, 255);
  265. pixelv[i + offset] = pv;
  266. }
  267. }
  268. printf("Done\n");
  269. }