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 /* Management of Siril's internal image format: unsigned 16-bit FITS */
22 #include <stdio.h>
23 #include <stdlib.h>
24 #include <math.h>
25 #include <ctype.h>
26 #include <string.h>
27 #include <float.h>
28 #include <gsl/gsl_statistics.h>
29 #ifdef _WIN32
30 #include <windows.h>
31 #endif
32
33 #include "core/siril.h"
34 #include "core/proto.h"
35 #include "core/siril_date.h"
36 #include "core/OS_utils.h"
37 #include "io/sequence.h"
38 #include "io/fits_sequence.h"
39 #include "gui/utils.h"
40 #include "gui/progress_and_log.h"
41 #include "algos/statistics.h"
42 #include "io/sequence.h"
43 #include "io/single_image.h"
44 #include "image_format_fits.h"
45 #include "algos/siril_wcs.h"
46
47 const char *fit_extension[] = {
48 ".fit",
49 ".fits",
50 ".fts",
51 ".fits.fz"
52 };
53
54 static char *MIPSHI[] = {"MIPS-HI", "CWHITE", "DATAMAX", NULL };
55 static char *MIPSLO[] = {"MIPS-LO", "CBLACK", "DATAMIN", NULL };
56 static char *PixSizeX[] = { "XPIXSZ", "XPIXELSZ", "PIXSIZE1", "PIXSIZEX", NULL };
57 static char *PixSizeY[] = { "YPIXSZ", "YPIXELSZ", "PIXSIZE2", "PIXSIZEY", NULL };
58 static char *BinX[] = { "XBINNING", "BINX", NULL };
59 static char *BinY[] = { "YBINNING", "BINY", NULL };
60 static char *Focal[] = { "FOCAL", "FOCALLEN", NULL };
61 static char *CCD_TEMP[] = { "CCD-TEMP", "CCD_TEMP", "CCDTEMP", "TEMPERAT", NULL };
62 static char *Exposure[] = { "EXPTIME", "EXPOSURE", NULL };
63 static char *Cvf[] = { "CVF", "EGAIN", NULL };
64 static char *OffsetLevel[] = { "OFFSET", "BLKLEVEL", NULL }; //Used for synthetic offset
65 static int CompressionMethods[] = { RICE_1, GZIP_1, GZIP_2, HCOMPRESS_1};
66
67 #define __tryToFindKeywords(fptr, type, keywords, value) \
68 { \
69 int __iter__ = 0; \
70 int __status__; \
71 do { \
72 __status__ = 0; \
73 fits_read_key(fptr, type, keywords[__iter__], value, NULL, &__status__); \
74 __iter__++; \
75 } while ((keywords[__iter__]) && (__status__ > 0)); \
76 }
77
read_fits_date_obs_header(fits * fit)78 static void read_fits_date_obs_header(fits *fit) {
79 int status = 0;
80 char ut_start[FLEN_VALUE] = { 0 };
81 char date_obs[FLEN_VALUE] = { 0 };
82
83 fits_read_key(fit->fptr, TSTRING, "DATE-OBS", &date_obs, NULL, &status);
84
85 /** Case seen in some FITS files. Needed to get date back in SER conversion **/
86 status = 0;
87 fits_read_key(fit->fptr, TSTRING, "UT-START", &ut_start, NULL, &status);
88 if (status == 0 && ut_start[0] != '\0' && date_obs[2] == '/') {
89 int year, month, day;
90 if (sscanf(date_obs, "%02d/%02d/%04d", &day, &month, &year) == 3) {
91 g_snprintf(date_obs, sizeof(date_obs), "%04d-%02d-%02dT%s", year, month, day, ut_start);
92 }
93 }
94 fit->date_obs = FITS_date_to_date_time(date_obs);
95 }
96
fit_get_photometry_data(fits * fit)97 void fit_get_photometry_data(fits *fit) {
98 read_fits_date_obs_header(fit);
99 __tryToFindKeywords(fit->fptr, TDOUBLE, Exposure, &fit->exposure);
100 }
101
fit_stats(fits * fit,float * mini,float * maxi)102 static int fit_stats(fits *fit, float *mini, float *maxi) {
103 int status = 0;
104 int ii;
105 long npixels = 1L;
106 long anaxes[3] = { 1L, 1L, 1L }, firstpix[3] = { 1L, 1L, 1L };
107 float *pix, sum = 0.f;
108 float meanval = 0.f, minval = 1.E33, maxval = -1.E33;
109
110 /* initialize value in case where it does not work */
111 *mini = 0;
112 *maxi = 0;
113
114 fits_get_img_size(fit->fptr, 3, anaxes, &status);
115
116 if (status) {
117 report_fits_error(status); /* print error message */
118 return(status);
119 }
120
121 npixels = anaxes[0]; /* no. of pixels to read in each row */
122 pix = malloc(npixels * sizeof(float)); /* memory for 1 row */
123 if (pix == NULL) {
124 PRINT_ALLOC_ERR;
125 return (1);
126 }
127
128 /* loop over all planes of the cube (2D images have 1 plane) */
129 for (firstpix[2] = 1; firstpix[2] <= anaxes[2]; firstpix[2]++) {
130 /* loop over all rows of the plane */
131 for (firstpix[1] = 1; firstpix[1] <= anaxes[1]; firstpix[1]++) {
132 /* give starting pixel coordinate and number of pixels to read */
133 if (fits_read_pix(fit->fptr, TFLOAT, firstpix, npixels, NULL, pix,
134 NULL, &status))
135 break; /* jump out of loop on error */
136
137 for (ii = 0; ii < npixels; ii++) {
138 sum += pix[ii]; /* accumulate sum */
139 if (pix[ii] < minval)
140 minval = pix[ii]; /* find min and */
141 if (pix[ii] > maxval)
142 maxval = pix[ii]; /* max values */
143 }
144 }
145 } /* end of loop over planes */
146 free(pix);
147
148 if (status) {
149 report_fits_error(status); /* print any error message */
150 } else {
151 if (npixels > 0)
152 meanval = sum / npixels;
153 siril_debug_print(" sum of pixels = %f\n", sum);
154 siril_debug_print(" mean value = %f\n", meanval);
155 siril_debug_print(" minimum value = %f\n", minval);
156 siril_debug_print(" maximum value = %f\n", maxval);
157 *maxi = maxval;
158 *mini = minval;
159 }
160 return status;
161 }
162
read_history_in_hdu(fitsfile * fptr,GSList ** list)163 static void read_history_in_hdu(fitsfile *fptr, GSList **list) {
164 int nkeys, status = 0;
165 char card[FLEN_CARD];
166 fits_get_hdrspace(fptr, &nkeys, NULL, &status);
167 for (int i = 1; i <= nkeys; i++) {
168 if (fits_read_record(fptr, i, card, &status))
169 break;
170 if (!strncmp(card, "HISTORY", 7)) {
171 *list = g_slist_prepend(*list, g_strdup(card + 8));
172 }
173 }
174 }
175
176 /* copy the header into a list of strings
177 * header is read from current HDU and the following HDU as long as they don't contain an image.
178 * Original active HDU is restored */
fits_read_history(fitsfile * fptr,GSList ** history)179 static void fits_read_history(fitsfile *fptr, GSList **history) {
180 GSList *list = NULL;
181 read_history_in_hdu(fptr, &list); // read from current HDU
182
183 // browse following HDU
184 int orig_hdu, type, status;
185 gboolean hdu_changed = FALSE;
186 fits_get_hdu_num(fptr, &orig_hdu);
187 do {
188 status = 0;
189 fits_movrel_hdu(fptr, 1, &type, &status);
190 if (status)
191 break;
192 hdu_changed = TRUE;
193 if (type == IMAGE_HDU)
194 break;
195 siril_debug_print("history read from another HDU (CHDU changed)\n");
196 read_history_in_hdu(fptr, &list);
197 } while (1);
198
199 // restore
200 if (hdu_changed) {
201 status = 0;
202 fits_movabs_hdu(fptr, orig_hdu, NULL, &status);
203 }
204
205 if (*history)
206 g_slist_free_full(*history, free);
207 list = g_slist_reverse(list);
208 *history = list;
209 }
210
try_read_float_lo_hi(fitsfile * fptr,WORD * lo,WORD * hi)211 static int try_read_float_lo_hi(fitsfile *fptr, WORD *lo, WORD *hi) {
212 float fhi, flo;
213 int status = 0;
214 fits_read_key(fptr, TFLOAT, "MIPS-FHI", &fhi, NULL, &status);
215 if (!status) {
216 *hi = float_to_ushort_range(fhi);
217 status = 0;
218 fits_read_key(fptr, TFLOAT, "MIPS-FLO", &flo, NULL, &status);
219 if (!status) {
220 *lo = float_to_ushort_range(flo);
221 }
222 }
223 return status;
224 }
225
226
227 /* reading the FITS header to get useful information
228 * stored in the fit, requires an opened file descriptor */
read_fits_header(fits * fit)229 void read_fits_header(fits *fit) {
230 /* about the status argument: http://heasarc.gsfc.nasa.gov/fitsio/c/c_user/node28.html */
231 int status = 0;
232 double scale, zero;
233 char str[FLEN_VALUE] = { 0 };
234
235 __tryToFindKeywords(fit->fptr, TUSHORT, MIPSLO, &fit->lo);
236 __tryToFindKeywords(fit->fptr, TUSHORT, MIPSHI, &fit->hi);
237 if (fit->orig_bitpix == FLOAT_IMG || fit->orig_bitpix == DOUBLE_IMG) {
238 try_read_float_lo_hi(fit->fptr, &fit->lo, &fit->hi);
239 }
240
241 if (fit->orig_bitpix == SHORT_IMG) {
242 if (fit->lo)
243 fit->lo += 32768;
244 if (fit->hi)
245 fit->hi += 32768;
246 }
247
248 status = 0;
249 fits_read_key(fit->fptr, TDOUBLE, "BSCALE", &scale, NULL, &status);
250 if (!status && 1.0 != scale) {
251 siril_log_message(_("Loaded FITS file has "
252 "a BSCALE different than 1 (%f)\n"), scale);
253 status = 0;
254 /* We reset the scaling factors as we don't use it */
255 fits_set_bscale(fit->fptr, 1.0, 0.0, &status);
256 }
257
258 status = 0;
259 fits_read_key(fit->fptr, TDOUBLE, "BZERO", &zero, NULL, &status);
260 if (!status && 0.0 != zero && fit->bitpix == FLOAT_IMG) {
261 fprintf(stdout, "ignoring BZERO\n");
262 fits_set_bscale(fit->fptr, 1.0, 0.0, &status);
263 }
264
265 /* check if file is created by Siril. If so, and if the file is FLOAT_IMG
266 * Then we are confident enough that the range is [0, 1]. So, no need to
267 * compute fit_stats
268 */
269 status = 0;
270 fits_read_key(fit->fptr, TSTRING, "PROGRAM", &str, NULL, &status);
271 gboolean not_from_siril = g_ascii_strncasecmp(str, PACKAGE, strlen(PACKAGE));
272
273 if ((fit->bitpix == FLOAT_IMG && not_from_siril) || (fit->bitpix == DOUBLE_IMG)) {
274 status = 0;
275 fits_read_key(fit->fptr, TDOUBLE, "DATAMAX", &(fit->data_max), NULL, &status);
276 if (status == KEY_NO_EXIST) {
277 float mini, maxi;
278 fit_stats(fit, &mini, &maxi);
279 fit->data_max = (double) maxi;
280 }
281 }
282
283 status = 0;
284 fits_read_key(fit->fptr, TSTRING, "ROWORDER", &(fit->row_order), NULL,
285 &status);
286
287 /*******************************************************************
288 * ************* CAMERA AND INSTRUMENT KEYWORDS ********************
289 * ****************************************************************/
290
291 __tryToFindKeywords(fit->fptr, TFLOAT, PixSizeX, &fit->pixel_size_x);
292 __tryToFindKeywords(fit->fptr, TFLOAT, PixSizeY, &fit->pixel_size_y);
293 __tryToFindKeywords(fit->fptr, TUINT, BinX, &fit->binning_x);
294 if (fit->binning_x <= 0) fit->binning_x = 1;
295 __tryToFindKeywords(fit->fptr, TUINT, BinY, &fit->binning_y);
296 if (fit->binning_y <= 0) fit->binning_y = 1;
297
298 status = 0;
299 fits_read_key(fit->fptr, TSTRING, "INSTRUME", &(fit->instrume), NULL,
300 &status);
301
302 status = 0;
303 fits_read_key(fit->fptr, TSTRING, "TELESCOP", &(fit->telescop), NULL,
304 &status);
305
306 status = 0;
307 fits_read_key(fit->fptr, TSTRING, "OBSERVER", &(fit->observer), NULL,
308 &status);
309
310 status = 0;
311 fits_read_key(fit->fptr, TSTRING, "BAYERPAT", &(fit->bayer_pattern), NULL,
312 &status);
313
314 status = 0;
315 fits_read_key(fit->fptr, TINT, "XBAYROFF", &(fit->bayer_xoffset), NULL,
316 &status);
317
318 status = 0;
319 fits_read_key(fit->fptr, TINT, "YBAYROFF", &(fit->bayer_yoffset), NULL,
320 &status);
321
322 read_fits_date_obs_header(fit);
323
324 status = 0;
325 char date[FLEN_VALUE];
326 fits_read_key(fit->fptr, TSTRING, "DATE", &date, NULL, &status);
327 fit->date = FITS_date_to_date_time(date);
328
329 __tryToFindKeywords(fit->fptr, TDOUBLE, Focal, &fit->focal_length);
330 if (fit->focal_length <= 0.0) {
331 /* this keyword is seen in some professional images, FLENGTH is in m. */
332 double flength;
333 status = 0;
334 fits_read_key(fit->fptr, TDOUBLE, "FLENGTH", &flength, NULL, &status);
335 if (!status) {
336 fit->focal_length = flength * 1000.0; // convert m to mm
337 }
338 }
339
340 __tryToFindKeywords(fit->fptr, TDOUBLE, CCD_TEMP, &fit->ccd_temp);
341 __tryToFindKeywords(fit->fptr, TDOUBLE, Exposure, &fit->exposure);
342
343 status = 0;
344 fits_read_key(fit->fptr, TDOUBLE, "APERTURE", &(fit->aperture), NULL, &status);
345
346 status = 0;
347 fits_read_key(fit->fptr, TDOUBLE, "ISOSPEED", &(fit->iso_speed), NULL, &status); // Non-standard keywords used in MaxIm DL
348
349 __tryToFindKeywords(fit->fptr, TDOUBLE, Cvf, &fit->cvf); // conversion gain in e-/ADU
350
351 status = 0;
352 fits_read_key(fit->fptr, TUSHORT, "GAIN", &(fit->key_gain), NULL, &status); // Gain setting from camera
353
354 __tryToFindKeywords(fit->fptr, TUSHORT, OffsetLevel, &fit->key_offset); // Offset setting from camera
355 /*******************************************************************
356 * ******************* PLATE SOLVING KEYWORDS **********************
357 * ****************************************************************/
358 status = 0;
359 fits_read_key(fit->fptr, TDOUBLE, "EQUINOX", &(fit->wcsdata.equinox), NULL, &status);
360
361 status = 0;
362 fits_read_key(fit->fptr, TSTRING, "OBJCTRA", &(fit->wcsdata.objctra), NULL, &status);
363
364 status = 0;
365 fits_read_key(fit->fptr, TSTRING, "OBJCTDEC", &(fit->wcsdata.objctdec), NULL, &status);
366
367 status = 0;
368 fits_read_key(fit->fptr, TDOUBLE, "CRPIX1", &(fit->wcsdata.crpix[0]), NULL, &status);
369
370 status = 0;
371 fits_read_key(fit->fptr, TDOUBLE, "CRPIX2", &(fit->wcsdata.crpix[1]), NULL, &status);
372
373 status = 0;
374 fits_read_key(fit->fptr, TDOUBLE, "CRVAL1", &(fit->wcsdata.crval[0]), NULL, &status);
375
376 status = 0;
377 fits_read_key(fit->fptr, TDOUBLE, "CRVAL2", &(fit->wcsdata.crval[1]), NULL, &status);
378
379 status = 0;
380 fits_read_key(fit->fptr, TDOUBLE, "CD1_1", &(fit->wcsdata.cd[0][0]), NULL, &status);
381
382 status = 0;
383 fits_read_key(fit->fptr, TDOUBLE, "CD1_2", &(fit->wcsdata.cd[0][1]), NULL, &status);
384
385 status = 0;
386 fits_read_key(fit->fptr, TDOUBLE, "CD2_1", &(fit->wcsdata.cd[1][0]), NULL, &status);
387
388 status = 0;
389 fits_read_key(fit->fptr, TDOUBLE, "CD2_2", &(fit->wcsdata.cd[1][1]), NULL, &status);
390
391 status = 0;
392 fits_read_key(fit->fptr, TDOUBLE, "CDELT1", &(fit->wcsdata.cdelt[0]), NULL, &status);
393
394 status = 0;
395 fits_read_key(fit->fptr, TDOUBLE, "CDELT2", &(fit->wcsdata.cdelt[1]), NULL, &status);
396
397 status = 0;
398 fits_read_key(fit->fptr, TDOUBLE, "CROTA1", &(fit->wcsdata.crota[0]), NULL, &status);
399
400 status = 0;
401 fits_read_key(fit->fptr, TDOUBLE, "CROTA2", &(fit->wcsdata.crota[1]), NULL, &status);
402
403 load_WCS_from_file(fit);
404
405 /*******************************************************************
406 * ************************* DFT KEYWORDS **************************
407 * ****************************************************************/
408
409 status = 0;
410 fits_read_key(fit->fptr, TDOUBLE, "DFTNORM0", &(fit->dft.norm[0]), NULL, &status);
411 status = 0;
412 fits_read_key(fit->fptr, TDOUBLE, "DFTNORM1", &(fit->dft.norm[1]), NULL, &status);
413 status = 0;
414 fits_read_key(fit->fptr, TDOUBLE, "DFTNORM2", &(fit->dft.norm[2]), NULL, &status);
415 status = 0;
416 fits_read_key(fit->fptr, TSTRING, "DFTORD", &(fit->dft.ord), NULL, &status);
417 status = 0;
418 fits_read_key(fit->fptr, TSTRING, "DFTTYPE", &(fit->dft.type), NULL, &status);
419
420 fits_read_history(fit->fptr, &(fit->history));
421 }
422
423 /* copy the header for the current HDU in a heap-allocated string */
copy_header_from_hdu(fitsfile * fptr,char ** header,int * strsize,int * strlength)424 static int copy_header_from_hdu(fitsfile *fptr, char **header, int *strsize, int *strlength) {
425 int nkeys, status = 0;
426 fits_get_hdrspace(fptr, &nkeys, NULL, &status);
427 if (status || nkeys < 0) {
428 free(*header);
429 return 1;
430 }
431 for (int i = 1; i <= nkeys; i++) {
432 int cardlen;
433 char *newstr;
434 char card[FLEN_CARD];
435 if (fits_read_record(fptr, i, card, &status))
436 break;
437 cardlen = strlen(card);
438 if (*strlength + cardlen + 1 >= *strsize) {
439 *strsize += 567;
440 newstr = realloc(*header, *strsize);
441 if (!newstr) {
442 PRINT_ALLOC_ERR;
443 free(*header);
444 return 1;
445 }
446 *header = newstr;
447 }
448 strcpy(*header + *strlength, card);
449 (*strlength) += cardlen;
450 strcpy(*header + *strlength, "\n");
451 (*strlength)++;
452 }
453 return 0;
454 }
455
456 /* copy the complete header in a heap-allocated string */
copy_header(fits * fit)457 char *copy_header(fits *fit) {
458 int strsize, strlength;
459 char *header;
460
461 /* each line in the FITS header is 80 character wide
462 * in our string, we also keep new lines, so that's 81.
463 * initial allocation is 20 times 81 = 1620
464 * reallocations are 7 lines, 567 */
465 strsize = 1620;
466 if (!(header = malloc(strsize))) {
467 PRINT_ALLOC_ERR;
468 return NULL;
469 }
470 header[0] = '\0';
471 strlength = 0;
472 if (copy_header_from_hdu(fit->fptr, &header, &strsize, &strlength))
473 return NULL;
474
475 int orig_hdu, type, status;
476 gboolean hdu_changed = FALSE;
477 fits_get_hdu_num(fit->fptr, &orig_hdu);
478 do {
479 status = 0;
480 fits_movrel_hdu(fit->fptr, 1, &type, &status);
481 if (status)
482 break;
483 hdu_changed = TRUE;
484 if (type == IMAGE_HDU)
485 break;
486 siril_debug_print("header read from another HDU (CHDU changed)\n");
487 if (copy_header_from_hdu(fit->fptr, &header, &strsize, &strlength))
488 break;
489 } while (1);
490
491 // restore
492 if (hdu_changed) {
493 status = 0;
494 fits_movabs_hdu(fit->fptr, orig_hdu, NULL, &status);
495 }
496
497 if (header[0] == '\0') {
498 free(header);
499 header = NULL;
500 }
501 if (!header)
502 return NULL;
503
504 /* we need to test if text is ut8
505 * indeed some header are not */
506 if (!g_utf8_validate(header, -1, NULL)) {
507 gchar *str = g_utf8_make_valid(header, -1);
508 free(header);
509 header = strdup(str);
510 g_free(str);
511 }
512
513 return header;
514 }
515
516 /***** data reading and transformation ******/
517
get_data_type(int bitpix)518 data_type get_data_type(int bitpix) {
519 if (bitpix == BYTE_IMG || bitpix == SHORT_IMG || bitpix == USHORT_IMG)
520 return DATA_USHORT;
521 if (bitpix == LONGLONG_IMG)
522 return DATA_UNSUPPORTED;
523 return DATA_FLOAT;
524 }
525
526
report_fits_error(int status)527 void report_fits_error(int status) {
528 if (status) {
529 char errmsg[FLEN_ERRMSG];
530 while (fits_read_errmsg(errmsg)) {
531 siril_log_message(_("FITS error: %s\n"),
532 errmsg[0] != '\0' ? errmsg : "unknown" );
533 }
534 }
535 }
536
conv_8_to_16(WORD * data,size_t nbdata)537 static void conv_8_to_16(WORD *data, size_t nbdata) {
538 for (size_t i = 0; i < nbdata; i++) {
539 double tmp = (double) data[i] / UCHAR_MAX_DOUBLE * USHRT_MAX_DOUBLE;
540 data[i] = round_to_WORD(tmp);
541 }
542 }
543
conv_16_to_32(WORD * udata,float * fdata,size_t nbdata)544 static void conv_16_to_32(WORD *udata, float *fdata, size_t nbdata) {
545 for (size_t i = 0; i < nbdata; i++) {
546 fdata[i] = (double) udata[i] / USHRT_MAX_DOUBLE;
547 }
548 }
549
550 /* convert FITS data formats to siril native.
551 * nbdata is the number of pixels, w * h.
552 * from is not freed, to must be allocated and can be the same as from */
convert_data_ushort(int bitpix,const void * from,WORD * to,size_t nbdata,gboolean values_above_1)553 static void convert_data_ushort(int bitpix, const void *from, WORD *to, size_t nbdata, gboolean values_above_1) {
554 BYTE *data8;
555 int16_t *data16;
556 size_t i;
557
558 switch (bitpix) {
559 case BYTE_IMG:
560 data8 = (BYTE *)from;
561 for (i = 0; i < nbdata; i++)
562 to[i] = (WORD)data8[i];
563 break;
564 case USHORT_IMG: // siril 0.9 native
565 // nothing to do
566 break;
567 case SHORT_IMG:
568 // add 2^15 to the read data to obtain unsigned
569 data16 = (int16_t *)from;
570 for (i = 0; i < nbdata; i++) {
571 int sum = 32768 + (int)data16[i];
572 to[i] = (WORD)sum;
573 }
574 break;
575 case LONGLONG_IMG: // 64-bit integer pixels
576 default:
577 siril_log_message(_("Unknown FITS data format in internal conversion\n"));
578 }
579 }
580
581 /* convert FITS data formats to siril native.
582 * nbdata is the number of pixels, w * h.
583 * from is not freed, to must be allocated and can be the same as from */
convert_data_float(int bitpix,const void * from,float * to,size_t nbdata)584 static void convert_data_float(int bitpix, const void *from, float *to, size_t nbdata) {
585 size_t i;
586 BYTE *data8;
587 WORD *ushort;
588 int16_t *data16;
589 double *pixels_double;
590 long *sdata32; // TO BE TESTED on 32-bit arch, seems to be a cfitsio bug
591 long mini = LONG_MAX;
592 long maxi = -LONG_MAX;
593 unsigned long *data32;
594 float *data32f;
595
596 switch (bitpix) {
597 case BYTE_IMG:
598 data8 = (BYTE *)from;
599 for (i = 0; i < nbdata; i++)
600 to[i] = (float)data8[i] * INV_UCHAR_MAX_SINGLE;
601 break;
602 case USHORT_IMG: // siril 0.9 native
603 ushort = (WORD *)from;
604 for (i = 0; i < nbdata; i++) {
605 to[i] = (float)ushort[i] * INV_USHRT_MAX_SINGLE;
606 }
607 break;
608 case SHORT_IMG:
609 // add 2^15 to the read data to obtain unsigned
610 data16 = (int16_t *)from;
611 for (i = 0; i < nbdata; i++) {
612 int sum = 32768 + (int)data16[i];
613 to[i] = (float)sum * INV_USHRT_MAX_SINGLE;
614 }
615 break;
616 case ULONG_IMG: // 32-bit unsigned integer pixels
617 data32 = (unsigned long *)from;
618 for (i = 0; i < nbdata; i++) {
619 mini = min(data32[i], mini);
620 maxi = max(data32[i], maxi);
621 }
622 for (i = 0; i < nbdata; i++)
623 to[i] = (float)((data32[i] - mini)) / (maxi - mini);
624 break;
625 case LONG_IMG: // 32-bit signed integer pixels
626 sdata32 = (long *)from;
627 for (i = 0; i < nbdata; i++) {
628 mini = min(sdata32[i], mini);
629 maxi = max(sdata32[i], maxi);
630 }
631 for (i = 0; i < nbdata; i++)
632 to[i] = (float)((sdata32[i] - mini)) / (maxi - mini);
633 break;
634 case FLOAT_IMG: // 32-bit floating point pixels, we use it only if float is not in the [0, 1] range
635 data32f = (float *)from;
636 for (i = 0; i < nbdata; i++)
637 to[i] = (data32f[i] / USHRT_MAX_SINGLE);
638 break;
639 case DOUBLE_IMG: // 64-bit floating point pixels
640 pixels_double = (double *)from;
641 for (i = 0; i < nbdata; i++) {
642 to[i] = (float)pixels_double[i];
643 }
644 break;
645 case LONGLONG_IMG: // 64-bit integer pixels
646 default:
647 siril_log_message(_("Unknown FITS data format in internal conversion\n"));
648 }
649 }
650
convert_floats(int bitpix,float * data,size_t nbdata)651 static void convert_floats(int bitpix, float *data, size_t nbdata) {
652 size_t i;
653 switch (bitpix) {
654 case BYTE_IMG:
655 for (i = 0; i < nbdata; i++)
656 data[i] = data[i] * INV_UCHAR_MAX_SINGLE;
657 break;
658 default:
659 case USHORT_IMG: // siril 0.9 native
660 for (i = 0; i < nbdata; i++)
661 data[i] = data[i] * INV_USHRT_MAX_SINGLE;
662 break;
663 case SHORT_IMG:
664 // add 2^15 to the read data to obtain unsigned
665 for (i = 0; i < nbdata; i++) {
666 data[i] = (32768.f + data[i]) * INV_USHRT_MAX_SINGLE;
667 }
668 break;
669 }
670 }
671
get_compression_type(int siril_compression_fits_method)672 static int get_compression_type(int siril_compression_fits_method) {
673 if (siril_compression_fits_method < 0
674 || siril_compression_fits_method > G_N_ELEMENTS(CompressionMethods)) {
675 return -1;
676 } else {
677 return CompressionMethods[siril_compression_fits_method];
678 }
679 }
680
681 /* Move to the first HDU of an opened FIT file
682 * with IMAGE type and with dimension > 0
683 */
siril_fits_move_first_image(fitsfile * fp)684 static int siril_fits_move_first_image(fitsfile* fp) {
685 int status = 0;
686 int naxis = 0;
687
688 /* Move to the first HDU of type image */
689 fits_movabs_hdu(fp, 1, IMAGE_HDU, &status);
690 if (status) {
691 siril_log_message(_("Cannot move to first image in FITS file\n"));
692 report_fits_error(status); /* print error message */
693 return status;
694 }
695
696 /* Find the first HDU of type image with dimension > 0 */
697 do {
698 //Check naxis for current image HDU
699 fits_get_img_dim(fp, &naxis, &status);
700 if (status) {
701 siril_log_message(_("Cannot get dimension of FITS file\n"));
702 report_fits_error(status); /* print error message */
703 return status;
704 }
705 if (naxis > 0) {
706 break;
707 }
708 //Check next image HDU
709 fits_movrel_hdu(fp, 1, IMAGE_HDU, &status);
710 if (status) {
711 siril_log_message(_("Cannot move to next image in FITS file\n"));
712 report_fits_error(status); /* print error message */
713 return status;
714 }
715 } while (!status);
716
717 siril_debug_print("Found image HDU (changed CHDU) with naxis %d (status %d)\n", naxis, status);
718 return status;
719 }
720
721 /* read buffer from an already open FITS file, fit should have all metadata
722 * correct, and convert the buffer to fit->data with the given type, which
723 * currently should be TBYTE or TUSHORT because fit doesn't contain other data.
724 * filename is for error reporting
725 */
read_fits_with_convert(fits * fit,const char * filename,gboolean force_float)726 int read_fits_with_convert(fits* fit, const char* filename, gboolean force_float) {
727 int status = 0, zero = 0, datatype;
728 BYTE *data8;
729 unsigned long *pixels_long;
730 // orig ^ gives the coordinate in each dimension of the first pixel to be read
731 size_t nbpix = fit->naxes[0] * fit->naxes[1];
732 size_t nbdata = nbpix * fit->naxes[2];
733 // with force_float, image is read as float data, type is stored as DATA_FLOAT
734 int fake_bitpix = force_float ? FLOAT_IMG : fit->bitpix;
735
736 switch (fake_bitpix) {
737 case BYTE_IMG:
738 case SHORT_IMG:
739 case USHORT_IMG:
740 /* we store these types as unsigned short */
741 if ((fit->data = malloc(nbdata * sizeof(WORD))) == NULL) {
742 PRINT_ALLOC_ERR;
743 return -1;
744 }
745 fit->pdata[RLAYER] = fit->data;
746 if (fit->naxis == 3) {
747 fit->pdata[GLAYER] = fit->data + nbpix;
748 fit->pdata[BLAYER] = fit->data + nbpix * 2;
749 } else {
750 fit->pdata[GLAYER] = fit->data;
751 fit->pdata[BLAYER] = fit->data;
752 }
753 fit->type = DATA_USHORT;
754 break;
755
756 case ULONG_IMG: // 32-bit unsigned integer pixels
757 case LONG_IMG: // 32-bit signed integer pixels
758 case DOUBLE_IMG: // 64-bit floating point pixels
759 case FLOAT_IMG: // 32-bit floating point pixels
760 /* we store these types as float */
761 if ((fit->fdata = malloc(nbdata * sizeof(float))) == NULL) {
762 PRINT_ALLOC_ERR;
763 return -1;
764 }
765 fit->fpdata[RLAYER] = fit->fdata;
766 if (fit->naxis == 3) {
767 fit->fpdata[GLAYER] = fit->fdata + nbpix;
768 fit->fpdata[BLAYER] = fit->fdata + nbpix * 2;
769 } else {
770 fit->fpdata[GLAYER] = fit->fdata;
771 fit->fpdata[BLAYER] = fit->fdata;
772 }
773 fit->type = DATA_FLOAT;
774 break;
775
776 case LONGLONG_IMG: // 64-bit integer pixels
777 default:
778 siril_log_message(_("FITS image format %d is not supported by Siril.\n"), fit->bitpix);
779 return -1;
780 }
781
782 status = 0;
783 switch (fake_bitpix) {
784 case BYTE_IMG:
785 data8 = malloc(nbdata * sizeof(BYTE));
786 datatype = fit->bitpix == BYTE_IMG ? TBYTE : TSBYTE;
787 fits_read_img(fit->fptr, datatype, 1, nbdata, &zero, data8, &zero, &status);
788 if (status) break;
789 convert_data_ushort(fit->bitpix, data8, fit->data, nbdata, FALSE);
790 free(data8);
791 break;
792 case SHORT_IMG:
793 fits_read_img(fit->fptr, TSHORT, 1, nbdata, &zero, fit->data, &zero, &status);
794 if (status) break;
795 convert_data_ushort(fit->bitpix, fit->data, fit->data, nbdata, FALSE);
796 fit->bitpix = USHORT_IMG;
797 break;
798 case USHORT_IMG:
799 // siril 0.9 native, no conversion required
800 fits_read_img(fit->fptr, TUSHORT, 1, nbdata, &zero, fit->data, &zero, &status);
801 if (status == NUM_OVERFLOW) {
802 // in case there are errors, we try short data
803 status = 0;
804 fits_read_img(fit->fptr, TSHORT, 1, nbdata, &zero, fit->data,
805 &zero, &status);
806 if (status)
807 break;
808 convert_data_ushort(SHORT_IMG, fit->data, fit->data, nbdata, FALSE);
809 if (fit->lo)
810 fit->lo += 32768;
811 if (fit->hi)
812 fit->hi += 32768;
813 fit->bitpix = USHORT_IMG;
814 }
815 break;
816
817 case ULONG_IMG: // 32-bit unsigned integer pixels
818 case LONG_IMG: // 32-bit signed integer pixels
819 pixels_long = malloc(nbdata * sizeof(unsigned long));
820 datatype = fit->bitpix == LONG_IMG ? TLONG : TULONG;
821 fits_read_img(fit->fptr, datatype, 1, nbdata, &zero, pixels_long, &zero, &status);
822 if (status) break;
823 convert_data_float(fit->bitpix, pixels_long, fit->fdata, nbdata);
824 free(pixels_long);
825 fit->bitpix = FLOAT_IMG;
826 break;
827 case FLOAT_IMG: // 32-bit floating point pixels
828 // siril 1.0 native, no conversion required
829 case DOUBLE_IMG: // 64-bit floating point pixels
830 // let cfitsio do the conversion
831 /* we assume we are in the range [0, 1]. But, for some images
832 * some values can be negative
833 */
834 fits_read_img(fit->fptr, TFLOAT, 1, nbdata, &zero, fit->fdata, &zero,
835 &status);
836 if ((fit->bitpix == USHORT_IMG || fit->bitpix == SHORT_IMG
837 || fit->bitpix == BYTE_IMG) || fit->data_max > 2.0) { // needed for some FLOAT_IMG
838 convert_floats(fit->bitpix, fit->fdata, nbdata);
839 }
840 fit->bitpix = FLOAT_IMG;
841 fit->orig_bitpix = FLOAT_IMG; // force this, to avoid problems saving the FITS if needed
842 break;
843 }
844
845 if (status) {
846 siril_log_message(_("Fitsio error reading data, file: %s.\n"), filename);
847 report_fits_error(status);
848 return -1;
849 }
850
851 return 0;
852 }
853
854 /* This function reads partial data on one layer from the opened FITS and
855 * convert it to siril's format (USHORT) */
internal_read_partial_fits(fitsfile * fptr,unsigned int ry,int bitpix,void * dest,int layer,const rectangle * area)856 int internal_read_partial_fits(fitsfile *fptr, unsigned int ry,
857 int bitpix, void *dest, int layer, const rectangle *area) {
858 int datatype;
859 BYTE *data8;
860 long *pixels_long;
861 long fpixel[3], lpixel[3], inc[3] = { 1L, 1L, 1L };
862 int zero = 0, status = 0;
863
864 /* fpixel is first pixel, lpixel is last pixel, starts with value 1 */
865 fpixel[0] = area->x + 1; // in siril, it starts with 0
866 fpixel[1] = ry - area->y - area->h + 1;
867 fpixel[2] = layer + 1;
868 lpixel[0] = area->x + area->w; // with w and h at least 1, we're ok
869 lpixel[1] = ry - area->y;
870 lpixel[2] = layer + 1;
871
872 size_t nbdata = area->w * area->h;
873
874 switch (bitpix) {
875 case BYTE_IMG:
876 data8 = malloc(nbdata * sizeof(BYTE));
877 datatype = bitpix == BYTE_IMG ? TBYTE : TSBYTE;
878 fits_read_subset(fptr, datatype, fpixel, lpixel, inc, &zero, data8,
879 &zero, &status);
880 if (status) break;
881 convert_data_ushort(bitpix, data8, dest, nbdata, FALSE);
882 free(data8);
883 break;
884 case SHORT_IMG:
885 fits_read_subset(fptr, TSHORT, fpixel, lpixel, inc, &zero, dest,
886 &zero, &status);
887 convert_data_ushort(bitpix, dest, dest, nbdata, FALSE);
888 break;
889 case USHORT_IMG:
890 fits_read_subset(fptr, TUSHORT, fpixel, lpixel, inc, &zero, dest,
891 &zero, &status);
892 break;
893
894 /* types below are read as a 32-bit float image: this breaks compatibility
895 * with functions working only with 16-bit int data */
896 case ULONG_IMG: // 32-bit unsigned integer pixels
897 case LONG_IMG: // 32-bit signed integer pixels
898 pixels_long = malloc(nbdata * sizeof(long));
899 status = 0;
900 datatype = bitpix == LONG_IMG ? TLONG : TULONG;
901 fits_read_subset(fptr, datatype, fpixel, lpixel, inc, &zero,
902 pixels_long, &zero, &status);
903 if (status) break;
904 convert_data_float(bitpix, pixels_long, dest, nbdata);
905 free(pixels_long);
906 break;
907 case DOUBLE_IMG: // 64-bit floating point pixels
908 case FLOAT_IMG: // 32-bit floating point pixels
909 fits_read_subset(fptr, TFLOAT, fpixel, lpixel, inc, &zero, dest,
910 &zero, &status);
911 break;
912 case LONGLONG_IMG: // 64-bit integer pixels
913 default:
914 siril_log_message(_("FITS image format %d is not supported by Siril.\n"), bitpix);
915 return -1;
916 }
917 return status;
918 }
919
siril_fits_create_diskfile(fitsfile ** fptr,const char * filename,int * status)920 int siril_fits_create_diskfile(fitsfile **fptr, const char *filename, int *status) {
921 gchar *localefilename = get_locale_filename(filename);
922 fits_create_diskfile(fptr, localefilename, status);
923 g_free(localefilename);
924 return *status;
925 }
926
save_wcs_keywords(fits * fit)927 static void save_wcs_keywords(fits *fit) {
928 int status = 0;
929
930 if (fit->wcsdata.equinox > 0.0) {
931 fits_update_key(fit->fptr, TSTRING, "CTYPE1", "RA---TAN", "Coordinate type for the first axis", &status);
932 status = 0;
933 fits_update_key(fit->fptr, TSTRING, "CTYPE2", "DEC--TAN", "Coordinate type for the second axis", &status);
934 status = 0;
935 fits_update_key(fit->fptr, TDOUBLE, "EQUINOX", &(fit->wcsdata.equinox), "Equatorial equinox", &status);
936 }
937 status = 0;
938
939 /* Needed for Aladin compatibility */
940 if (fit->naxes[2] == 3 && com.pref.rgb_aladin) {
941 status = 0;
942 fits_update_key(fit->fptr, TSTRING, "CTYPE3", "RGB", "RGB image", &status);
943 }
944
945 if (fit->wcsdata.objctra[0] != '\0') {
946 status = 0;
947 fits_update_key(fit->fptr, TSTRING, "OBJCTRA", &(fit->wcsdata.objctra), "Image center R.A. (hms)", &status);
948 status = 0;
949 fits_update_key(fit->fptr, TSTRING, "OBJCTDEC", &(fit->wcsdata.objctdec), "Image center declination (dms)", &status);
950 }
951 if (fit->wcsdata.crpix[0] != '\0') {
952 status = 0;
953 fits_update_key(fit->fptr, TDOUBLE, "CRPIX1", &(fit->wcsdata.crpix[0]), "Axis1 reference pixel", &status);
954 status = 0;
955 fits_update_key(fit->fptr, TDOUBLE, "CRPIX2", &(fit->wcsdata.crpix[1]), "Axis2 reference pixel", &status);
956 }
957 if (fit->wcsdata.crval[0] != '\0') {
958 status = 0;
959 fits_update_key(fit->fptr, TDOUBLE, "CRVAL1", &(fit->wcsdata.crval[0]), "Axis1 reference value", &status);
960 status = 0;
961 fits_update_key(fit->fptr, TDOUBLE, "CRVAL2", &(fit->wcsdata.crval[1]), "Axis2 reference value", &status);
962 }
963 if ((fit->wcsdata.cd[0][0] != 0.0) && (fit->wcsdata.cd[0][1] != 0.0) && (fit->wcsdata.cd[1][0] != 0.0) && (fit->wcsdata.cd[1][1] != 0.0)) {
964 status = 0;
965 fits_update_key(fit->fptr, TDOUBLE, "CD1_1", &(fit->wcsdata.cd[0][0]), "Scale matrix (1, 1)", &status);
966 status = 0;
967 fits_update_key(fit->fptr, TDOUBLE, "CD1_2", &(fit->wcsdata.cd[0][1]), "Scale matrix (1, 2)", &status);
968 status = 0;
969 fits_update_key(fit->fptr, TDOUBLE, "CD2_1", &(fit->wcsdata.cd[1][0]), "Scale matrix (2, 1)", &status);
970 status = 0;
971 fits_update_key(fit->fptr, TDOUBLE, "CD2_2", &(fit->wcsdata.cd[1][1]), "Scale matrix (2, 2)", &status);
972 status = 0;
973 fits_update_key(fit->fptr, TINT, "IMAGEW", &(fit->rx), "Image width, in pixels.", &status);
974 status = 0;
975 fits_update_key(fit->fptr, TINT, "IMAGEH", &(fit->ry), "Image height, in pixels.", &status);
976 }
977 }
978
save_fits_header(fits * fit)979 void save_fits_header(fits *fit) {
980 int i, status = 0;
981 double zero, scale;
982 char comment[FLEN_COMMENT];
983
984 if (fit->hi) { /* may not be initialized */
985 if (fit->type == DATA_USHORT) {
986 fits_update_key(fit->fptr, TUSHORT, "MIPS-HI", &(fit->hi),
987 "Upper visualization cutoff", &status);
988 fits_update_key(fit->fptr, TUSHORT, "MIPS-LO", &(fit->lo),
989 "Lower visualization cutoff", &status);
990 }
991 else if (fit->type == DATA_FLOAT) {
992 float fhi = ushort_to_float_range(fit->hi);
993 fits_update_key(fit->fptr, TFLOAT, "MIPS-FHI", &(fhi),
994 "Upper visualization cutoff", &status);
995 float flo = ushort_to_float_range(fit->lo);
996 fits_update_key(fit->fptr, TFLOAT, "MIPS-FLO", &(flo),
997 "Lower visualization cutoff", &status);
998 }
999 }
1000 // physical_value = BZERO + BSCALE * array_value
1001 // https://heasarc.gsfc.nasa.gov/docs/software/fitsio/c/c_user/node26.html
1002 switch (fit->bitpix) {
1003 case BYTE_IMG:
1004 case SHORT_IMG:
1005 zero = 0.0;
1006 scale = 1.0;
1007 break;
1008 case FLOAT_IMG:
1009 zero = 0.0;
1010 scale = 1.0;
1011 break;
1012 default:
1013 case USHORT_IMG:
1014 zero = 32768.0;
1015 scale = 1.0;
1016 break;
1017 }
1018 status = 0;
1019 fits_update_key(fit->fptr, TDOUBLE, "BZERO", &zero,
1020 "offset data range to that of unsigned short", &status);
1021
1022 status = 0;
1023 fits_update_key(fit->fptr, TDOUBLE, "BSCALE", &scale, "default scaling factor",
1024 &status);
1025
1026 status = 0;
1027 if (!g_strcmp0(fit->row_order, "BOTTOM-UP") || !g_strcmp0(fit->row_order, "TOP-DOWN")) {
1028 fits_update_key(fit->fptr, TSTRING, "ROWORDER", &fit->row_order,
1029 "Order of the rows in image array", &status);
1030 }
1031
1032 /*******************************************************************
1033 * ************* CAMERA AND INSTRUMENT KEYWORDS ********************
1034 * ******************** AND DATES **********************************
1035 * ****************************************************************/
1036
1037 status = 0;
1038 if (fit->instrume[0] != '\0')
1039 fits_update_key(fit->fptr, TSTRING, "INSTRUME", &(fit->instrume),
1040 "instrument name", &status);
1041 status = 0;
1042 if (fit->telescop[0] != '\0')
1043 fits_update_key(fit->fptr, TSTRING, "TELESCOP", &(fit->telescop),
1044 "telescope used to acquire this image", &status);
1045 status = 0;
1046 if (fit->observer[0] != '\0')
1047 fits_update_key(fit->fptr, TSTRING, "OBSERVER", &(fit->observer),
1048 "observer name", &status);
1049 status = 0;
1050 int itmp;
1051 char fit_date[40];
1052 fits_get_system_time(fit_date, &itmp, &status);
1053 fits_update_key(fit->fptr, TSTRING, "DATE", fit_date,
1054 "UTC date that FITS file was created", &status);
1055
1056 status = 0;
1057 if (fit->date_obs) {
1058 gchar *formatted_date = date_time_to_FITS_date(fit->date_obs);
1059 fits_update_key(fit->fptr, TSTRING, "DATE-OBS", formatted_date,
1060 "YYYY-MM-DDThh:mm:ss observation start, UT", &status);
1061
1062 g_free(formatted_date);
1063 }
1064
1065 status = 0;
1066 if (fit->exposure > 0.)
1067 fits_update_key(fit->fptr, TDOUBLE, "EXPTIME", &(fit->exposure),
1068 "Exposure time [s]", &status);
1069
1070 status = 0;
1071 if (fit->expstart > 0.)
1072 fits_update_key(fit->fptr, TDOUBLE, "EXPSTART", &(fit->expstart),
1073 "Exposure start time (standard Julian date)", &status);
1074
1075 status = 0;
1076 if (fit->expend > 0.)
1077 fits_update_key(fit->fptr, TDOUBLE, "EXPEND", &(fit->expend),
1078 "Exposure end time (standard Julian date)", &status);
1079
1080 /* all keywords below are non-standard */
1081 status = 0;
1082 if (fit->pixel_size_x > 0.)
1083 fits_update_key(fit->fptr, TFLOAT, "XPIXSZ", &(fit->pixel_size_x),
1084 "X pixel size microns", &status);
1085 if (fit->pixel_size_y > 0.)
1086 fits_update_key(fit->fptr, TFLOAT, "YPIXSZ", &(fit->pixel_size_y),
1087 "Y pixel size microns", &status);
1088
1089 status = 0;
1090 if (fit->binning_x)
1091 fits_update_key(fit->fptr, TUINT, "XBINNING", &(fit->binning_x),
1092 "Camera binning mode", &status);
1093 if (fit->binning_y)
1094 fits_update_key(fit->fptr, TUINT, "YBINNING", &(fit->binning_y),
1095 "Camera binning mode", &status);
1096
1097 status = 0;
1098 if (fit->focal_length > 0.)
1099 fits_update_key(fit->fptr, TDOUBLE, "FOCALLEN", &(fit->focal_length),
1100 "Camera focal length", &status);
1101
1102 status = 0;
1103 if (fit->ccd_temp)
1104 fits_update_key(fit->fptr, TDOUBLE, "CCD-TEMP", &(fit->ccd_temp),
1105 "CCD temp in C", &status);
1106
1107 status = 0;
1108 if (fit->aperture > 0.)
1109 fits_update_key(fit->fptr, TDOUBLE, "APERTURE", &(fit->aperture),
1110 "Aperture of the instrument", &status);
1111
1112 status = 0;
1113 if (fit->iso_speed > 0.)
1114 fits_update_key(fit->fptr, TDOUBLE, "ISOSPEED", &(fit->iso_speed),
1115 "ISO camera setting", &status);
1116
1117 status = 0;
1118 if (fit->bayer_pattern[0] != '\0') {
1119 fits_update_key(fit->fptr, TSTRING, "BAYERPAT", &(fit->bayer_pattern),
1120 "Bayer color pattern", &status);
1121
1122 status = 0;
1123 fits_update_key(fit->fptr, TINT, "XBAYROFF", &(fit->bayer_xoffset),
1124 "X offset of Bayer array", &status);
1125
1126 status = 0;
1127 fits_update_key(fit->fptr, TINT, "YBAYROFF", &(fit->bayer_yoffset),
1128 "Y offset of Bayer array", &status);
1129
1130 }
1131
1132 status = 0;
1133 if (fit->key_gain > 0.)
1134 fits_update_key(fit->fptr, TUSHORT, "GAIN", &(fit->key_gain),
1135 "Camera gain", &status);
1136
1137 status = 0;
1138 if (fit->key_offset > 0.)
1139 fits_update_key(fit->fptr, TUSHORT, "OFFSET", &(fit->key_offset),
1140 "Camera offset", &status);
1141
1142 status = 0;
1143 if (fit->cvf > 0.)
1144 fits_update_key(fit->fptr, TDOUBLE, "CVF", &(fit->cvf),
1145 "Conversion factor (e-/adu)", &status);
1146
1147 /*******************************************************************
1148 * ******************* PLATE SOLVING KEYWORDS **********************
1149 * ****************************************************************/
1150
1151 save_wcs_keywords(fit);
1152
1153 /*******************************************************************
1154 * ******************** PROGRAMM KEYWORDS **************************
1155 * ****************************************************************/
1156
1157 status = 0;
1158 char programm[32];
1159 sprintf(programm, "%s v%s", PACKAGE, VERSION);
1160 programm[0] = toupper(programm[0]); // convert siril to Siril
1161 fits_update_key(fit->fptr, TSTRING, "PROGRAM", programm,
1162 "Software that created this HDU", &status);
1163
1164 /*******************************************************************
1165 * ********************* HISTORY KEYWORDS **************************
1166 * ****************************************************************/
1167
1168 status = 0;
1169 if (fit->history) {
1170 GSList *list;
1171 for (list = fit->history; list; list = list->next) {
1172 fits_write_history(fit->fptr, (char *)list->data, &status);
1173 }
1174 }
1175
1176 status = 0;
1177 if (com.history) {
1178 for (i = 0; i < com.hist_display; i++) {
1179 if (com.history[i].history[0] != '\0')
1180 fits_write_history(fit->fptr, com.history[i].history, &status);
1181 }
1182 }
1183
1184 /*******************************************************************
1185 * ************************* DFT KEYWORDS **************************
1186 * ****************************************************************/
1187
1188 status = 0;
1189 if (fit->dft.type[0] != '\0') {
1190 if (fit->dft.type[0] == 'S')
1191 strcpy(comment, "Module of a Discrete Fourier Transform");
1192 else if (fit->dft.type[0] == 'P')
1193 strcpy(comment, "Phase of a Discrete Fourier Transform");
1194 else
1195 status = 1; // should not happen
1196 fits_update_key(fit->fptr, TSTRING, "DFTTYPE", &(fit->dft.type),
1197 comment, &status);
1198 }
1199
1200 status = 0;
1201 if (fit->dft.ord[0] != '\0') {
1202 if (fit->dft.ord[0] == 'C')
1203 strcpy(comment, "Low spatial freq. are located at image center");
1204 else if (fit->dft.ord[0] == 'R')
1205 strcpy(comment, "High spatial freq. are located at image center");
1206 else
1207 status = 1; // should not happen
1208 fits_update_key(fit->fptr, TSTRING, "DFTORD", &(fit->dft.ord), comment,
1209 &status);
1210 }
1211
1212 for (i = 0; i < fit->naxes[2]; i++) {
1213 if (fit->dft.norm[i] > 0.) {
1214 const char *str1 = "DFTNORM";
1215 const char *str2 = "Normalisation value for channel #";
1216 char key_str[FLEN_KEYWORD], comment_str[FLEN_VALUE];
1217 g_sprintf(key_str, "%s%d", str1, i);
1218 g_sprintf(comment_str, "%s%d", str2, i);
1219 status = 0;
1220 fits_update_key(fit->fptr, TDOUBLE, key_str, &(fit->dft.norm[i]),
1221 comment_str, &status);
1222 }
1223 }
1224 }
1225
1226 /********************** public functions ************************************/
1227
get_date_data_from_fitsfile(fitsfile * fptr,GDateTime ** dt,double * exposure)1228 void get_date_data_from_fitsfile(fitsfile *fptr, GDateTime **dt, double *exposure) {
1229 char date_obs[FLEN_VALUE];
1230
1231 __tryToFindKeywords(fptr, TDOUBLE, Exposure, exposure);
1232 int status = 0;
1233 fits_read_key(fptr, TSTRING, "DATE-OBS", &date_obs, NULL, &status);
1234 if (!status) {
1235 *dt = FITS_date_to_date_time(date_obs);
1236 }
1237 }
1238
import_metadata_from_fitsfile(fitsfile * fptr,fits * to)1239 int import_metadata_from_fitsfile(fitsfile *fptr, fits *to) {
1240 fits from = { 0 };
1241 from.fptr = fptr;
1242 read_fits_header(&from);
1243 copy_fits_metadata(&from, to);
1244 return 0;
1245 }
1246
1247 /* from bitpix, depending on BZERO, bitpix and orig_bitpix are set.
1248 *
1249 * since USHORT_IMG is a cfitsio trick and doesn't really exist in the
1250 * file, when we read it, the bitpix is given as SHORT_IMG. If BZERO is
1251 * 2^15, it means that it is USHORT data and we force the bitpix to it
1252 * in order to read it properly. Same thing for LONG and ULONG.
1253 * https://heasarc.gsfc.nasa.gov/docs/software/fitsio/c/c_user/node23.html
1254 */
manage_bitpix(fitsfile * fptr,int * bitpix,int * orig_bitpix)1255 void manage_bitpix(fitsfile *fptr, int *bitpix, int *orig_bitpix) {
1256 double offset;
1257 int status = 0;
1258 fits_read_key(fptr, TDOUBLE, "BZERO", &offset, NULL, &status);
1259 if (!status) {
1260 if (*bitpix == SHORT_IMG && offset != 0.0) {
1261 *bitpix = USHORT_IMG;
1262 }
1263 else if (*bitpix == LONG_IMG && offset != 0.0)
1264 *bitpix = ULONG_IMG;
1265 } else {
1266 /* but some software just put unsigned 16-bit data in the file
1267 * and don't set the BZERO keyword... */
1268 if (status == KEY_NO_EXIST && *bitpix == SHORT_IMG)
1269 *bitpix = USHORT_IMG;
1270 }
1271 // and we store the original bitpix to reuse it during later partial
1272 // reads when we have no access to the header or a fits * struct.
1273 *orig_bitpix = *bitpix;
1274 }
1275
1276 // return 0 on success, fills realname if not NULL with the opened file's name
readfits(const char * filename,fits * fit,char * realname,gboolean force_float)1277 int readfits(const char *filename, fits *fit, char *realname, gboolean force_float) {
1278 int status, retval = 1;
1279 char *name = NULL;
1280 gchar *basename;
1281 image_type imagetype;
1282
1283 if (stat_file(filename, &imagetype, &name)) {
1284 siril_log_message(_("%s.[any_allowed_extension] not found.\n"),
1285 filename);
1286 free(name);
1287 return 1;
1288 }
1289 if (imagetype != TYPEFITS) {
1290 siril_log_message(
1291 _("The file %s is not a FITS file or doesn't exists with FITS extensions.\n"),
1292 filename);
1293 free(name);
1294 return 1;
1295 }
1296
1297 if (realname)
1298 strcpy(realname, name);
1299
1300 status = 0;
1301 siril_fits_open_diskfile(&(fit->fptr), name, READONLY, &status);
1302 if (status) {
1303 report_fits_error(status);
1304 free(name);
1305 return status;
1306 }
1307 free(name);
1308
1309 if (siril_fits_move_first_image(fit->fptr)) {
1310 siril_log_message(_("Selecting the primary header failed, is the FITS file '%s' malformed?\n"), filename);
1311 goto close_readfits;
1312 }
1313
1314 status = read_fits_metadata(fit);
1315 if (status)
1316 goto close_readfits;
1317
1318 retval = read_fits_with_convert(fit, filename, force_float);
1319 fit->top_down = FALSE;
1320
1321 if (!retval) {
1322 // copy the entire header in memory
1323 if (fit->header)
1324 free(fit->header);
1325 fit->header = copy_header(fit);
1326
1327 basename = g_path_get_basename(filename);
1328 siril_log_message(_("Reading FITS: file %s, %ld layer(s), %ux%u pixels\n"),
1329 basename, fit->naxes[2], fit->rx, fit->ry) ;
1330 g_free(basename);
1331 }
1332
1333 close_readfits:
1334 status = 0;
1335 fits_close_file(fit->fptr, &status);
1336 return retval;
1337 }
1338
siril_fits_open_diskfile(fitsfile ** fptr,const char * filename,int iomode,int * status)1339 int siril_fits_open_diskfile(fitsfile **fptr, const char *filename, int iomode, int *status) {
1340 gchar *localefilename = get_locale_filename(filename);
1341 fits_open_diskfile(fptr, localefilename, iomode, status);
1342 if (!(*status)) {
1343 siril_fits_move_first_image(*fptr);
1344 }
1345 g_free(localefilename);
1346 return *status;
1347 }
1348
1349 // reset a fit data structure, deallocates everything in it and zero the data
clearfits(fits * fit)1350 void clearfits(fits *fit) {
1351 if (fit == NULL)
1352 return;
1353 if (fit->data)
1354 free(fit->data);
1355 if (fit->fdata)
1356 free(fit->fdata);
1357 if (fit->header)
1358 free(fit->header);
1359 if (fit->history)
1360 g_slist_free_full(fit->history, free);
1361 if (fit->date_obs)
1362 g_date_time_unref(fit->date_obs);
1363 if (fit->date)
1364 g_date_time_unref(fit->date);
1365 if (fit->stats) {
1366 for (int i = 0; i < fit->naxes[2]; i++)
1367 free_stats(fit->stats[i]);
1368 free(fit->stats);
1369 }
1370 free_wcs(fit);
1371
1372 memset(fit, 0, sizeof(fits));
1373 }
1374
1375 /* Read a rectangular section of a FITS image in Siril's format, pointed by its
1376 * exact filename. Only layer layer is read.
1377 * Returned fit->data is upside-down. */
readfits_partial(const char * filename,int layer,fits * fit,const rectangle * area,gboolean do_photometry)1378 int readfits_partial(const char *filename, int layer, fits *fit,
1379 const rectangle *area, gboolean do_photometry) {
1380 int status;
1381 size_t nbdata;
1382 double data_max = 0.0;
1383
1384 status = 0;
1385 if (siril_fits_open_diskfile(&(fit->fptr), filename, READONLY, &status)) {
1386 report_fits_error(status);
1387 return status;
1388 }
1389
1390 if (siril_fits_move_first_image(fit->fptr)) {
1391 siril_log_message(_("Selecting the primary header failed, is the FITS file '%s' malformed?\n"), filename);
1392 return -1;
1393 }
1394
1395 status = 0;
1396 fits_get_img_param(fit->fptr, 3, &(fit->bitpix), &(fit->naxis), fit->naxes,
1397 &status);
1398 if (status) {
1399 report_fits_error(status);
1400 status = 0;
1401 fits_close_file(fit->fptr, &status);
1402 return status;
1403 }
1404
1405 manage_bitpix(fit->fptr, &(fit->bitpix), &(fit->orig_bitpix));
1406
1407 if (do_photometry)
1408 fit_get_photometry_data(fit);
1409
1410 if (fit->naxis == 2 && fit->naxes[2] == 0)
1411 fit->naxes[2] = 1; // see readfits for the explanation
1412 if (layer > fit->naxes[2] - 1) {
1413 siril_log_message(_("FITS read partial: there is no layer %d in the image %s\n"),
1414 layer + 1, filename);
1415 status = 0;
1416 fits_close_file(fit->fptr, &status);
1417 return -1;
1418 }
1419
1420 if (fit->naxis == 3 && fit->naxes[2] != 3) {
1421 siril_log_message(_("Unsupported FITS image format (%ld axes).\n"),
1422 fit->naxes[2]);
1423 status = 0;
1424 fits_close_file(fit->fptr, &status);
1425 return -1;
1426 }
1427
1428 nbdata = area->w * area->h;
1429 fit->type = get_data_type(fit->bitpix);
1430 if (fit->type == DATA_UNSUPPORTED) {
1431 siril_log_message(_("Unknown FITS data format in internal conversion\n"));
1432 fits_close_file(fit->fptr, &status);
1433 return -1;
1434 }
1435
1436 if (fit->type == DATA_USHORT) {
1437 /* realloc fit->data to the image size */
1438 WORD *olddata = fit->data;
1439 if ((fit->data = realloc(fit->data, nbdata * sizeof(WORD))) == NULL) {
1440 PRINT_ALLOC_ERR;
1441 status = 0;
1442 fits_close_file(fit->fptr, &status);
1443 if (olddata)
1444 free(olddata);
1445 return -1;
1446 }
1447 fit->pdata[RLAYER] = fit->data;
1448 fit->pdata[GLAYER] = fit->data;
1449 fit->pdata[BLAYER] = fit->data;
1450
1451 status = internal_read_partial_fits(fit->fptr, fit->naxes[1],
1452 fit->bitpix, fit->data, layer, area);
1453 } else {
1454 if (fit->bitpix == FLOAT_IMG) {
1455 status = 0;
1456 fits_read_key(fit->fptr, TDOUBLE, "DATAMAX", &data_max, NULL, &status);
1457 if (status == KEY_NO_EXIST) {
1458 float mini, maxi;
1459 fit_stats(fit, &mini, &maxi);
1460 fit->data_max = (double) maxi;
1461 }
1462 }
1463
1464 /* realloc fit->fdata to the image size */
1465 float *olddata = fit->fdata;
1466 if ((fit->fdata = realloc(fit->fdata, nbdata * sizeof(float))) == NULL) {
1467 PRINT_ALLOC_ERR;
1468 status = 0;
1469 fits_close_file(fit->fptr, &status);
1470 if (olddata)
1471 free(olddata);
1472 return -1;
1473 }
1474 fit->fpdata[RLAYER] = fit->fdata;
1475 fit->fpdata[GLAYER] = fit->fdata;
1476 fit->fpdata[BLAYER] = fit->fdata;
1477
1478 status = internal_read_partial_fits(fit->fptr, fit->naxes[1],
1479 fit->bitpix, fit->fdata, layer, area);
1480 }
1481
1482 if (status) {
1483 report_fits_error(status);
1484 status = 0;
1485 fits_close_file(fit->fptr, &status);
1486 return -1;
1487 }
1488
1489 fit->naxes[0] = area->w;
1490 fit->naxes[1] = area->h;
1491 fit->rx = fit->naxes[0];
1492 fit->ry = fit->naxes[1];
1493 fit->naxes[2] = 1;
1494 fit->naxis = 2;
1495
1496 status = 0;
1497 fits_close_file(fit->fptr, &status);
1498 fprintf(stdout, _("Loaded partial FITS file %s\n"), filename);
1499 return 0;
1500 }
1501
read_fits_metadata(fits * fit)1502 int read_fits_metadata(fits *fit) {
1503 int status = 0;
1504 fit->naxes[2] = 1;
1505 fits_get_img_param(fit->fptr, 3, &(fit->bitpix), &(fit->naxis), fit->naxes, &status);
1506 if (status) {
1507 siril_log_message(_("FITSIO error getting image parameters.\n"));
1508 report_fits_error(status);
1509 status = 0;
1510 fits_close_file(fit->fptr, &status);
1511 return 1;
1512 }
1513
1514 manage_bitpix(fit->fptr, &(fit->bitpix), &(fit->orig_bitpix));
1515
1516 fit->rx = fit->naxes[0];
1517 fit->ry = fit->naxes[1];
1518
1519 if (fit->naxis == 3 && fit->naxes[2] != 3) {
1520 siril_log_color_message(_("The FITS image contains more than 3 channels (%ld). Opening only the three first.\n"), "salmon", fit->naxes[2]);
1521 if (fit->naxis == 3) fit->naxes[2] = 3;
1522 }
1523
1524 if (fit->naxis == 2 && fit->naxes[2] == 0) {
1525 fit->naxes[2] = 1;
1526 /* naxes[2] is set to 1 because:
1527 * - it doesn't matter, since naxis is 2, it's not used
1528 * - it's very convenient to use it in multiplications as the number of layers
1529 */
1530 }
1531 if (fit->bitpix == LONGLONG_IMG) {
1532 siril_log_message(
1533 _("FITS images with 64 bits signed integer per pixel.channel are not supported.\n"));
1534 status = 0;
1535 fits_close_file(fit->fptr, &status);
1536 return -1;
1537 }
1538
1539 read_fits_header(fit); // stores useful header data in fit
1540 return 0;
1541 }
1542
read_fits_metadata_from_path(const char * filename,fits * fit)1543 int read_fits_metadata_from_path(const char *filename, fits *fit) {
1544 int status = 0;
1545 siril_fits_open_diskfile(&(fit->fptr), filename, READONLY, &status);
1546 if (status) {
1547 report_fits_error(status);
1548 return status;
1549 }
1550
1551 if (siril_fits_move_first_image(fit->fptr)) {
1552 siril_log_message(_("Selecting the primary header failed, is the FITS file '%s' malformed?\n"), filename);
1553 return -1;
1554 }
1555
1556 read_fits_metadata(fit);
1557
1558 status = 0;
1559 fits_close_file(fit->fptr, &status);
1560 return status;
1561 }
1562
flip_buffer(int bitpix,void * buffer,const rectangle * area)1563 void flip_buffer(int bitpix, void *buffer, const rectangle *area) {
1564 /* reverse the read data, because it's stored upside-down */
1565 if (get_data_type(bitpix) == DATA_FLOAT) {
1566 int line_size = area->w * sizeof(float);
1567 void *swap = malloc(line_size);
1568 float *buf = (float *)buffer;
1569 int i;
1570 for (i = 0; i < area->h/2 ; i++) {
1571 memcpy(swap, buf + i*area->w, line_size);
1572 memcpy(buf + i*area->w, buf + (area->h - i - 1)*area->w, line_size);
1573 memcpy(buf + (area->h - i - 1)*area->w, swap, line_size);
1574 }
1575 free(swap);
1576 } else {
1577 int line_size = area->w * sizeof(WORD);
1578 void *swap = malloc(line_size);
1579 WORD *buf = (WORD *)buffer;
1580 int i;
1581 for (i = 0; i < area->h/2 ; i++) {
1582 memcpy(swap, buf + i*area->w, line_size);
1583 memcpy(buf + i*area->w, buf + (area->h - i - 1)*area->w, line_size);
1584 memcpy(buf + (area->h - i - 1)*area->w, swap, line_size);
1585 }
1586 free(swap);
1587 }
1588 }
1589
1590
1591 /* read subset of an opened fits file.
1592 * The rectangle's coordinates x,y start at 0,0 for first pixel in the image.
1593 * layer and index also start at 0.
1594 * buffer has to be allocated with enough space to store the area.
1595 */
read_opened_fits_partial(sequence * seq,int layer,int index,void * buffer,const rectangle * area)1596 int read_opened_fits_partial(sequence *seq, int layer, int index, void *buffer,
1597 const rectangle *area) {
1598 int status;
1599
1600 if (!seq || !seq->fptr || !seq->fptr[index]) {
1601 printf("data initialization error in read fits partial\n");
1602 return 1;
1603 }
1604 if (area->x < 0 || area->y < 0 || area->x >= seq->rx || area->y >= seq->ry
1605 || area->w <= 0 || area->h <= 0 || area->x + area->w > seq->rx
1606 || area->y + area->h > seq->ry) {
1607 fprintf(stderr, "partial read from FITS file has been requested outside image bounds or with invalid size\n");
1608 return 1;
1609 }
1610
1611 #ifdef _OPENMP
1612 g_assert(seq->fd_lock);
1613 omp_set_lock(&seq->fd_lock[index]);
1614 #endif
1615
1616 status = internal_read_partial_fits(seq->fptr[index], seq->ry, seq->bitpix, buffer, layer, area);
1617
1618 #ifdef _OPENMP
1619 omp_unset_lock(&seq->fd_lock[index]);
1620 #endif
1621 if (status)
1622 return 1;
1623
1624 flip_buffer(seq->bitpix, buffer, area);
1625 return 0;
1626 }
1627
siril_fits_compress(fits * f)1628 int siril_fits_compress(fits *f) {
1629 int status = 0;
1630 int comp_type = -1;
1631 siril_debug_print("Compressing FIT file with method %d and quantization %f\n",
1632 com.pref.comp.fits_method,
1633 com.pref.comp.fits_quantization);
1634 comp_type = get_compression_type(com.pref.comp.fits_method);
1635 siril_debug_print("cfitsio compression type %d\n",
1636 comp_type);
1637 if (comp_type < 0) {
1638 siril_log_message(_("Unknown FITS compression method in internal conversion\n"));
1639 return 1;
1640 }
1641
1642 /** In the doc: The default algorithm is called ``SUBTRACTIVE_DITHER_1''.
1643 * A second variation called ``SUBTRACTIVE_DITHER_2'' is also available,
1644 * which does the same thing except that any pixels with a value of 0.0 are not
1645 * dithered and instead the zero values are exactly preserved in the compressed image
1646 *
1647 * Here we need to use SUBTRACTIVE_DITHER in order to not have border artifacts at the
1648 * end of sracking with compressed images.
1649 */
1650 if (fits_set_quantize_dither(f->fptr, SUBTRACTIVE_DITHER_2, &status)) {
1651 report_fits_error(status);
1652 return 1;
1653 }
1654 status = 0;
1655
1656 if (fits_set_compression_type(f->fptr, comp_type, &status)) {
1657 report_fits_error(status);
1658 return 1;
1659 }
1660 status = 0;
1661
1662 if (fits_set_quantize_level(f->fptr, com.pref.comp.fits_quantization, &status)) {
1663 report_fits_error(status);
1664 return 1;
1665 }
1666
1667 status = 0;
1668
1669 /* Set the Hcompress scale factor if relevant */
1670 if (comp_type == HCOMPRESS_1) {
1671 if (fits_set_hcomp_scale(f->fptr, com.pref.comp.fits_hcompress_scale, &status)) {
1672 report_fits_error(status);
1673 return 1;
1674 }
1675 siril_debug_print("FITS HCompress scale factor %f\n",
1676 com.pref.comp.fits_hcompress_scale);
1677 status = 0;
1678 }
1679 return status;
1680 }
1681
1682 /* creates, saves and closes the file associated to f, overwriting previous */
savefits(const char * name,fits * f)1683 int savefits(const char *name, fits *f) {
1684 int status;
1685 char filename[256];
1686
1687 f->naxes[0] = f->rx;
1688 f->naxes[1] = f->ry;
1689
1690 if (f->naxis == 3 && f->naxes[2] != 3) {
1691 printf("Trying to save a FITS color file with more than 3 channels?");
1692 return 1;
1693 }
1694
1695 gboolean right_extension = FALSE;
1696 for (int i = 0; i < G_N_ELEMENTS(fit_extension); i++) {
1697 if (g_str_has_suffix(name, fit_extension[i])) {
1698 right_extension = TRUE;
1699 break;
1700 }
1701 }
1702
1703 if (!right_extension) {
1704 snprintf(filename, 255, "%s%s", name, com.pref.ext);
1705 } else {
1706 snprintf(filename, 255, "%s", name);
1707 }
1708
1709 g_unlink(filename); /* Delete old file if it already exists */
1710
1711 status = 0;
1712 if (siril_fits_create_diskfile(&(f->fptr), filename, &status)) { /* create new FITS file */
1713 report_fits_error(status);
1714 return 1;
1715 }
1716
1717 if (com.pref.comp.fits_enabled) {
1718 status = siril_fits_compress(f);
1719 if (status) {
1720 report_fits_error(status);
1721 return 1;
1722 }
1723 }
1724
1725 if (fits_create_img(f->fptr, f->bitpix, f->naxis, f->naxes, &status)) {
1726 report_fits_error(status);
1727 return 1;
1728 }
1729
1730 save_opened_fits(f);
1731
1732 status = 0;
1733 fits_close_file(f->fptr, &status);
1734 if (!status) {
1735 siril_log_message(_("Saving FITS: file %s, %ld layer(s), %ux%u pixels\n"),
1736 filename, f->naxes[2], f->rx, f->ry);
1737 }
1738 return 0;
1739 }
1740
save_opened_fits(fits * f)1741 int save_opened_fits(fits *f) {
1742 BYTE *data8;
1743 long orig[3] = { 1L, 1L, 1L };
1744 size_t i, pixel_count;
1745 int type, status = 0;
1746
1747 save_fits_header(f);
1748 pixel_count = f->naxes[0] * f->naxes[1] * f->naxes[2];
1749
1750 status = 0;
1751 switch (f->bitpix) {
1752 case BYTE_IMG:
1753 data8 = malloc(pixel_count * sizeof(BYTE));
1754 if (f->type == DATA_FLOAT) {
1755 for (i = 0; i < pixel_count; i++) {
1756 data8[i] = float_to_uchar_range(f->fdata[i]);
1757 }
1758 } else {
1759 double norm = get_normalized_value(f);
1760 for (i = 0; i < pixel_count; i++) {
1761 if (norm == USHRT_MAX_DOUBLE)
1762 data8[i] = conv_to_BYTE((double)f->data[i]);
1763 else
1764 data8[i] = round_to_BYTE(f->data[i]);
1765 }
1766 }
1767 if (fits_write_pix(f->fptr, TBYTE, orig, pixel_count, data8, &status)) {
1768 report_fits_error(status);
1769 free(data8);
1770 return 1;
1771 }
1772 f->lo >>= 8;
1773 f->hi >>= 8;
1774 free(data8);
1775 break;
1776 case SHORT_IMG:
1777 case USHORT_IMG:
1778 type = f->bitpix == SHORT_IMG ? TSHORT : TUSHORT;
1779 if (f->type == DATA_FLOAT) {
1780 WORD *data = float_buffer_to_ushort(f->fdata, f->naxes[0] * f->naxes[1] * f->naxes[2]);
1781 if (fits_write_pix(f->fptr, type, orig, pixel_count, data, &status)) {
1782 report_fits_error(status);
1783 free(data);
1784 return 1;
1785 }
1786 free(data);
1787 } else {
1788 if (f->orig_bitpix == BYTE_IMG) {
1789 conv_8_to_16(f->data, pixel_count);
1790 }
1791 if (fits_write_pix(f->fptr, type, orig, pixel_count, f->data,
1792 &status)) {
1793 report_fits_error(status);
1794 return 1;
1795 }
1796 }
1797 break;
1798 case FLOAT_IMG:
1799 if (f->type == DATA_USHORT) {
1800 if (f->orig_bitpix == BYTE_IMG) {
1801 conv_8_to_16(f->data, pixel_count);
1802 }
1803 f->fdata = malloc(pixel_count * sizeof(float));
1804 conv_16_to_32(f->data, f->fdata, pixel_count);
1805 fit_replace_buffer(f, f->fdata, DATA_FLOAT);
1806 }
1807 if (fits_write_pix(f->fptr, TFLOAT, orig, pixel_count, f->fdata, &status)) {
1808 report_fits_error(status);
1809 return 1;
1810 }
1811 break;
1812 case LONG_IMG:
1813 case LONGLONG_IMG:
1814 case DOUBLE_IMG:
1815 default:
1816 siril_log_message(_("ERROR: trying to save a FITS image "
1817 "with an unsupported format (%d).\n"), f->bitpix);
1818 fits_close_file(f->fptr, &status);
1819 return 1;
1820 }
1821
1822 if (!status) {
1823 // copy the entire header in memory
1824 if (f->header)
1825 free(f->header);
1826 f->header = copy_header(f);
1827 }
1828
1829 return 0;
1830 }
1831
1832 /* Duplicates some of a fits data into another, with various options; the third
1833 * parameter, oper, indicates with bits what operations will be done:
1834 *
1835 * - CP_ALLOC: allocates the to->data pointer to the size of from->data and
1836 * sets to->pdata; required if data is not already allocated with the
1837 * correct size or at all. No data is copied
1838 * - CP_INIT: initialize to->data with zeros, same size of the image in from,
1839 * but no other data is modified. Ignored if not used with CP_ALLOC.
1840 * - CP_COPYA: copies the actual data, from->data to to->data on all layers,
1841 * but no other information from the source. Should not be used with CP_INIT
1842 * - CP_FORMAT: copy all metadata and leaves data to null
1843 * - CP_EXPAND: forces the destination number of layers to be taken as 3, but
1844 * the other operations have no modifications, meaning that if the source
1845 * image has one layer, the output image will have only one actual layer
1846 * filled, and two filled with random data unless CP_INIT is used to fill it
1847 * with zeros.
1848 *
1849 * Example: to duplicate a fits from one to an unknown-allocated other, those
1850 * flags should be used: CP_ALLOC | CP_COPYA | CP_FORMAT
1851 *
1852 */
copyfits(fits * from,fits * to,unsigned char oper,int layer)1853 int copyfits(fits *from, fits *to, unsigned char oper, int layer) {
1854 int depth, i;
1855 size_t nbdata = from->naxes[0] * from->naxes[1];
1856
1857 if ((oper & CP_EXPAND))
1858 depth = 3;
1859 else depth = from->naxes[2];
1860
1861 if ((oper & CP_FORMAT)) {
1862 // copying metadata, not data or stats which are kept null
1863 memcpy(to, from, sizeof(fits));
1864 to->naxis = depth == 3 ? 3 : from->naxis;
1865 to->naxes[2] = depth;
1866 if (depth != from->naxes[2]) {
1867 to->maxi = -1.0;
1868 }
1869 to->stats = NULL;
1870 to->fptr = NULL;
1871 to->data = NULL;
1872 to->pdata[0] = NULL;
1873 to->pdata[1] = NULL;
1874 to->pdata[2] = NULL;
1875 to->fdata = NULL;
1876 to->fpdata[0] = NULL;
1877 to->fpdata[1] = NULL;
1878 to->fpdata[2] = NULL;
1879 to->header = NULL;
1880 to->history = NULL;
1881 to->date = NULL;
1882 to->date_obs = NULL;
1883 #ifdef HAVE_WCSLIB
1884 to->wcslib = NULL;
1885 #endif
1886 }
1887
1888 if ((oper & CP_ALLOC)) {
1889 // allocating to->data and assigning to->pdata
1890 if (from->type == DATA_USHORT) {
1891 WORD *olddata = to->data;
1892 if (!(to->data = realloc(to->data, nbdata * depth * sizeof(WORD)))) {
1893 PRINT_ALLOC_ERR;
1894 if (olddata)
1895 free(olddata);
1896 return -1;
1897 }
1898 to->type = DATA_USHORT;
1899 to->pdata[RLAYER] = to->data;
1900 if (depth == 3) {
1901 to->pdata[GLAYER] = to->data + nbdata;
1902 to->pdata[BLAYER] = to->data + 2 * nbdata;
1903 } else {
1904 to->pdata[GLAYER] = to->data;
1905 to->pdata[BLAYER] = to->data;
1906 }
1907
1908 if ((oper & CP_INIT)) {
1909 // clearing to->data allocated above
1910 memset(to->data, 0, nbdata * depth * sizeof(WORD));
1911 }
1912 }
1913 else if (from->type == DATA_FLOAT) {
1914 float *olddata = to->fdata;
1915 if (!(to->fdata = realloc(to->fdata, nbdata * depth * sizeof(float)))) {
1916 PRINT_ALLOC_ERR;
1917 if (olddata)
1918 free(olddata);
1919 return -1;
1920 }
1921 to->type = DATA_FLOAT;
1922 to->fpdata[RLAYER] = to->fdata;
1923 if (depth == 3) {
1924 to->fpdata[GLAYER] = to->fdata + nbdata;
1925 to->fpdata[BLAYER] = to->fdata + 2 * nbdata;
1926 } else {
1927 to->fpdata[GLAYER] = to->fdata;
1928 to->fpdata[BLAYER] = to->fdata;
1929 }
1930
1931 if ((oper & CP_INIT)) {
1932 // clearing to->fdata allocated above
1933 memset(to->fdata, 0, nbdata * depth * sizeof(WORD));
1934 }
1935 }
1936 else {
1937 fprintf(stderr, "unsupported copy\n");
1938 return -1;
1939 }
1940 }
1941
1942 if ((oper & CP_COPYA)) {
1943 // copying data
1944 if (to->type == DATA_USHORT)
1945 memcpy(to->data, from->data, nbdata * depth * sizeof(WORD));
1946 else if (to->type == DATA_FLOAT)
1947 memcpy(to->fdata, from->fdata, nbdata * depth * sizeof(float));
1948 else {
1949 fprintf(stderr, "unsupported copy\n");
1950 return -1;
1951 }
1952
1953 // copying stats
1954 if (from->stats) {
1955 for (i = 0; i < from->naxes[2]; i++) {
1956 if (from->stats[i])
1957 add_stats_to_fit(to, i, from->stats[i]);
1958 }
1959 } else {
1960 invalidate_stats_from_fit(to);
1961 }
1962 }
1963
1964 return 0;
1965 }
1966
extract_fits(fits * from,fits * to,int channel,gboolean to_float)1967 int extract_fits(fits *from, fits *to, int channel, gboolean to_float) {
1968 size_t nbdata = from->naxes[0] * from->naxes[1];
1969 // copying metadata, not data or stats which are kept null
1970 memcpy(to, from, sizeof(fits));
1971 to->naxis = 2;
1972 to->naxes[2] = 1L;
1973 to->maxi = -1.0;
1974 to->stats = NULL;
1975 to->fptr = NULL;
1976 to->header = NULL;
1977 to->history = NULL;
1978 to->date = NULL;
1979 to->date_obs = NULL;
1980 #ifdef HAVE_WCSLIB
1981 to->wcslib = NULL;
1982 #endif
1983
1984 if (from->type == DATA_USHORT)
1985 if (to_float) {
1986 /*if (from->bitpix == BYTE_IMG)
1987 to->fdata = ushort8_buffer_to_float(from->pdata[channel], nbdata);
1988 else to->fdata = ushort_buffer_to_float(from->pdata[channel], nbdata); */
1989 to->fdata = malloc(nbdata * sizeof(float));
1990 if (!to->fdata) {
1991 PRINT_ALLOC_ERR;
1992 return -1;
1993 }
1994 for (int i = 0; i < nbdata; i++)
1995 to->fdata[i] = (float)from->pdata[channel][i];
1996 to->type = DATA_FLOAT;
1997 to->bitpix = FLOAT_IMG;
1998 to->data = NULL;
1999 to->pdata[0] = NULL;
2000 to->pdata[1] = NULL;
2001 to->pdata[2] = NULL;
2002 to->fpdata[0] = to->fdata;
2003 to->fpdata[1] = to->fdata;
2004 to->fpdata[2] = to->fdata;
2005 }
2006 else {
2007 to->data = malloc(nbdata * sizeof(WORD));
2008 if (!to->data) {
2009 PRINT_ALLOC_ERR;
2010 return -1;
2011 }
2012 memcpy(to->data, from->pdata[channel], nbdata * sizeof(WORD));
2013 to->pdata[0] = to->data;
2014 to->pdata[1] = to->data;
2015 to->pdata[2] = to->data;
2016 }
2017 else if (from->type == DATA_FLOAT) {
2018 to->fdata = malloc(nbdata * sizeof(float));
2019 if (!to->fdata) {
2020 PRINT_ALLOC_ERR;
2021 return -1;
2022 }
2023 memcpy(to->fdata, from->fpdata[channel], nbdata * sizeof(float));
2024 to->fpdata[0] = to->fdata;
2025 to->fpdata[1] = to->fdata;
2026 to->fpdata[2] = to->fdata;
2027 }
2028 return 0;
2029 }
2030
2031 /* copy non-mandatory keywords from 'from' to 'to' */
copy_fits_metadata(fits * from,fits * to)2032 int copy_fits_metadata(fits *from, fits *to) {
2033 to->pixel_size_x = from->pixel_size_x;
2034 to->pixel_size_y = from->pixel_size_y;
2035 to->binning_x = from->binning_x;
2036 to->binning_y = from->binning_y;
2037
2038 if (from->date)
2039 to->date = g_date_time_ref(from->date);
2040 if (from->date_obs)
2041 to->date_obs = g_date_time_ref(from->date_obs);
2042 strncpy(to->instrume, from->instrume, FLEN_VALUE);
2043 strncpy(to->telescop, from->telescop, FLEN_VALUE);
2044 strncpy(to->observer, from->observer, FLEN_VALUE);
2045 strncpy(to->dft.type, from->dft.type, FLEN_VALUE);
2046 strncpy(to->dft.ord, from->dft.ord, FLEN_VALUE);
2047 strncpy(to->bayer_pattern, from->bayer_pattern, FLEN_VALUE);
2048 strncpy(to->row_order, from->row_order, FLEN_VALUE);
2049
2050 to->bayer_xoffset = from->bayer_xoffset;
2051 to->bayer_yoffset = from->bayer_yoffset;
2052 to->focal_length = from->focal_length;
2053 to->iso_speed = from->iso_speed;
2054 to->exposure = from->exposure;
2055 to->aperture = from->aperture;
2056 to->ccd_temp = from->ccd_temp;
2057 to->cvf = from->cvf;
2058 to->key_gain = from->key_gain;
2059 to->key_offset = from->key_offset;
2060 to->dft.norm[0] = from->dft.norm[0];
2061 to->dft.norm[1] = from->dft.norm[1];
2062 to->dft.norm[2] = from->dft.norm[2];
2063
2064 memcpy(&to->wcsdata, &from->wcsdata, sizeof(wcs_info));
2065
2066 return 0;
2067 }
2068
copy_fits_from_file(char * source,char * destination)2069 int copy_fits_from_file(char *source, char *destination) {
2070 fitsfile *infptr, *outfptr; /* FITS file pointers defined in fitsio.h */
2071 int status = 0; /* status must always be initialized = 0 */
2072
2073 /* Open the input file */
2074 if (!siril_fits_open_diskfile(&infptr, source, READONLY, &status)) {
2075 /* Create the output file */
2076 if (!siril_fits_create_diskfile(&outfptr, destination, &status)) {
2077
2078 /* copy the previous, current, and following HDUs */
2079 fits_copy_file(infptr, outfptr, 1, 1, 1, &status);
2080
2081 fits_close_file(outfptr, &status);
2082 }
2083 fits_close_file(infptr, &status);
2084 }
2085
2086 /* if error occured, print out error message */
2087 if (status)
2088 report_fits_error(status);
2089 return (status);
2090 }
2091
save1fits16(const char * filename,fits * fit,int layer)2092 int save1fits16(const char *filename, fits *fit, int layer) {
2093 if (layer != RLAYER) {
2094 size_t nbdata = fit->naxes[0] * fit->naxes[1];
2095 memcpy(fit->data, fit->data + layer * nbdata, nbdata * sizeof(WORD));
2096 }
2097 fit->naxis = 2;
2098 fit->naxes[2] = 1;
2099 return savefits(filename, fit);
2100 }
2101
save1fits32(const char * filename,fits * fit,int layer)2102 int save1fits32(const char *filename, fits *fit, int layer) {
2103 if (layer != RLAYER) {
2104 size_t nbdata = fit->naxes[0] * fit->naxes[1];
2105 memcpy(fit->fdata, fit->fdata + layer * nbdata, nbdata * sizeof(float));
2106 }
2107 fit->naxis = 2;
2108 fit->naxes[2] = 1;
2109 return savefits(filename, fit);
2110 }
2111
2112 /* this method converts 24-bit RGB or BGR data (no padding) to 48-bit FITS data.
2113 * order is RGB when inverted is FALSE, BGR when inverted is TRUE
2114 * fit->data has to be already allocated and fit->rx and fit->ry must be correct */
rgb24bit_to_fits48bit(unsigned char * rgbbuf,fits * fit,gboolean inverted)2115 void rgb24bit_to_fits48bit(unsigned char *rgbbuf, fits *fit, gboolean inverted) {
2116 size_t nbdata = fit->naxes[0] * fit->naxes[1];
2117 WORD *rdata, *gdata, *bdata;
2118 rdata = fit->pdata[RLAYER] = fit->data;
2119 gdata = fit->pdata[GLAYER] = fit->data + nbdata;
2120 bdata = fit->pdata[BLAYER] = fit->data + 2 * nbdata;
2121 for (size_t i = 0; i < nbdata; ++i) {
2122 if (inverted)
2123 *bdata++ = (WORD) *rgbbuf++;
2124 else
2125 *rdata++ = (WORD) *rgbbuf++;
2126 *gdata++ = (WORD) *rgbbuf++;
2127 if (inverted)
2128 *rdata++ = (WORD) *rgbbuf++;
2129 else
2130 *bdata++ = (WORD) *rgbbuf++;
2131 }
2132 }
2133
2134 /* this method converts 8-bit gray data to 16-bit FITS data.
2135 * fit->data has to be already allocated and fit->rx and fit->ry must be correct */
rgb8bit_to_fits16bit(unsigned char * graybuf,fits * fit)2136 void rgb8bit_to_fits16bit(unsigned char *graybuf, fits *fit) {
2137 WORD *data;
2138 size_t i, nbdata = fit->naxes[0] * fit->naxes[1];
2139 fit->pdata[0] = fit->data;
2140 fit->pdata[1] = fit->data;
2141 fit->pdata[2] = fit->data;
2142 data = fit->data;
2143 for (i = 0; i < nbdata; ++i) {
2144 *data++ = *graybuf++;
2145 }
2146 }
2147
2148 /* this method converts 48-bit RGB or BGR data (no padding) to 48-bit FITS data.
2149 * order is RGB when inverted is FALSE, BGR when inverted is TRUE
2150 * the endianness of the data, since we have two byte per value, may not match the endianness
2151 * of our FITS files, so the change_endian parameter allows to flip the endian.
2152 * fit->data has to be already allocated and fit->rx and fit->ry must be correct */
rgb48bit_to_fits48bit(WORD * rgbbuf,fits * fit,gboolean inverted,gboolean change_endian)2153 void rgb48bit_to_fits48bit(WORD *rgbbuf, fits *fit, gboolean inverted,
2154 gboolean change_endian) {
2155 size_t i, nbdata = fit->naxes[0] * fit->naxes[1];
2156 WORD *rdata, *gdata, *bdata, curval;
2157 rdata = fit->pdata[0] = fit->data;
2158 gdata = fit->pdata[1] = fit->data + nbdata;
2159 bdata = fit->pdata[2] = fit->data + 2 * nbdata;
2160 for (i = 0; i < nbdata; ++i) {
2161 curval = *rgbbuf++;
2162 if (change_endian)
2163 curval = change_endianness16(curval);
2164 if (inverted)
2165 *bdata++ = curval;
2166 else
2167 *rdata++ = curval;
2168
2169 curval = *rgbbuf++;
2170 if (change_endian)
2171 curval = change_endianness16(curval);
2172 *gdata++ = curval;
2173
2174 curval = *rgbbuf++;
2175 if (change_endian)
2176 curval = change_endianness16(curval);
2177 if (inverted)
2178 *rdata++ = curval;
2179 else
2180 *bdata++ = curval;
2181 }
2182 }
2183
2184 /* this method flips top-bottom of fit data.
2185 * fit->rx, fit->ry, fit->naxes[2] and fit->pdata[*] are required to be assigned correctly */
fits_flip_top_to_bottom_ushort(fits * fit)2186 static void fits_flip_top_to_bottom_ushort(fits *fit) {
2187 int line, axis, line_size;
2188 WORD *swapline, *src, *dst;
2189
2190 line_size = fit->rx * sizeof(WORD);
2191 swapline = malloc(line_size);
2192
2193 for (axis = 0; axis < fit->naxes[2]; axis++) {
2194 for (line = 0; line < fit->ry / 2; line++) {
2195 src = fit->pdata[axis] + line * fit->rx;
2196 dst = fit->pdata[axis] + (fit->ry - line - 1) * fit->rx;
2197
2198 memcpy(swapline, src, line_size);
2199 memcpy(src, dst, line_size);
2200 memcpy(dst, swapline, line_size);
2201 }
2202 }
2203 free(swapline);
2204 }
2205
fits_flip_top_to_bottom_float(fits * fit)2206 static void fits_flip_top_to_bottom_float(fits *fit) {
2207 int line, axis, line_size;
2208 float *swapline, *src, *dst;
2209
2210 line_size = fit->rx * sizeof(float);
2211 swapline = malloc(line_size);
2212
2213 for (axis = 0; axis < fit->naxes[2]; axis++) {
2214 for (line = 0; line < fit->ry / 2; line++) {
2215 src = fit->fpdata[axis] + line * fit->rx;
2216 dst = fit->fpdata[axis] + (fit->ry - line - 1) * fit->rx;
2217
2218 memcpy(swapline, src, line_size);
2219 memcpy(src, dst, line_size);
2220 memcpy(dst, swapline, line_size);
2221 }
2222 }
2223 free(swapline);
2224 }
2225
fits_flip_top_to_bottom(fits * fit)2226 void fits_flip_top_to_bottom(fits *fit) {
2227 if (fit->type == DATA_USHORT)
2228 fits_flip_top_to_bottom_ushort(fit);
2229 else if (fit->type == DATA_FLOAT)
2230 fits_flip_top_to_bottom_float(fit);
2231 }
2232
2233 /* This function copies an area from the fits 'from' on layer 'layer' into
2234 * another and initializes all relevant data */
2235 /* the crop function does the same but in place and for all channels without
2236 * reallocating */
extract_region_from_fits_ushort(fits * from,int layer,fits * to,const rectangle * area)2237 static void extract_region_from_fits_ushort(fits *from, int layer, fits *to,
2238 const rectangle *area) {
2239 int x, y, d, ystart, yend;
2240 clearfits(to);
2241 to->data = malloc(area->w * area->h * sizeof(WORD));
2242
2243 d = 0;
2244 ystart = from->ry - area->y - area->h;
2245 yend = from->ry - area->y;
2246 for (y = ystart; y < yend; y++) {
2247 for (x = area->x; x < area->x + area->w; x++) {
2248 to->data[d++] = from->pdata[layer][x + y * from->rx];
2249 }
2250 }
2251
2252 to->rx = area->w;
2253 to->ry = area->h;
2254 to->naxes[0] = area->w;
2255 to->naxes[1] = area->h;
2256 to->naxes[2] = 1;
2257 to->naxis = 2;
2258 to->pdata[0] = to->data;
2259 to->pdata[1] = to->data;
2260 to->pdata[2] = to->data;
2261 to->bitpix = from->bitpix;
2262 to->type = DATA_USHORT;
2263 }
2264
extract_region_from_fits_float(fits * from,int layer,fits * to,const rectangle * area)2265 static void extract_region_from_fits_float(fits *from, int layer, fits *to,
2266 const rectangle *area) {
2267 int x, y, d, ystart, yend;
2268 clearfits(to);
2269 to->fdata = malloc(area->w * area->h * sizeof(float));
2270
2271 d = 0;
2272 ystart = from->ry - area->y - area->h;
2273 yend = from->ry - area->y;
2274 for (y = ystart; y < yend; y++) {
2275 for (x = area->x; x < area->x + area->w; x++) {
2276 to->fdata[d++] = from->fpdata[layer][x + y * from->rx];
2277 }
2278 }
2279
2280 to->rx = area->w;
2281 to->ry = area->h;
2282 to->naxes[0] = area->w;
2283 to->naxes[1] = area->h;
2284 to->naxes[2] = 1;
2285 to->naxis = 2;
2286 to->fpdata[0] = to->fdata;
2287 to->fpdata[1] = to->fdata;
2288 to->fpdata[2] = to->fdata;
2289 to->bitpix = from->bitpix;
2290 to->type = DATA_FLOAT;
2291 }
2292
extract_region_from_fits(fits * from,int layer,fits * to,const rectangle * area)2293 void extract_region_from_fits(fits *from, int layer, fits *to,
2294 const rectangle *area) {
2295 if (from->type == DATA_USHORT)
2296 extract_region_from_fits_ushort(from, layer, to, area);
2297 else if (from->type == DATA_FLOAT)
2298 extract_region_from_fits_float(from, layer, to, area);
2299 }
2300
new_fit_image(fits ** fit,int width,int height,int nblayer,data_type type)2301 int new_fit_image(fits **fit, int width, int height, int nblayer, data_type type) {
2302 return new_fit_image_with_data(fit, width, height, nblayer, type, NULL);
2303 }
2304
2305 /* creates a new fit image from scratch (NULL fit) or into a fits * previously
2306 * allocated, with provided or non-cleared (random) data. */
new_fit_image_with_data(fits ** fit,int width,int height,int nblayer,data_type type,void * data)2307 int new_fit_image_with_data(fits **fit, int width, int height, int nblayer, data_type type, void *data) {
2308 size_t npixels, data_size;
2309 gboolean data_is_local = FALSE;
2310 g_assert(width > 0);
2311 g_assert(height > 0);
2312 g_assert(nblayer == 1 || nblayer == 3);
2313
2314 npixels = width * height;
2315 data_size = type == DATA_USHORT ? sizeof(WORD) : sizeof(float);
2316
2317 if (!data) {
2318 data = malloc(npixels * nblayer * data_size);
2319 if (!data) {
2320 PRINT_ALLOC_ERR;
2321 return -1;
2322 }
2323 data_is_local = TRUE;
2324 }
2325
2326 if (*fit)
2327 clearfits(*fit);
2328 else {
2329 *fit = calloc(1, sizeof(fits));
2330 if (!*fit) {
2331 PRINT_ALLOC_ERR;
2332 if (data_is_local)
2333 free(data);
2334 return -1;
2335 }
2336 }
2337
2338 if (nblayer > 1)
2339 (*fit)->naxis = 3;
2340 else (*fit)->naxis = 2;
2341 (*fit)->rx = width;
2342 (*fit)->ry = height;
2343 (*fit)->naxes[0] = width;
2344 (*fit)->naxes[1] = height;
2345 (*fit)->naxes[2] = nblayer;
2346 (*fit)->type = type;
2347
2348 if (type == DATA_USHORT) {
2349 (*fit)->bitpix = USHORT_IMG;
2350 (*fit)->orig_bitpix = USHORT_IMG;
2351 (*fit)->data = (WORD *)data;
2352 (*fit)->pdata[RLAYER] = (*fit)->data;
2353 if ((*fit)->naxis == 3) {
2354 (*fit)->pdata[GLAYER] = (*fit)->data + npixels;
2355 (*fit)->pdata[BLAYER] = (*fit)->data + npixels* 2;
2356 } else {
2357 (*fit)->pdata[GLAYER] = (*fit)->data;
2358 (*fit)->pdata[BLAYER] = (*fit)->data;
2359 }
2360 } else if (type == DATA_FLOAT) {
2361 (*fit)->bitpix = FLOAT_IMG;
2362 (*fit)->orig_bitpix = FLOAT_IMG;
2363 (*fit)->fdata = (float *)data;
2364 (*fit)->fpdata[RLAYER] = (*fit)->fdata;
2365 if ((*fit)->naxis == 3) {
2366 (*fit)->fpdata[GLAYER] = (*fit)->fdata + npixels;
2367 (*fit)->fpdata[BLAYER] = (*fit)->fdata + npixels * 2;
2368 } else {
2369 (*fit)->fpdata[GLAYER] = (*fit)->fdata;
2370 (*fit)->fpdata[BLAYER] = (*fit)->fdata;
2371 }
2372 }
2373 return 0;
2374 }
2375
fit_replace_buffer(fits * fit,void * newbuf,data_type newtype)2376 void fit_replace_buffer(fits *fit, void *newbuf, data_type newtype) {
2377 fit->type = newtype;
2378 invalidate_stats_from_fit(fit);
2379
2380 size_t nbdata = fit->naxes[0] * fit->naxes[1];
2381 if (newtype == DATA_USHORT) {
2382 fit->bitpix = USHORT_IMG;
2383 fit->orig_bitpix = USHORT_IMG;
2384 fit->data = (WORD *)newbuf;
2385 fit->pdata[RLAYER] = fit->data;
2386 if (fit->naxis == 3) {
2387 fit->pdata[GLAYER] = fit->data + nbdata;
2388 fit->pdata[BLAYER] = fit->data + nbdata * 2;
2389 } else {
2390 fit->pdata[GLAYER] = fit->data;
2391 fit->pdata[BLAYER] = fit->data;
2392 }
2393 if (fit->fdata) {
2394 free(fit->fdata);
2395 fit->fdata = NULL;
2396 }
2397 fit->fpdata[0] = NULL;
2398 fit->fpdata[1] = NULL;
2399 fit->fpdata[2] = NULL;
2400 siril_debug_print("Changed a fit data (WORD)\n");
2401 } else if (newtype == DATA_FLOAT) {
2402 fit->bitpix = FLOAT_IMG;
2403 fit->orig_bitpix = FLOAT_IMG;
2404 fit->fdata = (float *)newbuf;
2405 fit->fpdata[RLAYER] = fit->fdata;
2406 if (fit->naxis == 3) {
2407 fit->fpdata[GLAYER] = fit->fdata + nbdata;
2408 fit->fpdata[BLAYER] = fit->fdata + nbdata * 2;
2409 } else {
2410 fit->fpdata[GLAYER] = fit->fdata;
2411 fit->fpdata[BLAYER] = fit->fdata;
2412 }
2413 if (fit->data) {
2414 free(fit->data);
2415 fit->data = NULL;
2416 }
2417 fit->pdata[0] = NULL;
2418 fit->pdata[1] = NULL;
2419 fit->pdata[2] = NULL;
2420 siril_debug_print("Changed a fit data (FLOAT)\n");
2421 }
2422 }
2423
fit_debayer_buffer(fits * fit,void * newbuf)2424 void fit_debayer_buffer(fits *fit, void *newbuf) {
2425 size_t nbdata = fit->naxes[0] * fit->naxes[1];
2426
2427 /* before changing naxis, we clear fit->stats that was allocated for
2428 * one channel */
2429 full_stats_invalidation_from_fit(fit);
2430
2431 fit->naxis = 3;
2432 fit->naxes[2] = 3;
2433 if (fit->type == DATA_USHORT) {
2434 if (fit->data)
2435 free(fit->data);
2436 fit->data = (WORD *)newbuf;
2437 fit->pdata[RLAYER] = fit->data;
2438 fit->pdata[GLAYER] = fit->data + nbdata;
2439 fit->pdata[BLAYER] = fit->data + nbdata * 2;
2440 }
2441 else if (fit->type == DATA_FLOAT) {
2442 if (fit->fdata)
2443 free(fit->fdata);
2444 fit->fdata = (float *)newbuf;
2445 fit->fpdata[RLAYER] = fit->fdata;
2446 fit->fpdata[GLAYER] = fit->fdata + nbdata;
2447 fit->fpdata[BLAYER] = fit->fdata + nbdata * 2;
2448 }
2449 }
2450
gray2rgb(float gray,guchar * rgb)2451 static void gray2rgb(float gray, guchar *rgb) {
2452 *rgb++ = (guchar) round_to_BYTE(255. * gray);
2453 *rgb++ = (guchar) round_to_BYTE(255. * gray);
2454 *rgb++ = (guchar) round_to_BYTE(255. * gray);
2455 }
2456
free_preview_data(guchar * pixels,gpointer data)2457 static GdkPixbufDestroyNotify free_preview_data(guchar *pixels, gpointer data) {
2458 free(pixels);
2459 free(data);
2460 return FALSE;
2461 }
2462
logviz(double arg)2463 static double logviz(double arg) {
2464 return log1p(arg);
2465 } // for PREVIEW_LOG
2466
2467 #define TRYFITS(f, ...) \
2468 do{ \
2469 status = FALSE; \
2470 f(__VA_ARGS__, &status); \
2471 if(status){ \
2472 free(ima_data); \
2473 fits_close_file(fp, &status); \
2474 return NULL; \
2475 } \
2476 } while(0)
2477
2478 /**
2479 * Create a monochrome preview of a FITS file in a GdkPixbuf
2480 * @param filename
2481 * @return a GdkPixbuf containing the preview or NULL
2482 */
get_thumbnail_from_fits(char * filename,gchar ** descr)2483 GdkPixbuf* get_thumbnail_from_fits(char *filename, gchar **descr) {
2484 fitsfile *fp;
2485 gchar *description;
2486 const int MAX_SIZE = com.pref.thumbnail_size;
2487 float nullval = 0.;
2488 int naxis, dtype, stat, status, frames;
2489
2490 long naxes[4];
2491 float *ima_data = NULL;
2492
2493 TRYFITS(fits_open_diskfile, &fp, filename, READONLY);
2494
2495 if (siril_fits_move_first_image(fp)) {
2496 siril_log_message(_("Selecting the primary header failed, is the FITS file '%s' malformed?\n"), filename);
2497 return NULL;
2498 }
2499
2500 TRYFITS(fits_get_img_param, fp, 4, &dtype, &naxis, naxes);
2501
2502 const int w = naxes[0];
2503 const int h = naxes[1];
2504 size_t sz = w * h;
2505 ima_data = malloc(sz * sizeof(float));
2506
2507 TRYFITS(fits_read_img, fp, TFLOAT, 1, sz, &nullval, ima_data, &stat);
2508
2509 const int x = (int) ceil((float) w / MAX_SIZE);
2510 const int y = (int) ceil((float) h / MAX_SIZE);
2511 const int pixScale = (x > y) ? x : y; // picture scale factor
2512 const int Ws = w / pixScale; // picture width in pixScale blocks
2513 const int Hs = h / pixScale; // -//- height pixScale
2514
2515 const int n_channels = naxis == 3 ? naxis : 1;
2516
2517 if (fitseq_is_fitseq(filename, &frames)) { // FIXME: we reopen the file in this function
2518 description = g_strdup_printf("%d x %d %s\n%d %s (%d bits)\n%d %s\n%s", w,
2519 h, ngettext("pixel", "pixels", h), n_channels,
2520 ngettext("channel", "channels", n_channels), abs(dtype), frames,
2521 ngettext("frame", "frames", frames), _("(Monochrome Preview)"));
2522 } else {
2523 description = g_strdup_printf("%d x %d %s\n%d %s (%d bits)\n%s", w,
2524 h, ngettext("pixel", "pixels", h), n_channels,
2525 ngettext("channel", "channels", n_channels), abs(dtype),
2526 _("(Monochrome Preview)"));
2527 }
2528
2529 #ifdef _OPENMP
2530 #pragma omp parallel
2531 #endif
2532 {
2533 // array for preview picture line
2534 float pix[MAX_SIZE];
2535 #ifdef _OPENMP
2536 #pragma omp for
2537 #endif
2538 for (int i = 0; i < Hs; i++) { // cycle through a blocks by lines
2539 int M = i * pixScale;
2540 for (int j = 0; j < MAX_SIZE; j++) { // zero line buffer
2541 pix[j] = 0;
2542 }
2543 unsigned int m = 0; // amount of strings read in block
2544 for (int l = 0; l < pixScale; l++, m++) { // cycle through a block lines
2545 float *ptr = &ima_data[M * w];
2546 int N = 0; // number of column
2547 for (int j = 0; j < Ws; j++) { // cycle through a blocks by columns
2548 unsigned int n = 0; // amount of columns read in block
2549 float sum = 0.f; // average intensity in block
2550 for (int k = 0; k < pixScale; k++, n++) { // cycle through block pixels
2551 if (N++ < w) // row didn't end
2552 sum += *ptr++; // sum[(pix-min)/wd]/n = [sum(pix)/n-min]/wd
2553 else
2554 break;
2555 }
2556 pix[j] += sum / n; //(byte / n - min)/wd;
2557 }
2558 if (++M >= h)
2559 break;
2560 }
2561 // fill unused picture pixels
2562 float *ptr = &ima_data[i * Ws];
2563 for (int l = 0; l < Ws; l++)
2564 *ptr++ = pix[l] / m;
2565 }
2566 }
2567
2568 float *ptr = ima_data;
2569 sz = Ws * Hs;
2570 float max = *ptr;
2571 float min = max;
2572 float avr = 0.f;
2573 for (size_t i = 0; i < sz; i++, ptr++) {
2574 const float val = *ptr;
2575 max = max(max, val);
2576 min = min(min, val);
2577 avr += val;
2578 }
2579 avr /= (float) sz;
2580
2581 /* use FITS keyword if available for a better visualization */
2582 float lo = 0.f;
2583 float hi = 0.f;
2584 __tryToFindKeywords(fp, TFLOAT, MIPSLO, &lo);
2585 __tryToFindKeywords(fp, TFLOAT, MIPSHI, &hi);
2586
2587 if (hi != lo && hi != 0.f && abs(dtype) <= USHORT_IMG) {
2588 min = lo;
2589 max = hi;
2590 } else if (dtype <= FLOAT_IMG) { // means float or double image
2591 WORD wlo, whi;
2592 if (!try_read_float_lo_hi(fp, &wlo, &whi)) {
2593 min = (float) wlo / USHRT_MAX_SINGLE;
2594 max = (float) whi / USHRT_MAX_SINGLE;
2595 }
2596 }
2597
2598 float wd = max - min;
2599 avr = (avr - min) / wd; // normal average by preview
2600 avr = -logf(avr); // scale factor
2601 if (avr > 1.) {
2602 wd /= avr;
2603 }
2604
2605 guchar *pixbuf_data = malloc(3 * MAX_SIZE * MAX_SIZE * sizeof(guchar));
2606
2607 #ifdef _OPENMP
2608 #pragma omp parallel for num_threads(com.max_thread)
2609 #endif
2610 for (int i = Hs - 1; i > -1; i--) { // fill pixbuf mirroring image by vertical
2611 guchar *pptr = &pixbuf_data[Ws * i * 3];
2612 float *p = &ima_data[(Hs - i - 1) * Ws];
2613 for (int j = 0; j < Ws; j++) {
2614 gray2rgb(logviz((*p++ - min) / wd), pptr);
2615 pptr += 3;
2616 }
2617 }
2618 fits_close_file(fp, &status);
2619 free(ima_data);
2620 GdkPixbuf *pixbuf = gdk_pixbuf_new_from_data(pixbuf_data, // guchar* data
2621 GDK_COLORSPACE_RGB, // only this supported
2622 FALSE, // no alpha
2623 8, // number of bits
2624 Ws, Hs, // size
2625 Ws * 3, // line length in bytes
2626 (GdkPixbufDestroyNotify) free_preview_data, // function (*GdkPixbufDestroyNotify) (guchar *pixels, gpointer data);
2627 NULL);
2628 *descr = description;
2629 return pixbuf;
2630 }
2631
2632