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 <stdlib.h>
22 #include <stdio.h>
23 #include <assert.h>
24 #include "core/siril.h"
25 #include "core/proto.h"
26 #include "core/processing.h"
27 #include "io/sequence.h"
28 #include "io/ser.h"
29 #include "io/image_format_fits.h"
30 #include "stacking.h"
31 #include "gui/progress_and_log.h"
32
33 struct sum_stacking_data {
34 guint64 *sum[3]; // the new image's channels
35 double *fsum[3]; // the new image's channels, for float input image
36 double exposure; // sum of the exposures
37 int reglayer; // layer used for registration data
38 int ref_image; // reference image index in the stacked sequence
39 gboolean input_32bits; // input is a sequence of 32-bit float images
40 gboolean output_32bits; // output a 32-bit float image instead of the default ushort
41 };
42
sum_stacking_prepare_hook(struct generic_seq_args * args)43 static int sum_stacking_prepare_hook(struct generic_seq_args *args) {
44 struct sum_stacking_data *ssdata = args->user;
45 size_t nbdata = args->seq->ry * args->seq->rx;
46
47 if (ssdata->input_32bits) {
48 ssdata->fsum[0] = calloc(nbdata, sizeof(double) * args->seq->nb_layers);
49 if (ssdata->fsum[0] == NULL){
50 PRINT_ALLOC_ERR;
51 return ST_ALLOC_ERROR;
52 }
53 if(args->seq->nb_layers == 3){
54 ssdata->fsum[1] = ssdata->fsum[0] + nbdata;
55 ssdata->fsum[2] = ssdata->fsum[0] + nbdata*2;
56 } else {
57 ssdata->fsum[1] = NULL;
58 ssdata->fsum[2] = NULL;
59 }
60 ssdata->sum[0] = NULL;
61 } else {
62 ssdata->sum[0] = calloc(nbdata, sizeof(guint64) * args->seq->nb_layers);
63 if (ssdata->sum[0] == NULL){
64 PRINT_ALLOC_ERR;
65 return ST_ALLOC_ERROR;
66 }
67 if(args->seq->nb_layers == 3){
68 ssdata->sum[1] = ssdata->sum[0] + nbdata; // index of green layer in sum[0]
69 ssdata->sum[2] = ssdata->sum[0] + nbdata*2; // index of blue layer in sum[0]
70 } else {
71 ssdata->sum[1] = NULL;
72 ssdata->sum[2] = NULL;
73 }
74 ssdata->fsum[0] = NULL;
75 }
76
77 ssdata->exposure = 0.0;
78 return ST_OK;
79 }
80
sum_stacking_image_hook(struct generic_seq_args * args,int o,int i,fits * fit,rectangle * _)81 static int sum_stacking_image_hook(struct generic_seq_args *args, int o, int i, fits *fit, rectangle *_) {
82 struct sum_stacking_data *ssdata = args->user;
83 int shiftx, shifty, nx, ny, x, y, layer;
84 size_t ii, pixel = 0; // index in sum[0]
85
86 #ifdef _OPENMP
87 #pragma omp atomic
88 #endif
89 ssdata->exposure += fit->exposure;
90
91 if (ssdata->reglayer != -1 && args->seq->regparam[ssdata->reglayer]) {
92 shiftx = round_to_int(args->seq->regparam[ssdata->reglayer][i].shiftx * (float)args->seq->upscale_at_stacking);
93 shifty = round_to_int(args->seq->regparam[ssdata->reglayer][i].shifty * (float)args->seq->upscale_at_stacking);
94 } else {
95 shiftx = 0;
96 shifty = 0;
97 }
98
99 for (y = 0; y < fit->ry; ++y) {
100 for (x = 0; x < fit->rx; ++x) {
101 nx = x - shiftx;
102 ny = y - shifty;
103 if (nx >= 0 && nx < fit->rx && ny >= 0 && ny < fit->ry) {
104 // we have data for this pixel
105 ii = ny * fit->rx + nx; // index in source image
106 if (ii < fit->rx * fit->ry) {
107 for (layer = 0; layer < args->seq->nb_layers; ++layer) {
108 if (ssdata->input_32bits) {
109 #ifdef _OPENMP
110 #pragma omp atomic
111 #endif
112 ssdata->fsum[layer][pixel] += (double)fit->fpdata[layer][ii];
113 }
114 else
115 #ifdef _OPENMP
116 #pragma omp atomic
117 #endif
118 ssdata->sum[layer][pixel] += fit->pdata[layer][ii];
119 }
120 }
121 }
122 ++pixel;
123 }
124 }
125 return ST_OK;
126 }
127
128 // convert the result and store it into gfit
sum_stacking_finalize_hook(struct generic_seq_args * args)129 static int sum_stacking_finalize_hook(struct generic_seq_args *args) {
130 struct sum_stacking_data *ssdata = args->user;
131 guint64 max = 0L; // max value of the image's channels
132 double fmax = 0.0;
133 size_t i, nbdata;
134 int layer;
135
136 nbdata = args->seq->ry * args->seq->rx * args->seq->nb_layers;
137 // find the max first
138 if (ssdata->input_32bits) {
139 #ifdef _OPENMP
140 #pragma omp parallel for reduction(max:fmax)
141 #endif
142 for (i = 0; i < nbdata; ++i) {
143 if (ssdata->fsum[0][i] > fmax)
144 fmax = ssdata->fsum[0][i];
145 }
146 } else {
147 #ifdef _OPENMP
148 #pragma omp parallel for reduction(max:max)
149 #endif
150 for (i = 0; i < nbdata; ++i)
151 if (ssdata->sum[0][i] > max)
152 max = ssdata->sum[0][i];
153 }
154
155 clearfits(&gfit);
156 fits *fit = &gfit;
157 if (new_fit_image(&fit, args->seq->rx, args->seq->ry, args->seq->nb_layers, ssdata->output_32bits ? DATA_FLOAT : DATA_USHORT))
158 return ST_GENERIC_ERROR;
159
160 /* We copy metadata from reference to the final fit */
161 int ref = ssdata->ref_image;
162 if (args->seq->type == SEQ_REGULAR) {
163 if (!seq_open_image(args->seq, ref)) {
164 import_metadata_from_fitsfile(args->seq->fptr[ref], &gfit);
165 seq_close_image(args->seq, ref);
166 }
167 } else if (args->seq->type == SEQ_FITSEQ) {
168 if (!fitseq_set_current_frame(args->seq->fitseq_file, ref))
169 import_metadata_from_fitsfile(args->seq->fitseq_file->fptr, &gfit);
170 } else if (args->seq->type == SEQ_SER) {
171 import_metadata_from_serfile(args->seq->ser_file, &gfit);
172 }
173
174 gfit.exposure = ssdata->exposure;
175 nbdata = args->seq->ry * args->seq->rx;
176
177 if (ssdata->output_32bits) {
178 if (ssdata->input_32bits) {
179 double ratio = 1.0 / fmax;
180 for (layer=0; layer<args->seq->nb_layers; ++layer){
181 double *from = ssdata->fsum[layer];
182 float *to = gfit.fpdata[layer];
183 for (i=0; i < nbdata; ++i) {
184 *to++ = (float)((double)(*from++) * ratio);
185 }
186 }
187 } else {
188 double ratio = 1.0 / (double)max;
189 for (layer=0; layer<args->seq->nb_layers; ++layer){
190 guint64 *from = ssdata->sum[layer];
191 float *to = gfit.fpdata[layer];
192 for (i=0; i < nbdata; ++i) {
193 *to++ = (float)((double)(*from++) * ratio);
194 }
195 }
196 }
197 } else {
198 double ratio = 1.0;
199 if (max > USHRT_MAX) {
200 ratio = USHRT_MAX_DOUBLE / (double)max;
201 siril_log_color_message(_("Reducing the stacking output to a 16-bit image will result in precision loss\n"), "salmon");
202 }
203
204 for (layer=0; layer<args->seq->nb_layers; ++layer){
205 guint64 *from = ssdata->sum[layer];
206 WORD *to = gfit.pdata[layer];
207 for (i=0; i < nbdata; ++i) {
208 if (ratio == 1.0)
209 *to++ = round_to_WORD(*from++);
210 else *to++ = round_to_WORD((double)(*from++) * ratio);
211 }
212 }
213 }
214
215 if (ssdata->sum[0]) free(ssdata->sum[0]);
216 if (ssdata->fsum[0]) free(ssdata->fsum[0]);
217 free(ssdata);
218 args->user = NULL;
219
220 return ST_OK;
221 }
222
stack_summing_generic(struct stacking_args * stackargs)223 int stack_summing_generic(struct stacking_args *stackargs) {
224 struct generic_seq_args *args = create_default_seqargs(stackargs->seq);
225 args->filtering_criterion = stackargs->filtering_criterion;
226 args->filtering_parameter = stackargs->filtering_parameter;
227 args->nb_filtered_images = stackargs->nb_images_to_stack;
228 args->prepare_hook = sum_stacking_prepare_hook;
229 args->image_hook = sum_stacking_image_hook;
230 args->finalize_hook = sum_stacking_finalize_hook;
231 args->description = _("Sum stacking");
232 args->already_in_a_thread = TRUE;
233
234 struct sum_stacking_data *ssdata = malloc(sizeof(struct sum_stacking_data));
235 ssdata->reglayer = stackargs->reglayer;
236 ssdata->ref_image = stackargs->ref_image;
237 assert(ssdata->ref_image >= 0 && ssdata->ref_image < args->seq->number);
238 ssdata->input_32bits = get_data_type(args->seq->bitpix) == DATA_FLOAT;
239 ssdata->output_32bits = stackargs->use_32bit_output;
240 if (ssdata->input_32bits)
241 assert(ssdata->output_32bits);
242 args->user = ssdata;
243
244 generic_sequence_worker(args);
245 return args->retval;
246 }
247
248