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.
 
 
 
 
 

404 lines
10 KiB

  1. /*
  2. * This file is part of Aptdec.
  3. * Copyright (c) 2004-2009 Thierry Leconte (F4DWV), Xerbo (xerbo@protonmail.com) 2019-2020
  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 <stdlib.h>
  24. #include "common.h"
  25. #include "offsets.h"
  26. #define REGORDER 3
  27. typedef struct {
  28. double cf[REGORDER + 1];
  29. } rgparam_t;
  30. extern void polyreg(const int m, const int n, const double x[], const double y[], double c[]);
  31. // Compute regression
  32. static void rgcomp(double x[16], rgparam_t * rgpr) {
  33. // { 0.106, 0.215, 0.324, 0.433, 0.542, 0.652, 0.78, 0.87, 0.0 }
  34. const double y[9] = { 31.07, 63.02, 94.96, 126.9, 158.86, 191.1, 228.62, 255.0, 0.0 };
  35. polyreg(REGORDER, 9, x, y, rgpr->cf);
  36. }
  37. // Convert a value to 0-255 based off the provided regression curve
  38. static double rgcal(float x, rgparam_t *rgpr) {
  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. 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)CLIP(prow[y][x+offset], 0, 255)]++;
  55. // Calculate cumulative frequency
  56. long sum = 0, cf[256] = { 0 };
  57. for(int i = 0; i < 255; i++){
  58. sum += histogram[i];
  59. cf[i] = sum;
  60. }
  61. // Apply histogram
  62. int area = nrow * width;
  63. for(int y = 0; y < nrow; y++){
  64. for(int x = 0; x < width; x++){
  65. int k = prow[y][x+offset];
  66. prow[y][x+offset] = (256.0/area) * cf[k];
  67. }
  68. }
  69. }
  70. void linearEnhance(float **prow, int nrow, int offset, int width){
  71. // Plot histogram
  72. int histogram[256] = { 0 };
  73. for(int y = 0; y < nrow; y++)
  74. for(int x = 0; x < width; x++)
  75. histogram[(int)CLIP(prow[y][x+offset], 0, 255)]++;
  76. // Find min/max points
  77. int min = -1, max = -1;
  78. for(int i = 5; i < 250; i++){
  79. if(histogram[i]/width/(nrow/255.0) > 0.1){
  80. if(min == -1) min = i;
  81. max = i;
  82. }
  83. }
  84. // Stretch the brightness into the new range
  85. for(int y = 0; y < nrow; y++){
  86. for(int x = 0; x < width; x++){
  87. prow[y][x+offset] = (prow[y][x+offset]-min) / (max-min) * 255.0;
  88. prow[y][x+offset] = CLIP(prow[y][x+offset], 0.0, 255.0);
  89. }
  90. }
  91. }
  92. // Brightness calibrate, including telemetry
  93. void calibrateImage(float **prow, int nrow, int offset, int width, rgparam_t regr){
  94. offset -= SYNC_WIDTH+SPC_WIDTH;
  95. for (int y = 0; y < nrow; y++) {
  96. for (int x = 0; x < width+SYNC_WIDTH+SPC_WIDTH+TELE_WIDTH; x++) {
  97. float pv = rgcal(prow[y][x + offset], &regr);
  98. prow[y][x + offset] = CLIP(pv, 0, 255);
  99. }
  100. }
  101. }
  102. double teleNoise(double wedges[16]){
  103. double pattern[9] = { 31.07, 63.02, 94.96, 126.9, 158.86, 191.1, 228.62, 255.0, 0.0 };
  104. double noise = 0;
  105. for(int i = 0; i < 9; i++)
  106. noise += fabs(wedges[i] - pattern[i]);
  107. return noise;
  108. }
  109. // Get telemetry data for thermal calibration/equalization
  110. int calibrate(float **prow, int nrow, int offset, int width) {
  111. double teleline[MAX_HEIGHT] = { 0.0 };
  112. double wedge[16];
  113. rgparam_t regr[MAX_HEIGHT/FRAME_LEN + 1];
  114. int telestart, mtelestart = 0;
  115. int channel = -1;
  116. // The minimum rows required to decode a full frame
  117. if (nrow < 192) {
  118. fprintf(stderr, "Telemetry decoding error, not enough rows\n");
  119. return 0;
  120. }
  121. // Calculate average of a row of telemetry
  122. for (int y = 0; y < nrow; y++) {
  123. for (int x = 3; x < 43; x++)
  124. teleline[y] += prow[y][x + offset + width];
  125. teleline[y] /= 40.0;
  126. }
  127. /* Wedge 7 is white and 8 is black, this will have the largest
  128. * difference in brightness, this can be used to find the current
  129. * position within the telemetry.
  130. */
  131. float max = 0.0f;
  132. for (int 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. telestart = (mtelestart + 64) % FRAME_LEN;
  144. // Make sure that theres at least one full frame in the image
  145. if (nrow < telestart + FRAME_LEN) {
  146. fprintf(stderr, "Telemetry decoding error, not enough rows\n");
  147. return 0;
  148. }
  149. // Find the least noisy frame
  150. double minNoise = -1;
  151. int bestFrame = -1;
  152. for (int n = telestart, k = 0; n < nrow - FRAME_LEN; n += FRAME_LEN, k++) {
  153. int j;
  154. for (j = 0; j < 16; j++) {
  155. int i;
  156. wedge[j] = 0.0;
  157. for (i = 1; i < 7; i++)
  158. wedge[j] += teleline[n + j * 8 + i];
  159. wedge[j] /= 6;
  160. }
  161. double noise = teleNoise(wedge);
  162. if(noise < minNoise || minNoise == -1){
  163. minNoise = noise;
  164. bestFrame = k;
  165. // Compute & apply regression on the wedges
  166. rgcomp(wedge, &regr[k]);
  167. for (int j = 0; j < 16; j++)
  168. tele[j] = rgcal(wedge[j], &regr[k]);
  169. /* Compare the channel ID wedge to the reference
  170. * wedges, the wedge with the closest match will
  171. * be the channel ID
  172. */
  173. float min = -1;
  174. for (int j = 0; j < 6; j++) {
  175. float df = tele[15] - tele[j];
  176. df *= df;
  177. if (df < min || min == -1) {
  178. channel = j;
  179. min = df;
  180. }
  181. }
  182. // Find the brightness of the minute marker, I don't really know what for
  183. Cs = 0.0;
  184. int i, j = n;
  185. for (i = 0, j = n; j < n + FRAME_LEN; j++) {
  186. float csline = 0.0;
  187. for (int l = 3; l < 43; l++)
  188. csline += prow[n][l + offset - SPC_WIDTH];
  189. csline /= 40.0;
  190. if (csline > 50.0) {
  191. Cs += csline;
  192. i++;
  193. }
  194. }
  195. Cs = rgcal(Cs / i, &regr[k]);
  196. }
  197. }
  198. if(bestFrame == -1){
  199. fprintf(stderr, "Something has gone very wrong, please file a bug report.");
  200. return 0;
  201. }
  202. calibrateImage(prow, nrow, offset, width, regr[bestFrame]);
  203. return channel + 1;
  204. }
  205. void distrib(options_t *opts, image_t *img, char chid) {
  206. int max = 0;
  207. // Options
  208. options_t options;
  209. options.path = opts->path;
  210. options.effects = "";
  211. options.map = "";
  212. // Image options
  213. image_t distrib;
  214. strcpy(distrib.name, img->name);
  215. distrib.nrow = 256;
  216. // Assign memory
  217. for(int i = 0; i < 256; i++)
  218. distrib.prow[i] = (float *) malloc(sizeof(float) * 256);
  219. for(int n = 0; n < img->nrow; n++) {
  220. float *pixelv = img->prow[n];
  221. for(int i = 0; i < CH_WIDTH; i++) {
  222. int y = CLIP((int)pixelv[i + CHA_OFFSET], 0, 255);
  223. int x = CLIP((int)pixelv[i + CHB_OFFSET], 0, 255);
  224. distrib.prow[y][x]++;
  225. if(distrib.prow[y][x] > max)
  226. max = distrib.prow[y][x];
  227. }
  228. }
  229. // Scale to 0-255
  230. for(int x = 0; x < 256; x++)
  231. for(int y = 0; y < 256; y++)
  232. distrib.prow[y][x] = distrib.prow[y][x] / max * 255.0;
  233. extern int ImageOut(options_t *opts, image_t *img, int offset, int width, char *desc, char chid, char *palette);
  234. ImageOut(&options, &distrib, 0, 256, "Distribution", chid, NULL);
  235. }
  236. extern float quick_select(float arr[], int n);
  237. // Biased median denoise, pretyt ugly
  238. #define TRIG_LEVEL 40
  239. void denoise(float **prow, int nrow, int offset, int width){
  240. for(int y = 2; y < nrow-2; y++){
  241. for(int x = offset+1; x < offset+width-1; x++){
  242. if(prow[y][x+1] - prow[y][x] > TRIG_LEVEL ||
  243. prow[y][x-1] - prow[y][x] > TRIG_LEVEL ||
  244. prow[y+1][x] - prow[y][x] > TRIG_LEVEL ||
  245. prow[y-1][x] - prow[y][x] > TRIG_LEVEL){
  246. prow[y][x] = quick_select((float[]){
  247. prow[y+2][x-1], prow[y+2][x], prow[y+2][x+1],
  248. prow[y+1][x-1], prow[y+1][x], prow[y+1][x+1],
  249. prow[y-1][x-1], prow[y-1][x], prow[y-1][x+1],
  250. prow[y-2][x-1], prow[y-2][x], prow[y-2][x+1]
  251. }, 12);
  252. }
  253. }
  254. }
  255. }
  256. #undef TRIG_LEVEL
  257. // Flips a channe, for southbound passes
  258. void flipImage(image_t *img, int width, int offset){
  259. for(int y = 1; y < img->nrow; y++){
  260. for(int x = 1; x < ceil(width / 2.0); x++){
  261. // Flip top-left & bottom-right
  262. float buffer = img->prow[img->nrow - y][offset + x];
  263. img->prow[img->nrow - y][offset + x] = img->prow[y][offset + (width - x)];
  264. img->prow[y][offset + (width - x)] = buffer;
  265. }
  266. }
  267. }
  268. // --- Temperature Calibration --- //
  269. #include "satcal.h"
  270. typedef struct {
  271. double Nbb;
  272. double Cs;
  273. double Cb;
  274. int ch;
  275. } tempparam_t;
  276. // IR channel temperature compensation
  277. static void tempcomp(double t[16], int ch, int satnum, tempparam_t *tpr) {
  278. double Tbb, T[4];
  279. double C;
  280. tpr->ch = ch - 4;
  281. // Compute equivalent T blackbody temperature
  282. for (int n = 0; n < 4; n++) {
  283. float d0, d1, d2;
  284. C = t[9 + n] * 4.0;
  285. d0 = satcal[satnum].d[n][0];
  286. d1 = satcal[satnum].d[n][1];
  287. d2 = satcal[satnum].d[n][2];
  288. T[n] = d0;
  289. T[n] += d1 * C;
  290. C *= C;
  291. T[n] += d2 * C;
  292. }
  293. Tbb = (T[0] + T[1] + T[2] + T[3]) / 4.0;
  294. Tbb = satcal[satnum].rad[tpr->ch].A + satcal[satnum].rad[tpr->ch].B * Tbb;
  295. // Compute blackbody radiance temperature
  296. C = satcal[satnum].rad[tpr->ch].vc;
  297. tpr->Nbb = c1 * C * C * C / (expm1(c2 * C / Tbb));
  298. // Store blackbody count and space
  299. tpr->Cs = Cs * 4.0;
  300. tpr->Cb = t[14] * 4.0;
  301. }
  302. // IR channel temperature calibration
  303. static double tempcal(float Ce, int satnum, tempparam_t * rgpr) {
  304. double Nl, Nc, Ns, Ne;
  305. double T, vc;
  306. Ns = satcal[satnum].cor[rgpr->ch].Ns;
  307. Nl = Ns + (rgpr->Nbb - Ns) * (rgpr->Cs - Ce * 4.0) / (rgpr->Cs - rgpr->Cb);
  308. Nc = satcal[satnum].cor[rgpr->ch].b[0] +
  309. satcal[satnum].cor[rgpr->ch].b[1] * Nl +
  310. satcal[satnum].cor[rgpr->ch].b[2] * Nl * Nl;
  311. Ne = Nl + Nc;
  312. vc = satcal[satnum].rad[rgpr->ch].vc;
  313. T = c2 * vc / log(c1 * vc * vc * vc / Ne + 1.0);
  314. T = (T - satcal[satnum].rad[rgpr->ch].A) / satcal[satnum].rad[rgpr->ch].B;
  315. // Convert to celsius
  316. T -= 273.15;
  317. // Rescale to 0-255 for -100°C to +60°C
  318. T = (T + 100.0) / 160.0 * 255.0;
  319. return T;
  320. }
  321. // Temperature calibration wrapper
  322. void temperature(options_t *opts, image_t *img, int offset, int width){
  323. tempparam_t temp;
  324. printf("Temperature... ");
  325. fflush(stdout);
  326. tempcomp(tele, img->chB, opts->satnum - 15, &temp);
  327. for (int y = 0; y < img->nrow; y++) {
  328. for (int x = 0; x < width; x++) {
  329. img->prow[y][x + offset] = tempcal(img->prow[y][x + offset], opts->satnum - 15, &temp);
  330. }
  331. }
  332. printf("Done\n");
  333. }