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.
 
 
 
 
 

413 lines
11 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 <math.h>
  22. #include <stdlib.h>
  23. #include "apt.h"
  24. #include "libs/reg.h"
  25. #include "image.h"
  26. #define REGORDER 3
  27. typedef struct {
  28. double cf[REGORDER + 1];
  29. } rgparam_t;
  30. // Compute regression
  31. static void rgcomp(double x[16], rgparam_t * 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_t *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. void apt_histogramEqualise(float **prow, int nrow, int offset, int width){
  49. // Plot histogram
  50. int histogram[256] = { 0 };
  51. for(int y = 0; y < nrow; y++)
  52. for(int x = 0; x < width; x++)
  53. histogram[(int)CLIP(prow[y][x+offset], 0, 255)]++;
  54. // Calculate cumulative frequency
  55. long sum = 0, cf[256] = { 0 };
  56. for(int i = 0; i < 255; i++){
  57. sum += histogram[i];
  58. cf[i] = sum;
  59. }
  60. // Apply histogram
  61. int area = nrow * width;
  62. for(int y = 0; y < nrow; y++){
  63. for(int x = 0; x < width; x++){
  64. int k = (int)prow[y][x+offset];
  65. prow[y][x+offset] = (256.0f/area) * cf[k];
  66. }
  67. }
  68. }
  69. void apt_linearEnhance(float **prow, int nrow, int offset, int width){
  70. // Plot histogram
  71. int histogram[256] = { 0 };
  72. for(int y = 0; y < nrow; y++)
  73. for(int x = 0; x < width; x++)
  74. histogram[(int)CLIP(prow[y][x+offset], 0, 255)]++;
  75. // Find min/max points
  76. int min = -1, max = -1;
  77. for(int i = 5; i < 250; i++){
  78. if(histogram[i]/width/(nrow/255.0) > 0.1){
  79. if(min == -1) min = i;
  80. max = i;
  81. }
  82. }
  83. // Stretch the brightness into the new range
  84. for(int y = 0; y < nrow; y++){
  85. for(int x = 0; x < width; x++){
  86. prow[y][x+offset] = (prow[y][x+offset]-min) / (max-min) * 255.0f;
  87. prow[y][x+offset] = CLIP(prow[y][x+offset], 0.0f, 255.0f);
  88. }
  89. }
  90. }
  91. // Brightness calibrate, including telemetry
  92. void calibrateImage(float **prow, int nrow, int offset, int width, rgparam_t regr){
  93. offset -= APT_SYNC_WIDTH+APT_SPC_WIDTH;
  94. for (int y = 0; y < nrow; y++) {
  95. for (int x = 0; x < width+APT_SYNC_WIDTH+APT_SPC_WIDTH+APT_TELE_WIDTH; x++) {
  96. float pv = (float)rgcal(prow[y][x + offset], &regr);
  97. prow[y][x + offset] = CLIP(pv, 0, 255);
  98. }
  99. }
  100. }
  101. double teleNoise(double wedges[16]){
  102. double pattern[9] = { 31.07, 63.02, 94.96, 126.9, 158.86, 191.1, 228.62, 255.0, 0.0 };
  103. double noise = 0;
  104. for(int i = 0; i < 9; i++)
  105. noise += fabs(wedges[i] - pattern[i]);
  106. return noise;
  107. }
  108. // Get telemetry data for thermal calibration
  109. int apt_calibrate(float **prow, int nrow, int offset, int width) {
  110. double teleline[APT_MAX_HEIGHT] = { 0.0 };
  111. double wedge[16];
  112. rgparam_t regr[APT_MAX_HEIGHT/APT_FRAME_LEN + 1];
  113. int telestart, mtelestart = 0;
  114. int channel = -1;
  115. // The minimum rows required to decode a full frame
  116. if (nrow < 192) {
  117. fprintf(stderr, "Telemetry decoding error, not enough rows\n");
  118. return 0;
  119. }
  120. // Calculate average of a row of telemetry
  121. for (int y = 0; y < nrow; y++) {
  122. for (int x = 3; x < 43; x++)
  123. teleline[y] += prow[y][x + offset + width];
  124. teleline[y] /= 40.0;
  125. }
  126. /* Wedge 7 is white and 8 is black, this will have the largest
  127. * difference in brightness, this can be used to find the current
  128. * position within the telemetry.
  129. */
  130. float max = 0.0f;
  131. for (int n = nrow / 3 - 64; n < 2 * nrow / 3 - 64; n++) {
  132. float df;
  133. // (sum 4px below) - (sum 4px above)
  134. df = (float)((teleline[n - 4] + teleline[n - 3] + teleline[n - 2] + teleline[n - 1]) -
  135. (teleline[n + 0] + teleline[n + 1] + teleline[n + 2] + teleline[n + 3]));
  136. // Find the maximum difference
  137. if (df > max) {
  138. mtelestart = n;
  139. max = df;
  140. }
  141. }
  142. telestart = (mtelestart + 64) % APT_FRAME_LEN;
  143. // Make sure that theres at least one full frame in the image
  144. if (nrow < telestart + APT_FRAME_LEN) {
  145. fprintf(stderr, "Telemetry decoding error, not enough rows\n");
  146. return 0;
  147. }
  148. // Find the least noisy frame
  149. double minNoise = -1;
  150. int bestFrame = -1;
  151. for (int n = telestart, k = 0; n < nrow - APT_FRAME_LEN; n += APT_FRAME_LEN, k++) {
  152. int j;
  153. for (j = 0; j < 16; j++) {
  154. int i;
  155. wedge[j] = 0.0;
  156. for (i = 1; i < 7; i++)
  157. wedge[j] += teleline[n + j * 8 + i];
  158. wedge[j] /= 6;
  159. }
  160. double noise = teleNoise(wedge);
  161. if(noise < minNoise || minNoise == -1){
  162. minNoise = noise;
  163. bestFrame = k;
  164. // Compute & apply regression on the wedges
  165. rgcomp(wedge, &regr[k]);
  166. for (int j = 0; j < 16; j++)
  167. tele[j] = (float)rgcal((float)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 = -1;
  173. for (int j = 0; j < 6; j++) {
  174. float df = (float)(tele[15] - tele[j]);
  175. df *= df;
  176. if (df < min || min == -1) {
  177. channel = j;
  178. min = df;
  179. }
  180. }
  181. // Find the brightness of the minute marker, I don't really know what for
  182. Cs = 0.0;
  183. int i, j = n;
  184. for (i = 0, j = n; j < n + APT_FRAME_LEN; j++) {
  185. float csline = 0.0;
  186. for (int l = 3; l < 43; l++)
  187. csline += prow[n][l + offset - APT_SPC_WIDTH];
  188. csline /= 40.0;
  189. if (csline > 50.0) {
  190. Cs += csline;
  191. i++;
  192. }
  193. }
  194. Cs = rgcal((float)(Cs / i), &regr[k]);
  195. }
  196. }
  197. if(bestFrame == -1){
  198. fprintf(stderr, "Something has gone very wrong, please file a bug report.\n");
  199. return 0;
  200. }
  201. calibrateImage(prow, nrow, offset, width, regr[bestFrame]);
  202. return channel + 1;
  203. }
  204. extern float quick_select(float arr[], int n);
  205. // Biased median denoise, pretyt ugly
  206. #define TRIG_LEVEL 40
  207. void apt_denoise(float **prow, int nrow, int offset, int width){
  208. for(int y = 2; y < nrow-2; y++){
  209. for(int x = offset+1; x < offset+width-1; x++){
  210. if(prow[y][x+1] - prow[y][x] > TRIG_LEVEL ||
  211. prow[y][x-1] - prow[y][x] > TRIG_LEVEL ||
  212. prow[y+1][x] - prow[y][x] > TRIG_LEVEL ||
  213. prow[y-1][x] - prow[y][x] > TRIG_LEVEL){
  214. prow[y][x] = quick_select((float[]){
  215. prow[y+2][x-1], prow[y+2][x], prow[y+2][x+1],
  216. prow[y+1][x-1], prow[y+1][x], prow[y+1][x+1],
  217. prow[y-1][x-1], prow[y-1][x], prow[y-1][x+1],
  218. prow[y-2][x-1], prow[y-2][x], prow[y-2][x+1]
  219. }, 12);
  220. }
  221. }
  222. }
  223. }
  224. #undef TRIG_LEVEL
  225. // Flips a channel, for northbound passes
  226. void apt_flipImage(apt_image_t *img, int width, int offset){
  227. for(int y = 1; y < img->nrow; y++){
  228. for(int x = 1; x < ceil(width / 2.0); x++){
  229. // Flip top-left & bottom-right
  230. float buffer = img->prow[img->nrow - y][offset + x];
  231. img->prow[img->nrow - y][offset + x] = img->prow[y][offset + (width - x)];
  232. img->prow[y][offset + (width - x)] = buffer;
  233. }
  234. }
  235. }
  236. // Calculate crop to reomve noise from the start and end of an image
  237. int apt_cropNoise(apt_image_t *img){
  238. #define NOISE_THRESH 180.0
  239. // Average value of minute marker
  240. float spc_rows[APT_MAX_HEIGHT] = { 0.0 };
  241. int startCrop = 0; int endCrop = img->nrow;
  242. for(int y = 0; y < img->nrow; y++) {
  243. for(int x = 0; x < APT_SPC_WIDTH; x++) {
  244. spc_rows[y] += img->prow[y][x + (APT_CHB_OFFSET - APT_SPC_WIDTH)];
  245. }
  246. spc_rows[y] /= APT_SPC_WIDTH;
  247. // Skip minute markings
  248. if(spc_rows[y] < 10) {
  249. spc_rows[y] = spc_rows[y-1];
  250. }
  251. }
  252. // 3 row average
  253. for(int y = 0; y < img->nrow; y++){
  254. spc_rows[y] = (spc_rows[y+1] + spc_rows[y+2] + spc_rows[y+3])/3;
  255. //img.prow[y][0] = spc_rows[y];
  256. }
  257. // Find ends
  258. for(int y = 0; y < img->nrow-1; y++) {
  259. if(spc_rows[y] > NOISE_THRESH){
  260. endCrop = y;
  261. }
  262. }
  263. for(int y = img->nrow; y > 0; y--) {
  264. if(spc_rows[y] > NOISE_THRESH) {
  265. startCrop = y;
  266. }
  267. }
  268. //printf("Crop rows: %i -> %i\n", startCrop, endCrop);
  269. // Remove the noisy rows at start
  270. for(int y = 0; y < img->nrow-startCrop; y++) {
  271. memmove(img->prow[y], img->prow[y+startCrop], sizeof(float)*APT_PROW_WIDTH);
  272. }
  273. // Ignore the noisy rows at the end
  274. img->nrow = (endCrop - startCrop);
  275. return startCrop;
  276. }
  277. // --- Temperature Calibration --- //
  278. #include "satcal.h"
  279. typedef struct {
  280. double Nbb;
  281. double Cs;
  282. double Cb;
  283. int ch;
  284. } tempparam_t;
  285. // IR channel temperature compensation
  286. static void tempcomp(double t[16], int ch, int satnum, tempparam_t *tpr) {
  287. double Tbb, T[4];
  288. double C;
  289. tpr->ch = ch - 4;
  290. // Compute equivalent T blackbody temperature
  291. for (int n = 0; n < 4; n++) {
  292. float d0, d1, d2;
  293. C = t[9 + n] * 4.0;
  294. d0 = satcal[satnum].d[n][0];
  295. d1 = satcal[satnum].d[n][1];
  296. d2 = satcal[satnum].d[n][2];
  297. T[n] = d0;
  298. T[n] += d1 * C;
  299. C *= C;
  300. T[n] += d2 * C;
  301. }
  302. Tbb = (T[0] + T[1] + T[2] + T[3]) / 4.0;
  303. Tbb = satcal[satnum].rad[tpr->ch].A + satcal[satnum].rad[tpr->ch].B * Tbb;
  304. // Compute blackbody radiance temperature
  305. C = satcal[satnum].rad[tpr->ch].vc;
  306. tpr->Nbb = c1 * C * C * C / (expm1(c2 * C / Tbb));
  307. // Store blackbody count and space
  308. tpr->Cs = Cs * 4.0;
  309. tpr->Cb = t[14] * 4.0;
  310. }
  311. // IR channel temperature calibration
  312. static double tempcal(float Ce, int satnum, tempparam_t * rgpr) {
  313. double Nl, Nc, Ns, Ne;
  314. double T, vc;
  315. Ns = satcal[satnum].cor[rgpr->ch].Ns;
  316. Nl = Ns + (rgpr->Nbb - Ns) * (rgpr->Cs - Ce * 4.0) / (rgpr->Cs - rgpr->Cb);
  317. Nc = satcal[satnum].cor[rgpr->ch].b[0] +
  318. satcal[satnum].cor[rgpr->ch].b[1] * Nl +
  319. satcal[satnum].cor[rgpr->ch].b[2] * Nl * Nl;
  320. Ne = Nl + Nc;
  321. vc = satcal[satnum].rad[rgpr->ch].vc;
  322. T = c2 * vc / log(c1 * vc * vc * vc / Ne + 1.0);
  323. T = (T - satcal[satnum].rad[rgpr->ch].A) / satcal[satnum].rad[rgpr->ch].B;
  324. // Convert to celsius
  325. T -= 273.15;
  326. // Rescale to 0-255 for -120°C to +40°C
  327. T = (T + 100.0) / 160.0 * 255.0;
  328. return T;
  329. }
  330. // Temperature calibration wrapper
  331. void temperature(options_t *opts, apt_image_t *img, int offset, int width){
  332. tempparam_t temp;
  333. printf("Temperature... ");
  334. fflush(stdout);
  335. tempcomp(tele, img->chB, opts->satnum - 15, &temp);
  336. for (int y = 0; y < img->nrow; y++) {
  337. for (int x = 0; x < width; x++) {
  338. img->prow[y][x + offset] = (float)tempcal(img->prow[y][x + offset], opts->satnum - 15, &temp);
  339. }
  340. }
  341. printf("Done\n");
  342. }