1 /*
2 * This file is part of Siril, an astronomy image processor.
3 * Copyright (C) 2005-2011 Francois Meyer (dulle at free.fr)
4 * Copyright (C) 2012-2021 team free-astro (see more in AUTHORS file)
5 * Reference site is https://free-astro.org/index.php/Siril
6 *
7 * Siril is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * Siril 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 General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with Siril. If not, see <http://www.gnu.org/licenses/>.
19 */
20
21 #include <string.h>
22 #include <math.h>
23 #include <stdint.h>
24 #include <assert.h>
25
26 #include "core/siril.h"
27 #include "core/proto.h"
28 #include "core/processing.h"
29 #include "core/command.h"
30 #include "gui/utils.h"
31 #include "gui/progress_and_log.h"
32 #include "gui/message_dialog.h"
33 #include "gui/dialogs.h"
34 #include "io/sequence.h"
35 #include "io/fits_sequence.h"
36 #include "io/seqwriter.h"
37 #include "io/image_format_fits.h"
38 #include "io/conversion.h"
39 #include "algos/demosaicing.h"
40 #include "algos/statistics.h"
41
42 #define USE_SIRIL_DEBAYER FALSE
43
44 const char *filter_pattern[] = {
45 "RGGB",
46 "BGGR",
47 "GBRG",
48 "GRBG",
49
50 /* XTRANS */
51 "GGRGGB" // ----> XTRANS_1
52 "GGBGGR"
53 "BRGRBG"
54 "GGBGGR"
55 "GGRGGB"
56 "RBGBRG",
57
58 "RBGBRG" // ----> XTRANS_2
59 "GGRGGB"
60 "GGBGGR"
61 "BRGRBG"
62 "GGBGGR"
63 "GGRGGB",
64
65 "GRGGBG"
66 "BGBRGR"
67 "GRGGBG"
68 "GBGGRG"
69 "RGTBGB"
70 "GBGGRG",
71
72 "GBGGRG"
73 "RGRBGB"
74 "GBGGRG"
75 "GRGGBG"
76 "BGBRGR"
77 "GRGGBG"
78 };
79
80 /** Calculate the bayer pattern color from the row and column **/
FC(const size_t row,const size_t col,const uint32_t filters)81 static inline int FC(const size_t row, const size_t col, const uint32_t filters) {
82 return filters >> (((row << 1 & 14) + (col & 1)) << 1) & 3;
83 }
84
85 /* width and height are sizes of the original image */
super_pixel_ushort(const WORD * buf,WORD * newbuf,int width,int height,sensor_pattern pattern)86 static void super_pixel_ushort(const WORD *buf, WORD *newbuf, int width, int height,
87 sensor_pattern pattern) {
88 long i = 0;
89 for (int row = 0; row < height - 1; row += 2) {
90 for (int col = 0; col < width - 1; col += 2) {
91 float tmp;
92 switch (pattern) {
93 default:
94 case BAYER_FILTER_RGGB:
95 newbuf[i + 0] = buf[col + row * width];
96 tmp = buf[1 + col + row * width];
97 tmp += buf[col + (1 + row) * width];
98 newbuf[i + 1] = round_to_WORD(tmp * 0.5f);
99 newbuf[i + 2] = buf[1 + col + (1 + row) * width];
100 break;
101 case BAYER_FILTER_BGGR:
102 newbuf[i + 2] = buf[col + row * width];
103 tmp = buf[1 + col + row * width];
104 tmp += buf[(col + row * width) + width];
105 newbuf[i + 1] = round_to_WORD(tmp * 0.5f);
106 newbuf[i + 0] = buf[(1 + col + row * width) + width];
107 break;
108 case BAYER_FILTER_GBRG:
109 newbuf[i + 2] = buf[1 + col + row * width];
110 newbuf[i + 0] = buf[(col + row * width) + width];
111 tmp = buf[col + row * width];
112 tmp += buf[(1 + col + row * width) + width];
113 newbuf[i + 1] = round_to_WORD(tmp * 0.5f);
114 break;
115 case BAYER_FILTER_GRBG:
116 newbuf[i + 0] = buf[1 + col + row * width];
117 newbuf[i + 2] = buf[(col + row * width) + width];
118 tmp = buf[col + row * width];
119 tmp += buf[(1 + col + row * width) + width];
120 newbuf[i + 1] = round_to_WORD(tmp * 0.5f);
121 break;
122 }
123 i += 3;
124 }
125 }
126 }
127
128 /* width and height are sizes of the original image */
super_pixel_float(const float * buf,float * newbuf,int width,int height,sensor_pattern pattern)129 static void super_pixel_float(const float *buf, float *newbuf, int width, int height,
130 sensor_pattern pattern) {
131 long i = 0;
132 for (int row = 0; row < height - 1; row += 2) {
133 for (int col = 0; col < width - 1; col += 2) {
134 float tmp;
135 switch (pattern) {
136 default:
137 case BAYER_FILTER_RGGB:
138 newbuf[i + 0] = buf[col + row * width];
139 tmp = buf[1 + col + row * width];
140 tmp += buf[col + (1 + row) * width];
141 newbuf[i + 1] = tmp * 0.5f;
142 newbuf[i + 2] = buf[1 + col + (1 + row) * width];
143 break;
144 case BAYER_FILTER_BGGR:
145 newbuf[i + 2] = buf[col + row * width];
146 tmp = buf[1 + col + row * width];
147 tmp += buf[(col + row * width) + width];
148 newbuf[i + 1] = tmp * 0.5f;
149 newbuf[i + 0] = buf[(1 + col + row * width) + width];
150 break;
151 case BAYER_FILTER_GBRG:
152 newbuf[i + 2] = buf[1 + col + row * width];
153 newbuf[i + 0] = buf[(col + row * width) + width];
154 tmp = buf[col + row * width];
155 tmp += buf[(1 + col + row * width) + width];
156 newbuf[i + 1] = tmp * 0.5f;
157 break;
158 case BAYER_FILTER_GRBG:
159 newbuf[i + 0] = buf[1 + col + row * width];
160 newbuf[i + 2] = buf[(col + row * width) + width];
161 tmp = buf[col + row * width];
162 tmp += buf[(1 + col + row * width) + width];
163 newbuf[i + 1] = tmp * 0.5f;
164 break;
165 }
166 i += 3;
167 }
168 }
169 }
170
171 /***************************************************
172 *
173 * Written by Damien Douxchamps and Frederic Devernay
174 * The original VNG and AHD Bayer decoding are from Dave Coffin's DCRAW.
175 * https://code.google.com/p/gst-plugins-elphel/
176 *
177 * *************************************************/
178
ClearBorders(WORD * rgb,int sx,int sy,int w)179 static void ClearBorders(WORD *rgb, int sx, int sy, int w) {
180 int i, j;
181
182 /* black edges: */
183 i = 3 * sx * w - 1;
184 j = 3 * sx * sy - 1;
185 while (i >= 0) {
186 rgb[i--] = 0;
187 rgb[j--] = 0;
188 }
189
190 int low = sx * (w - 1) * 3 - 1 + w * 3;
191 i = low + sx * (sy - w * 2 + 1) * 3;
192 while (i > low) {
193 j = 6 * w;
194 while (j > 0) {
195 rgb[i--] = 0;
196 j--;
197 }
198 i -= (sx - 2 * w) * 3;
199 }
200 }
201
202 /* OpenCV's Bayer decoding */
bayer_Bilinear(const WORD * bayer,WORD * rgb,int sx,int sy,sensor_pattern tile)203 static int bayer_Bilinear(const WORD *bayer, WORD *rgb, int sx, int sy,
204 sensor_pattern tile) {
205 const int bayerStep = sx;
206 const int rgbStep = 3 * sx;
207 int width = sx;
208 int height = sy;
209 int blue = tile == BAYER_FILTER_BGGR || tile == BAYER_FILTER_GBRG ? -1 : 1;
210 int start_with_green = tile == BAYER_FILTER_GBRG
211 || tile == BAYER_FILTER_GRBG;
212
213 if (tile > BAYER_FILTER_MAX || tile < BAYER_FILTER_MIN)
214 return -1;
215
216 ClearBorders(rgb, sx, sy, 1);
217 rgb += rgbStep + 3 + 1;
218 height -= 2;
219 width -= 2;
220
221 for (; height--; bayer += bayerStep, rgb += rgbStep) {
222 int t0, t1;
223 const WORD *bayerEnd = bayer + width;
224
225 if (start_with_green) {
226 t0 = (bayer[1] + bayer[bayerStep * 2 + 1] + 1) >> 1;
227 t1 = (bayer[bayerStep] + bayer[bayerStep + 2] + 1) >> 1;
228 rgb[-blue] = round_to_WORD(t0);
229 rgb[0] = bayer[bayerStep + 1];
230 rgb[blue] = round_to_WORD(t1);
231 bayer++;
232 rgb += 3;
233 }
234
235 if (blue > 0) {
236 for (; bayer <= bayerEnd - 2; bayer += 2, rgb += 6) {
237 t0 = (bayer[0] + bayer[2] + bayer[bayerStep * 2]
238 + bayer[bayerStep * 2 + 2] + 2) >> 2;
239 t1 = (bayer[1] + bayer[bayerStep] + bayer[bayerStep + 2]
240 + bayer[bayerStep * 2 + 1] + 2) >> 2;
241 rgb[-1] = round_to_WORD(t0);
242 rgb[0] = round_to_WORD(t1);
243 rgb[1] = bayer[bayerStep + 1];
244
245 t0 = (bayer[2] + bayer[bayerStep * 2 + 2] + 1) >> 1;
246 t1 = (bayer[bayerStep + 1] + bayer[bayerStep + 3] + 1) >> 1;
247 rgb[2] = round_to_WORD(t0);
248 rgb[3] = bayer[bayerStep + 2];
249 rgb[4] = round_to_WORD(t1);
250 }
251 } else {
252 for (; bayer <= bayerEnd - 2; bayer += 2, rgb += 6) {
253 t0 = (bayer[0] + bayer[2] + bayer[bayerStep * 2]
254 + bayer[bayerStep * 2 + 2] + 2) >> 2;
255 t1 = (bayer[1] + bayer[bayerStep] + bayer[bayerStep + 2]
256 + bayer[bayerStep * 2 + 1] + 2) >> 2;
257 rgb[1] = round_to_WORD(t0);
258 rgb[0] = round_to_WORD(t1);
259 rgb[-1] = bayer[bayerStep + 1];
260
261 t0 = (bayer[2] + bayer[bayerStep * 2 + 2] + 1) >> 1;
262 t1 = (bayer[bayerStep + 1] + bayer[bayerStep + 3] + 1) >> 1;
263 rgb[4] = round_to_WORD(t0);
264 rgb[3] = bayer[bayerStep + 2];
265 rgb[2] = round_to_WORD(t1);
266 }
267 }
268
269 if (bayer < bayerEnd) {
270 t0 = (bayer[0] + bayer[2] + bayer[bayerStep * 2]
271 + bayer[bayerStep * 2 + 2] + 2) >> 2;
272 t1 = (bayer[1] + bayer[bayerStep] + bayer[bayerStep + 2]
273 + bayer[bayerStep * 2 + 1] + 2) >> 2;
274 rgb[-blue] = round_to_WORD(t0);
275 rgb[0] = round_to_WORD(t1);
276 rgb[blue] = bayer[bayerStep + 1];
277 bayer++;
278 rgb += 3;
279 }
280
281 bayer -= width;
282 rgb -= width * 3;
283
284 blue = -blue;
285 start_with_green = !start_with_green;
286 }
287
288 return 0;
289 }
290
291 #define ABSOLU(x) (((int)(x) ^ ((int)(x) >> 31)) - ((int)(x) >> 31))
292
bayer_VNG(const WORD * bayer,WORD * dst,int sx,int sy,sensor_pattern pattern)293 static int bayer_VNG(const WORD *bayer, WORD *dst, int sx, int sy,
294 sensor_pattern pattern) {
295 const signed char bayervng_terms[] = { -2, -2, +0, -1, 0, 0x01, -2, -2,
296 +0, +0, 1, 0x01, -2, -1, -1, +0, 0, 0x01, -2, -1, +0, -1, 0, 0x02, -2,
297 -1, +0, +0, 0, 0x03, -2, -1, +0, +1, 1, 0x01, -2, +0, +0, -1, 0, 0x06,
298 -2, +0, +0, +0, 1, 0x02, -2, +0, +0, +1, 0, 0x03, -2, +1, -1, +0, 0,
299 0x04, -2, +1, +0, -1, 1, 0x04, -2, +1, +0, +0, 0, 0x06, -2, +1, +0, +1,
300 0, 0x02, -2, +2, +0, +0, 1, 0x04, -2, +2, +0, +1, 0, 0x04, -1, -2, -1,
301 +0, 0, 0x80, -1, -2, +0, -1, 0, 0x01, -1, -2, +1, -1, 0, 0x01, -1, -2,
302 +1, +0, 1, 0x01, -1, -1, -1, +1, 0, 0x88, -1, -1, +1, -2, 0, 0x40, -1,
303 -1, +1, -1, 0, 0x22, -1, -1, +1, +0, 0, 0x33, -1, -1, +1, +1, 1, 0x11,
304 -1, +0, -1, +2, 0, 0x08, -1, +0, +0, -1, 0, 0x44, -1, +0, +0, +1, 0,
305 0x11, -1, +0, +1, -2, 1, 0x40, -1, +0, +1, -1, 0, 0x66, -1, +0, +1, +0,
306 1, 0x22, -1, +0, +1, +1, 0, 0x33, -1, +0, +1, +2, 1, 0x10, -1, +1, +1,
307 -1, 1, 0x44, -1, +1, +1, +0, 0, 0x66, -1, +1, +1, +1, 0, 0x22, -1, +1,
308 +1, +2, 0, 0x10, -1, +2, +0, +1, 0, 0x04, -1, +2, +1, +0, 1, 0x04, -1,
309 +2, +1, +1, 0, 0x04, +0, -2, +0, +0, 1, 0x80, +0, -1, +0, +1, 1, 0x88,
310 +0, -1, +1, -2, 0, 0x40, +0, -1, +1, +0, 0, 0x11, +0, -1, +2, -2, 0,
311 0x40, +0, -1, +2, -1, 0, 0x20, +0, -1, +2, +0, 0, 0x30, +0, -1, +2, +1,
312 1, 0x10, +0, +0, +0, +2, 1, 0x08, +0, +0, +2, -2, 1, 0x40, +0, +0, +2,
313 -1, 0, 0x60, +0, +0, +2, +0, 1, 0x20, +0, +0, +2, +1, 0, 0x30, +0, +0,
314 +2, +2, 1, 0x10, +0, +1, +1, +0, 0, 0x44, +0, +1, +1, +2, 0, 0x10, +0,
315 +1, +2, -1, 1, 0x40, +0, +1, +2, +0, 0, 0x60, +0, +1, +2, +1, 0, 0x20,
316 +0, +1, +2, +2, 0, 0x10, +1, -2, +1, +0, 0, 0x80, +1, -1, +1, +1, 0,
317 0x88, +1, +0, +1, +2, 0, 0x08, +1, +0, +2, -1, 0, 0x40, +1, +0, +2, +1,
318 0, 0x10 }, bayervng_chood[] = { -1, -1, -1, 0, -1, +1, 0, +1, +1, +1,
319 +1, 0, +1, -1, 0, -1 };
320 const int height = sy, width = sx;
321 const signed char *cp;
322 WORD (*brow[5])[3], *pix; /* [FD] */
323 int code[8][2][320], *ip, gval[8], gmin, gmax, sum[4];
324 int row, col, x, y, x1, x2, y1, y2, t, weight, grads, color, diag;
325 int g, diff, thold, num, c;
326 unsigned long int filters; /* [FD] */
327
328 /* first, use bilinear bayer decoding */
329
330 bayer_Bilinear(bayer, dst, sx, sy, pattern);
331
332 switch (pattern) {
333 case BAYER_FILTER_BGGR:
334 filters = 0x16161616;
335 break;
336 case BAYER_FILTER_GRBG:
337 filters = 0x61616161;
338 break;
339 case BAYER_FILTER_RGGB:
340 filters = 0x94949494;
341 break;
342 case BAYER_FILTER_GBRG:
343 filters = 0x49494949;
344 break;
345 default:
346 return -1;
347 }
348
349 for (row = 0; row < 8; row++) { /* Precalculate for VNG */
350 for (col = 0; col < 2; col++) {
351 ip = code[row][col];
352 for (cp = bayervng_terms, t = 0; t < 64; t++) {
353 y1 = *cp++;
354 x1 = *cp++;
355 y2 = *cp++;
356 x2 = *cp++;
357 weight = *cp++;
358 grads = *cp++;
359 color = FC(row + y1, col + x1, filters);
360 if (FC(row+y2,col+x2, filters) != (unsigned int) color)
361 continue;
362 diag = (FC(row,col+1, filters) == (unsigned int) color
363 && FC(row+1,col, filters) == (unsigned int) color) ? 2 : 1;
364 if (abs(y1 - y2) == diag && abs(x1 - x2) == diag)
365 continue;
366 *ip++ = (y1 * width + x1) * 3 + color; /* [FD] */
367 *ip++ = (y2 * width + x2) * 3 + color; /* [FD] */
368 *ip++ = weight;
369 for (g = 0; g < 8; g++)
370 if (grads & 1 << g)
371 *ip++ = g;
372 *ip++ = -1;
373 }
374 *ip++ = INT_MAX;
375 for (cp = bayervng_chood, g = 0; g < 8; g++) {
376 y = *cp++;
377 x = *cp++;
378 *ip++ = (y * width + x) * 3; /* [FD] */
379 color = FC(row, col, filters);
380 if (FC(row+y,col+x, filters) != (unsigned int) color
381 && FC(row+y*2,col+x*2, filters) == (unsigned int) color)
382 *ip++ = (y * width + x) * 6 + color; /* [FD] */
383 else
384 *ip++ = 0;
385 }
386 }
387 }
388 brow[4] = calloc(width * 3, sizeof **brow);
389 //merror (brow[4], "vng_interpolate()");
390 for (row = 0; row < 3; row++)
391 brow[row] = brow[4] + row * width;
392 for (row = 2; row < height - 2; row++) { /* Do VNG interpolation */
393 for (col = 2; col < width - 2; col++) {
394 pix = dst + (row * width + col) * 3; /* [FD] */
395 ip = code[row & 7][col & 1];
396 memset(gval, 0, sizeof gval);
397 while ((g = ip[0]) != INT_MAX) { /* Calculate gradients */
398 diff = ABSOLU(pix[g] - pix[ip[1]]) << ip[2];
399 gval[ip[3]] += diff;
400 ip += 5;
401 if ((g = ip[-1]) == -1)
402 continue;
403 gval[g] += diff;
404 while ((g = *ip++) != -1)
405 gval[g] += diff;
406 }
407 ip++;
408 gmin = gmax = gval[0]; /* Choose a threshold */
409 for (g = 1; g < 8; g++) {
410 if (gmin > gval[g])
411 gmin = gval[g];
412 if (gmax < gval[g])
413 gmax = gval[g];
414 }
415 if (gmax == 0) {
416 memcpy(brow[2][col], pix, 3 * sizeof *dst); /* [FD] */
417 continue;
418 }
419 thold = gmin + (gmax >> 1);
420 memset(sum, 0, sizeof sum);
421 color = FC(row, col, filters);
422 for (num = g = 0; g < 8; g++, ip += 2) { /* Average the neighbors */
423 if (gval[g] <= thold) {
424 for (c = 0; c < 3; c++) /* [FD] */
425 if (c == color && ip[1])
426 sum[c] += (pix[c] + pix[ip[1]]) >> 1;
427 else
428 sum[c] += pix[ip[0] + c];
429 num++;
430 }
431 }
432 for (c = 0; c < 3; c++) { /* [FD] Save to buffer */
433 t = pix[color];
434 if (c != color)
435 t += (sum[c] - sum[color]) / num;
436 //~ CLIP16(t,brow[2][col][c], 16); /* [FD] */
437 brow[2][col][c] = round_to_WORD(t); /* [FD] */
438 }
439 }
440 if (row > 3) /* Write buffer to image */
441 memcpy(dst + 3 * ((row - 2) * width + 2), brow[0] + 2,
442 (width - 4) * 3 * sizeof *dst); /* [FD] */
443 for (g = 0; g < 4; g++)
444 brow[(g - 1) & 3] = brow[g];
445 }
446 memcpy(dst + 3 * ((row - 2) * width + 2), brow[0] + 2,
447 (width - 4) * 3 * sizeof *dst);
448 memcpy(dst + 3 * ((row - 1) * width + 2), brow[1] + 2,
449 (width - 4) * 3 * sizeof *dst);
450 free(brow[4]);
451
452 return 0;
453 }
454
455 /* AHD interpolation ported from dcraw to libdc1394 by Samuel Audet */
456 static gboolean ahd_inited = FALSE; /* WARNING: not multi-processor safe */
457
458 #define LIM(x,min,max) MAX(min,MIN(x,max))
459 #define ULIM(x,y,z) ((y) < (z) ? LIM(x,y,z) : LIM(x,z,y))
460
461 static const float xyz_rgb[3][3] = { /* XYZ from RGB */
462 { 0.412453f, 0.357580f, 0.180423f },
463 { 0.212671f, 0.715160f, 0.072169f },
464 { 0.019334f, 0.119193f, 0.950227f } };
465 /* TODO: is it wise to use a D65 here? */
466 static const float d65_white[3] = { 0.950456f, 1.0f, 1.088754f };
467 /* TODO: store the precomputation of xyz_rgb * d65 instead of running a silly init */
468
cam_to_cielab(uint16_t cam[3],float lab[3])469 static void cam_to_cielab(uint16_t cam[3], float lab[3]) /* [SA] */
470 {
471 float xyz[3];
472 static float cbrt[0x10000], xyz_cam[3][4];
473
474 if (cam == NULL) {
475 int i, j;
476
477 for (i = 0; i < 0x10000; i++) {
478 float r = i / 65535.0;
479 cbrt[i] = r > 0.008856 ? pow(r, 1 / 3.0) : 7.787 * r + 16 / 116.0;
480 }
481 for (i = 0; i < 3; i++)
482 for (j = 0; j < 3; j++) /* [SA] */
483 xyz_cam[i][j] = xyz_rgb[i][j] / d65_white[i]; /* [SA] */
484 } else {
485 int c;
486
487 xyz[0] = xyz[1] = xyz[2] = 0.5;
488 for (c = 0; c < 3; c++)
489 { /* [SA] */
490 xyz[0] += xyz_cam[0][c] * cam[c];
491 xyz[1] += xyz_cam[1][c] * cam[c];
492 xyz[2] += xyz_cam[2][c] * cam[c];
493 }
494 xyz[0] = cbrt[round_to_WORD(xyz[0])]; /* [SA] */
495 xyz[1] = cbrt[round_to_WORD(xyz[1])]; /* [SA] */
496 xyz[2] = cbrt[round_to_WORD(xyz[2])]; /* [SA] */
497 lab[0] = 116 * xyz[1] - 16;
498 lab[1] = 500 * (xyz[0] - xyz[1]);
499 lab[2] = 200 * (xyz[1] - xyz[2]);
500 }
501 }
502
503 /*
504 Adaptive Homogeneity-Directed interpolation is based on
505 the work of Keigo Hirakawa, Thomas Parks, and Paul Lee.
506 */
507 #define TS 256 /* Tile Size */
508
bayer_AHD(const WORD * bayer,WORD * dst,int sx,int sy,sensor_pattern pattern)509 static int bayer_AHD(const WORD *bayer, WORD *dst, int sx, int sy,
510 sensor_pattern pattern) {
511 int i, j, top, left, row, col, tr, tc, fc, c, d, val, hm[2];
512 /* the following has the same type as the image */
513 uint16_t (*pix)[3], (*rix)[3]; /* [SA] */
514 static const int dir[4] = { -1, 1, -TS, TS };
515 unsigned ldiff[2][4], abdiff[2][4], leps, abeps;
516 float flab[3];
517 uint16_t (*rgb)[TS][TS][3]; /* [SA] */
518 short (*lab)[TS][TS][3];
519 char (*homo)[TS][TS], *buffer;
520 /* start - new code for libdc1394 */
521 uint32_t filters;
522 const int height = sy, width = sx;
523 int x, y;
524
525 if (ahd_inited == FALSE) {
526 /* WARNING: this might not be multi-processor safe */
527 cam_to_cielab(NULL, NULL);
528 ahd_inited = TRUE;
529 }
530
531 switch (pattern) {
532 case BAYER_FILTER_BGGR:
533 filters = 0x16161616;
534 break;
535 case BAYER_FILTER_GRBG:
536 filters = 0x61616161;
537 break;
538 case BAYER_FILTER_RGGB:
539 filters = 0x94949494;
540 break;
541 case BAYER_FILTER_GBRG:
542 filters = 0x49494949;
543 break;
544 default:
545 return -1;
546 }
547
548 /* fill-in destination with known exact values */
549 for (y = 0; y < height; y++) {
550 for (x = 0; x < width; x++) {
551 int channel = FC(y, x, filters);
552 dst[(y * width + x) * 3 + channel] = bayer[y * width + x];
553 }
554 }
555 /* end - new code for libdc1394 */
556
557 /* start - code from border_interpolate(int border) */
558 {
559 int border = 3;
560 unsigned row, col, y, x, f, c, sum[8];
561
562 for (row = 0; row < (unsigned int) height; row++)
563 for (col = 0; col < (unsigned int) width; col++) {
564 if (col == (unsigned int) border && row >= (unsigned int) border
565 && row < (unsigned int) height - (unsigned int) border)
566 col = width - border;
567 memset(sum, 0, sizeof sum);
568 for (y = row - 1; y != row + 2; y++)
569 for (x = col - 1; x != col + 2; x++)
570 if (y < (unsigned int) height
571 && x < (unsigned int) width) {
572 f = FC(y, x, filters);
573 sum[f] += dst[(y * width + x) * 3 + f]; /* [SA] */
574 sum[f + 4]++;
575 }
576 f = FC(row, col, filters);
577 for (c = 0; c < 3; c++)
578 if (c != f && sum[c + 4]) /* [SA] */
579 dst[(row * width + col) * 3 + c] = sum[c] / sum[c + 4]; /* [SA] */
580 }
581 }
582 /* end - code from border_interpolate(int border) */
583
584 buffer = (char *) malloc(26 * TS * TS); /* 1664 kB */
585 /* merror (buffer, "ahd_interpolate()"); */
586 rgb = (uint16_t (*)[TS][TS][3]) buffer; /* [SA] */
587 lab = (short (*)[TS][TS][3]) (buffer + 12 * TS * TS);
588 homo = (char (*)[TS][TS]) (buffer + 24 * TS * TS);
589
590 for (top = 0; top < height; top += TS - 6)
591 for (left = 0; left < width; left += TS - 6) {
592 memset(rgb, 0, 12 * TS * TS);
593
594 /* Interpolate green horizontally and vertically: */
595 for (row = ((top < 2) ? 2 : top);
596 row < top + TS && row < height - 2; row++) {
597 col = left + (FC(row,left, filters) == 1);
598 if (col < 2)
599 col += 2;
600 for (fc = FC(row, col, filters); col < left + TS && col < width - 2;
601 col += 2) {
602 pix = (uint16_t (*)[3]) dst + (row * width + col); /* [SA] */
603 val = ((pix[-1][1] + pix[0][fc] + pix[1][1]) * 2
604 - pix[-2][fc] - pix[2][fc]) >> 2;
605 rgb[0][row - top][col - left][1] = ULIM(val, pix[-1][1],
606 pix[1][1]);
607 val = ((pix[-width][1] + pix[0][fc] + pix[width][1]) * 2
608 - pix[-2 * width][fc] - pix[2 * width][fc]) >> 2;
609 rgb[1][row - top][col - left][1] = ULIM(val, pix[-width][1],
610 pix[width][1]);
611 }
612 }
613 /* Interpolate red and blue, and convert to CIELab: */
614 for (d = 0; d < 2; d++)
615 for (row = top + 1; row < top + TS - 1 && row < height - 1;
616 row++)
617 for (col = left + 1; col < left + TS - 1 && col < width - 1;
618 col++) {
619 pix = (uint16_t (*)[3]) dst + (row * width + col); /* [SA] */
620 rix = &rgb[d][row - top][col - left];
621 if ((c = 2 - FC(row, col, filters)) == 1) {
622 c = FC(row + 1, col, filters);
623 val = pix[0][1]
624 + ((pix[-1][2 - c] + pix[1][2 - c]
625 - rix[-1][1] - rix[1][1]) >> 1);
626 rix[0][2 - c] = round_to_WORD(val); /* [SA] */
627 val = pix[0][1]
628 + ((pix[-width][c] + pix[width][c]
629 - rix[-TS][1] - rix[TS][1]) >> 1);
630 } else
631 val = rix[0][1]
632 + ((pix[-width - 1][c] + pix[-width + 1][c]
633 + pix[+width - 1][c]
634 + pix[+width + 1][c]
635 - rix[-TS - 1][1] - rix[-TS + 1][1]
636 - rix[+TS - 1][1] - rix[+TS + 1][1]
637 + 1) >> 2);
638 rix[0][c] = round_to_WORD((double) val); /* [SA] */
639 c = FC(row, col, filters);
640 rix[0][c] = pix[0][c];
641 cam_to_cielab(rix[0], flab);
642 for (c = 0; c < 3; c++)
643 lab[d][row - top][col - left][c] = 64 * flab[c];
644 }
645 /* Build homogeneity maps from the CIELab images: */
646 memset(homo, 0, 2 * TS * TS);
647 for (row = top + 2; row < top + TS - 2 && row < height; row++) {
648 tr = row - top;
649 for (col = left + 2; col < left + TS - 2 && col < width;
650 col++) {
651 tc = col - left;
652 for (d = 0; d < 2; d++)
653 for (i = 0; i < 4; i++)
654 ldiff[d][i] = ABSOLU(
655 lab[d][tr][tc][0]
656 - lab[d][tr][tc + dir[i]][0]);
657 leps = MIN(MAX(ldiff[0][0],ldiff[0][1]),
658 MAX(ldiff[1][2],ldiff[1][3]));
659 for (d = 0; d < 2; d++)
660 for (i = 0; i < 4; i++)
661 if (i >> 1 == d || ldiff[d][i] <= leps)
662 abdiff[d][i] =
663 SQR(
664 lab[d][tr][tc][1]
665 - lab[d][tr][tc + dir[i]][1]) + SQR(lab[d][tr][tc][2] -lab[d][tr][tc+dir[i]][2]);
666 abeps = MIN(MAX(abdiff[0][0],abdiff[0][1]),
667 MAX(abdiff[1][2],abdiff[1][3]));
668 for (d = 0; d < 2; d++)
669 for (i = 0; i < 4; i++)
670 if (ldiff[d][i] <= leps && abdiff[d][i] <= abeps)
671 homo[d][tr][tc]++;
672 }
673 }
674 /* Combine the most homogenous pixels for the final result: */
675 for (row = top + 3; row < top + TS - 3 && row < height - 3; row++) {
676 tr = row - top;
677 for (col = left + 3; col < left + TS - 3 && col < width - 3;
678 col++) {
679 tc = col - left;
680 for (d = 0; d < 2; d++)
681 for (hm[d] = 0, i = tr - 1; i <= tr + 1; i++)
682 for (j = tc - 1; j <= tc + 1; j++)
683 hm[d] += homo[d][i][j];
684 if (hm[0] != hm[1])
685 for (c = 0; c < 3; c++)
686 dst[(row * width + col) * 3 + c] = round_to_WORD(
687 rgb[hm[1] > hm[0]][tr][tc][c]); /* [SA] */
688 else
689 for (c = 0; c < 3; c++)
690 dst[(row * width + col) * 3 + c] = round_to_WORD(
691 (rgb[0][tr][tc][c] + rgb[1][tr][tc][c])
692 >> 1); /* [SA] */
693 }
694 }
695 }
696 free(buffer);
697
698 return 0;
699 }
700
701 #define fcol(row, col) xtrans[(row) % 6][(col) % 6]
702
703 /* Code from RAWTherapee:
704 * It is a simple algorithm. Certainly not the best (probably the worst) but it works yet. */
fast_xtrans_interpolate(const WORD * bayer,WORD * dst,int sx,int sy,unsigned int xtrans[6][6])705 static int fast_xtrans_interpolate(const WORD *bayer, WORD *dst, int sx, int sy,
706 unsigned int xtrans[6][6]) {
707 uint32_t filters = 9;
708 const int height = sy, width = sx;
709 int row;
710
711 /* start - code from border_interpolate(int border) */
712 {
713 int border = 1;
714 unsigned row, col, y, x, f, c, sum[8];
715
716 for (row = 0; row < (unsigned int) height; row++)
717 for (col = 0; col < (unsigned int) width; col++) {
718 if (col == (unsigned int) border && row >= (unsigned int) border
719 && row < (unsigned int) height - (unsigned int) border)
720 col = width - border;
721 memset(sum, 0, sizeof sum);
722 for (y = row - 1; y != row + 2; y++)
723 for (x = col - 1; x != col + 2; x++)
724 if (y < (unsigned int) height
725 && x < (unsigned int) width) {
726 f = FC(y, x, filters);
727 sum[f] += dst[(y * width + x) * 3 + f]; /* [SA] */
728 sum[f + 4]++;
729 }
730 f = FC(row, col, filters);
731 for (c = 0; c < 3; c++)
732 if (c != f && sum[c + 4]) /* [SA] */
733 dst[(row * width + col) * 3 + c] = sum[c] / sum[c + 4]; /* [SA] */
734 }
735 }
736 /* end - code from border_interpolate(int border) */
737 #ifdef _OPENMP
738 #pragma omp parallel for num_threads(com.max_thread) schedule(static)
739 #endif
740 for (row = 1; row < (height - 1); row ++) {
741 int col;
742 for (col = 1; col < (width - 1); col ++) {
743 float sum[3] = { 0.f };
744
745 int v, h;
746 for (v = -1; v <= 1; v++) {
747 for (h = -1; h <= 1; h++) {
748 sum[fcol(row + v, col + h)] += bayer[(col + h) + (row + v) * width];
749 }
750 }
751
752 switch (fcol(row, col)) {
753 case 0: /* Red */
754 dst[(row * width + col) * 3 + 0] = bayer[col + row * width];
755 dst[(row * width + col) * 3 + 1] = sum[1] * 0.2f;
756 dst[(row * width + col) * 3 + 2] = sum[2] * 0.33333333f;
757 break;
758
759 case 1: /* Green */
760 dst[(row * width + col) * 3 + 0] = sum[0] * 0.5f;
761 dst[(row * width + col) * 3 + 1] = bayer[col + row * width];
762 dst[(row * width + col) * 3 + 2] = sum[2] * 0.5f;
763 break;
764
765 case 2: /* Blue */
766 dst[(row * width + col) * 3 + 0] = sum[0] * 0.33333333f;
767 dst[(row * width + col) * 3 + 1] = sum[1] * 0.2f;
768 dst[(row * width + col) * 3 + 2] = bayer[col + row * width];
769 break;
770 }
771
772 }
773 }
774 return 0;
775 }
776
777 #undef fcol
778
debayer_buffer_siril(WORD * buf,int * width,int * height,interpolation_method interpolation,sensor_pattern pattern,unsigned int xtrans[6][6],int bit_depth)779 static WORD *debayer_buffer_siril(WORD *buf, int *width, int *height,
780 interpolation_method interpolation, sensor_pattern pattern, unsigned int xtrans[6][6], int bit_depth) {
781 WORD *newbuf;
782 size_t npixels;
783 int retval;
784 switch (interpolation) {
785 default:
786 npixels = (*width) * (*height);
787 break;
788 case BAYER_SUPER_PIXEL:
789 npixels = (*width / 2 + *width % 2) * (*height / 2 + *height % 2);
790 break;
791 }
792 newbuf = calloc(3, npixels * sizeof(WORD));
793 if (!newbuf) {
794 PRINT_ALLOC_ERR;
795 return NULL;
796 }
797
798 switch (interpolation) {
799 case BAYER_BILINEAR:
800 retval = bayer_Bilinear(buf, newbuf, *width, *height, pattern);
801 break;
802 // case BAYER_NEARESTNEIGHBOR:
803 // retval = bayer_NearestNeighbor(buf, newbuf, *width, *height, pattern);
804 // break;
805 default:
806 case BAYER_VNG:
807 retval = bayer_VNG(buf, newbuf, *width, *height, pattern);
808 break;
809 case BAYER_AHD:
810 retval = bayer_AHD(buf, newbuf, *width, *height, pattern);
811 break;
812 case BAYER_SUPER_PIXEL:
813 super_pixel_ushort(buf, newbuf, *width, *height, pattern);
814 *width = *width / 2 + *width % 2;
815 *height = *height / 2 + *height % 2;
816 retval = 0;
817 break;
818 case XTRANS:
819 if (!xtrans)
820 return NULL;
821 retval = fast_xtrans_interpolate(buf, newbuf, *width, *height, xtrans);
822 break;
823 }
824 if (retval) {
825 free(newbuf);
826 return NULL;
827 }
828 return newbuf;
829 }
830
retrieve_Bayer_pattern(fits * fit,sensor_pattern * pattern)831 int retrieve_Bayer_pattern(fits *fit, sensor_pattern *pattern) {
832 int xbayeroff = 0, ybayeroff = 0;
833 gboolean top_down = com.pref.debayer.top_down;
834
835 if (com.pref.debayer.use_bayer_header) {
836 if (!g_strcmp0(fit->row_order, "TOP-DOWN")) {
837 top_down = TRUE;
838 } else if (!g_strcmp0(fit->row_order, "BOTTOM-UP")) {
839 top_down = FALSE;
840 }
841 }
842
843 siril_debug_print("Debayer will be done %s\n", top_down ? "top-down" : "bottom-up");
844
845 if (!com.pref.debayer.use_bayer_header) {
846 xbayeroff = com.pref.debayer.xbayeroff;
847 ybayeroff = com.pref.debayer.ybayeroff;
848 } else {
849 xbayeroff = fit->bayer_xoffset;
850 ybayeroff = fit->bayer_yoffset;
851 }
852
853 if (xbayeroff == 1) {
854 switch (*pattern) {
855 case BAYER_FILTER_RGGB:
856 *pattern = BAYER_FILTER_GRBG;
857 break;
858 case BAYER_FILTER_BGGR:
859 *pattern = BAYER_FILTER_GBRG;
860 break;
861 case BAYER_FILTER_GBRG:
862 *pattern = BAYER_FILTER_BGGR;
863 break;
864 case BAYER_FILTER_GRBG:
865 *pattern = BAYER_FILTER_RGGB;
866 break;
867 default:
868 return 1;
869 }
870 }
871
872 if (ybayeroff == 1) {
873 switch (*pattern) {
874 case BAYER_FILTER_RGGB:
875 *pattern = BAYER_FILTER_GBRG;
876 break;
877 case BAYER_FILTER_BGGR:
878 *pattern = BAYER_FILTER_GRBG;
879 break;
880 case BAYER_FILTER_GBRG:
881 *pattern = BAYER_FILTER_RGGB;
882 break;
883 case BAYER_FILTER_GRBG:
884 *pattern = BAYER_FILTER_BGGR;
885 break;
886 default:
887 return 1;
888 }
889 }
890
891 /* read bottom-up */
892 if (!top_down && !(fit->ry % 2)) {
893 switch (*pattern) {
894 case BAYER_FILTER_RGGB:
895 *pattern = BAYER_FILTER_GBRG;
896 break;
897 case BAYER_FILTER_BGGR:
898 *pattern = BAYER_FILTER_GRBG;
899 break;
900 case BAYER_FILTER_GBRG:
901 *pattern = BAYER_FILTER_RGGB;
902 break;
903 case BAYER_FILTER_GRBG:
904 *pattern = BAYER_FILTER_BGGR;
905 break;
906 default:
907 return 1;
908 }
909 }
910 return 0;
911 }
912
913 /* This function retrieve the xtrans matrix from the FITS header */
retrieve_XTRANS_pattern(char * bayer,unsigned int xtrans[6][6])914 static int retrieve_XTRANS_pattern(char *bayer, unsigned int xtrans[6][6]) {
915 int i = 0;
916
917 if (strlen(bayer) != 36) {
918 siril_log_color_message(_("FITS header does not contain a proper XTRANS pattern, demosaicing cannot be done"), "red");
919 return 1;
920 }
921
922 for (int x = 0; x < 6; x++) {
923 for (int y = 0; y < 6; y++) {
924 switch (bayer[i]) {
925 case 'R':
926 xtrans[x][y] = 0;
927 break;
928 case 'G':
929 xtrans[x][y] = 1;
930 break;
931 case 'B':
932 xtrans[x][y] = 2;
933 break;
934 default:
935 return 1;
936 }
937 i++;
938 }
939 }
940 return 0;
941 }
942
debayer_buffer_superpixel_ushort(WORD * buf,int * width,int * height,sensor_pattern pattern)943 WORD *debayer_buffer_superpixel_ushort(WORD *buf, int *width, int *height, sensor_pattern pattern) {
944 int new_rx = *width / 2 + *width % 2;
945 int new_ry = *height / 2 + *height % 2;
946 size_t npixels = new_rx * new_ry;
947 WORD *newbuf = malloc(3 * npixels * sizeof(WORD));
948 if (!newbuf) {
949 PRINT_ALLOC_ERR;
950 return NULL;
951 }
952 super_pixel_ushort(buf, newbuf, *width, *height, pattern);
953 *width = new_rx;
954 *height = new_ry;
955 return newbuf;
956 }
957
958 /**
959 * debayer a buffer of a given size into a newly allocated and returned buffer,
960 * using the given bayer pattern and interpolation (only used for SER demosaicing)
961 *
962 * @param buf original RAW data
963 * @param width width of image
964 * @param height height of image
965 * @param interpolation type of interpolation used for demosaicing algorithm
966 * @param pattern type of pattern used for demosaicing algorithm
967 * @return a new buffer of demosaiced data
968 */
debayer_buffer(WORD * buf,int * width,int * height,interpolation_method interpolation,sensor_pattern pattern,int bit_depth)969 WORD *debayer_buffer(WORD *buf, int *width, int *height,
970 interpolation_method interpolation, sensor_pattern pattern, int bit_depth) {
971 if (USE_SIRIL_DEBAYER)
972 return debayer_buffer_siril(buf, width, height, interpolation, pattern, NULL, bit_depth);
973 return debayer_buffer_new_ushort(buf, width, height, interpolation, pattern, NULL, bit_depth);
974 }
975
debayer_buffer_superpixel_float(float * buf,int * width,int * height,sensor_pattern pattern)976 float *debayer_buffer_superpixel_float(float *buf, int *width, int *height, sensor_pattern pattern) {
977 int new_rx = *width / 2 + *width % 2;
978 int new_ry = *height / 2 + *height % 2;
979 size_t npixels = new_rx * new_ry;
980 float *newbuf = malloc(3 * npixels * sizeof(float));
981 if (!newbuf) {
982 PRINT_ALLOC_ERR;
983 return NULL;
984 }
985 super_pixel_float(buf, newbuf, *width, *height, pattern);
986 *width = new_rx;
987 *height = new_ry;
988 return newbuf;
989 }
990
991 /* From an area, get the area corresponding to the debayer data for all colors,
992 * the dashed area below.
993 * 0 1 2 3 4 5
994 * - - - - - -
995 * - - - - - -
996 * - - G R - -
997 * - - B G - -
998 * - - - - - -
999 * - - - - - -
1000 *
1001 * area is the requested area of an image (simplified as GRBG above)
1002 * debayer_area is the result of this function, the area with enough pixels to
1003 * have a valid debayer
1004 * image_area is the size of the image, to avoid going out of bounds
1005 * debayer_offset_x and y are the offset that need to be applied to the debayer
1006 * data to find the original area (between 0 and 3).
1007 */
get_debayer_area(const rectangle * area,rectangle * debayer_area,const rectangle * image_area,int * debayer_offset_x,int * debayer_offset_y)1008 void get_debayer_area(const rectangle *area, rectangle *debayer_area,
1009 const rectangle *image_area, int *debayer_offset_x,
1010 int *debayer_offset_y) {
1011 int right, bottom; // temp debayer negative offsets
1012
1013 /* left side */
1014 if (area->x & 1)
1015 *debayer_offset_x = 11;
1016 else
1017 *debayer_offset_x = 10;
1018 if (area->x - *debayer_offset_x < 0) {
1019 debayer_area->x = 0;
1020 *debayer_offset_x = area->x;
1021 } else {
1022 debayer_area->x = area->x - *debayer_offset_x;
1023 }
1024
1025 /* right side */
1026 int xend = area->x + area->w - 1;
1027 if (xend & 1)
1028 right = 10;
1029 else
1030 right = 11;
1031 if (xend + right >= image_area->w) {
1032 right = image_area->w - xend - 1;
1033 }
1034 debayer_area->w = area->w + (area->x - debayer_area->x) + right;
1035
1036 /* top */
1037 if (area->y & 1)
1038 *debayer_offset_y = 11;
1039 else
1040 *debayer_offset_y = 10;
1041 if (area->y - *debayer_offset_y < 0) {
1042 debayer_area->y = 0;
1043 *debayer_offset_y = area->y;
1044 } else {
1045 debayer_area->y = area->y - *debayer_offset_y;
1046 }
1047
1048 /* bottom */
1049 int yend = area->y + area->h - 1;
1050 if (yend & 1)
1051 bottom = 10;
1052 else
1053 bottom = 11;
1054 if (yend + bottom >= image_area->h) {
1055 bottom = image_area->h - yend - 1;
1056 }
1057 debayer_area->h = area->h + (area->y - debayer_area->y) + bottom;
1058
1059 assert(debayer_area->x < image_area->w);
1060 assert(debayer_area->y < image_area->h);
1061 assert(debayer_area->h > 2);
1062 assert(debayer_area->w > 2);
1063 }
1064
clear_Bayer_information(fits * fit)1065 void clear_Bayer_information(fits *fit) {
1066 memset(fit->bayer_pattern, 0, FLEN_VALUE);
1067 }
1068
debayer_ushort(fits * fit,interpolation_method interpolation,sensor_pattern pattern)1069 static int debayer_ushort(fits *fit, interpolation_method interpolation, sensor_pattern pattern) {
1070 size_t i, j, npixels = fit->naxes[0] * fit->naxes[1];
1071 int width = fit->rx;
1072 int height = fit->ry;
1073 WORD *buf = fit->data;
1074 gboolean top_down = com.pref.debayer.top_down;
1075
1076 unsigned int xtrans[6][6];
1077 if (interpolation == XTRANS) {
1078 if (com.pref.debayer.use_bayer_header) {
1079 if (!g_strcmp0(fit->row_order, "TOP-DOWN")) {
1080 top_down = TRUE;
1081 } else if (!g_strcmp0(fit->row_order, "BOTTOM-UP")) {
1082 top_down = FALSE;
1083 }
1084 }
1085 if (!top_down)
1086 fits_flip_top_to_bottom(fit); // TODO: kind of ugly but not easy with xtrans
1087 retrieve_XTRANS_pattern(fit->bayer_pattern, xtrans);
1088 } else {
1089 retrieve_Bayer_pattern(fit, &pattern);
1090 }
1091
1092 if (USE_SIRIL_DEBAYER) {
1093 WORD *newbuf = debayer_buffer_siril(buf, &width, &height, interpolation, pattern, xtrans, fit->bitpix);
1094 if (!newbuf)
1095 return 1;
1096
1097 fit_debayer_buffer(fit, newbuf);
1098 // size might have changed, in case of superpixel
1099 fit->naxes[0] = width;
1100 fit->naxes[1] = height;
1101 fit->rx = width;
1102 fit->ry = height;
1103 fit->bitpix = fit->orig_bitpix;
1104 // color RGBRGB format to fits RRGGBB format
1105 for (i = 0, j = 0; j < npixels; i += 3, j++) {
1106 double r = (double) newbuf[i + RLAYER];
1107 double g = (double) newbuf[i + GLAYER];
1108 double b = (double) newbuf[i + BLAYER];
1109 fit->pdata[RLAYER][j] =
1110 (fit->bitpix == 8) ? round_to_BYTE(r) : round_to_WORD(r);
1111 fit->pdata[GLAYER][j] =
1112 (fit->bitpix == 8) ? round_to_BYTE(g) : round_to_WORD(g);
1113 fit->pdata[BLAYER][j] =
1114 (fit->bitpix == 8) ? round_to_BYTE(b) : round_to_WORD(b);
1115 }
1116 } else {
1117 // use librtprocess debayer
1118 WORD *newbuf = debayer_buffer_new_ushort(buf, &width, &height, interpolation, pattern, xtrans, fit->bitpix);
1119 if (!newbuf)
1120 return 1;
1121
1122 fit_debayer_buffer(fit, newbuf);
1123 if (interpolation == XTRANS && !top_down) {
1124 fits_flip_top_to_bottom(fit);
1125 }
1126 }
1127 /* we remove Bayer header because not needed now */
1128 clear_Bayer_information(fit);
1129 return 0;
1130 }
1131
debayer_float(fits * fit,interpolation_method interpolation,sensor_pattern pattern)1132 static int debayer_float(fits* fit, interpolation_method interpolation, sensor_pattern pattern) {
1133 int width = fit->rx;
1134 int height = fit->ry;
1135 float *buf = fit->fdata;
1136 gboolean top_down = com.pref.debayer.top_down;
1137
1138 unsigned int xtrans[6][6];
1139 if (interpolation == XTRANS) {
1140 if (com.pref.debayer.use_bayer_header) {
1141 if (!g_strcmp0(fit->row_order, "TOP-DOWN")) {
1142 top_down = TRUE;
1143 } else if (!g_strcmp0(fit->row_order, "BOTTOM-UP")) {
1144 top_down = FALSE;
1145 }
1146 }
1147 if (!top_down)
1148 fits_flip_top_to_bottom(fit); // TODO: kind of ugly but not easy with xtrans
1149 retrieve_XTRANS_pattern(fit->bayer_pattern, xtrans);
1150 } else {
1151 retrieve_Bayer_pattern(fit, &pattern);
1152 }
1153
1154 float *newbuf = debayer_buffer_new_float(buf, &width, &height, interpolation, pattern, xtrans);
1155 if (!newbuf)
1156 return 1;
1157
1158 fit_debayer_buffer(fit, newbuf);
1159 if (interpolation == XTRANS && !top_down) {
1160 fits_flip_top_to_bottom(fit);
1161 }
1162 /* we remove Bayer header because not needed now */
1163 clear_Bayer_information(fit);
1164 return 0;
1165 }
1166
debayer(fits * fit,interpolation_method interpolation,sensor_pattern pattern)1167 int debayer(fits *fit, interpolation_method interpolation, sensor_pattern pattern) {
1168 if (fit->type == DATA_USHORT)
1169 return debayer_ushort(fit, interpolation, pattern);
1170 else if (fit->type == DATA_FLOAT)
1171 return debayer_float(fit, interpolation, pattern);
1172 else return -1;
1173 }
1174
retrieveBayerPatternFromChar(char * bayer)1175 int retrieveBayerPatternFromChar(char *bayer) {
1176 for (int i = 0; i < G_N_ELEMENTS(filter_pattern); i++) {
1177 if (g_ascii_strcasecmp(bayer, filter_pattern[i]) == 0) {
1178 return i;
1179 }
1180 }
1181 return BAYER_FILTER_NONE;
1182 }
1183
1184 // debayers the image if it's a FITS image and if debayer is activated globally
1185 // or if the force argument is passed
debayer_if_needed(image_type imagetype,fits * fit,gboolean force_debayer)1186 int debayer_if_needed(image_type imagetype, fits *fit, gboolean force_debayer) {
1187 if (imagetype != TYPEFITS || (!com.pref.debayer.open_debayer && !force_debayer))
1188 return 0;
1189
1190 /* Siril's FITS are stored bottom-up, debayering will give wrong results.
1191 * So before demosacaing we need to flip the image with fits_flip_top_to_bottom().
1192 * But sometimes FITS are created by acquisition software top-down, in that case
1193 * the user can indicate it ('compatibility') and we don't flip for debayer.
1194 */
1195 if (fit->naxes[2] != 1) {
1196 siril_log_message(_("Cannot perform debayering on image with more than one channel\n"));
1197 return 0;
1198 }
1199
1200 /* Get Bayer informations from header if available */
1201 sensor_pattern tmp_pattern = com.pref.debayer.bayer_pattern;
1202 interpolation_method tmp_algo = com.pref.debayer.bayer_inter;
1203 if (com.pref.debayer.use_bayer_header) {
1204 sensor_pattern bayer;
1205 bayer = retrieveBayerPatternFromChar(fit->bayer_pattern);
1206
1207 if (bayer <= BAYER_FILTER_MAX) {
1208 if (bayer != tmp_pattern) {
1209 if (bayer == BAYER_FILTER_NONE) {
1210 siril_log_color_message(_("No Bayer pattern found in the header file.\n"), "red");
1211 }
1212 else {
1213 siril_log_color_message(_("Bayer pattern found in header (%s) is different"
1214 " from Bayer pattern in settings (%s). Overriding settings.\n"),
1215 "salmon", filter_pattern[bayer], filter_pattern[com.pref.debayer.bayer_pattern]);
1216 tmp_pattern = bayer;
1217 }
1218 }
1219 } else {
1220 tmp_pattern = bayer;
1221 tmp_algo = XTRANS;
1222 siril_log_color_message(_("XTRANS Sensor detected. Using special algorithm.\n"), "green");
1223 }
1224 }
1225 if (tmp_pattern >= BAYER_FILTER_MIN && tmp_pattern <= BAYER_FILTER_MAX) {
1226 siril_log_message(_("Filter Pattern: %s\n"),
1227 filter_pattern[tmp_pattern]);
1228 }
1229
1230 int retval = 0;
1231 if (debayer(fit, tmp_algo, tmp_pattern)) {
1232 siril_log_message(_("Cannot perform debayering\n"));
1233 retval = -1;
1234 }
1235 return retval;
1236 }
1237