1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
|
/*
* This file is part of Siril, an astronomy image processor.
* Copyright (C) 2005-2011 Francois Meyer (dulle at free.fr)
* Copyright (C) 2012-2019 team free-astro (see more in AUTHORS file)
* Reference site is https://free-astro.org/index.php/Siril
*
* Siril is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Siril is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Siril. If not, see <http://www.gnu.org/licenses/>.
*/
#include "core/siril.h"
#include "core/processing.h"
#include "core/proto.h" // FITS functions
#include "io/sequence.h"
#include "stacking.h"
struct sum_stacking_data {
unsigned long *sum[3]; // the new image's channels
double exposure; // sum of the exposures
int reglayer; // layer used for registration data
};
static int sum_stacking_prepare_hook(struct generic_seq_args *args) {
struct sum_stacking_data *ssdata = args->user;
unsigned int nbdata = args->seq->ry * args->seq->rx;
ssdata->sum[0] = calloc(nbdata, sizeof(unsigned long)*args->seq->nb_layers);
if (ssdata->sum[0] == NULL){
printf("Stacking: memory allocation failure\n");
return -1;
}
if(args->seq->nb_layers == 3){
ssdata->sum[1] = ssdata->sum[0] + nbdata; // index of green layer in sum[0]
ssdata->sum[2] = ssdata->sum[0] + nbdata*2; // index of blue layer in sum[0]
} else {
ssdata->sum[1] = NULL;
ssdata->sum[2] = NULL;
}
ssdata->exposure = 0.0;
return 0;
}
static int sum_stacking_image_hook(struct generic_seq_args *args, int o, int i, fits *fit, rectangle *_) {
struct sum_stacking_data *ssdata = args->user;
int shiftx, shifty, nx, ny, x, y, ii, layer;
int pixel = 0; // index in sum[0]
#ifdef _OPENMP
#pragma omp atomic
#endif
ssdata->exposure += fit->exposure;
if (ssdata->reglayer != -1 && args->seq->regparam[ssdata->reglayer]) {
shiftx = roundf_to_int(args->seq->regparam[ssdata->reglayer][i].shiftx * args->seq->upscale_at_stacking);
shifty = roundf_to_int(args->seq->regparam[ssdata->reglayer][i].shifty * args->seq->upscale_at_stacking);
} else {
shiftx = 0;
shifty = 0;
}
for (y=0; y < fit->ry; ++y){
for (x=0; x < fit->rx; ++x){
nx = x - shiftx;
ny = y - shifty;
if (nx >= 0 && nx < fit->rx && ny >= 0 && ny < fit->ry) {
// we have data for this pixel
ii = ny * fit->rx + nx; // index in source image
if (ii >= 0 && ii < fit->rx * fit->ry){
for(layer=0; layer<args->seq->nb_layers; ++layer) {
#ifdef _OPENMP
#pragma omp atomic
#endif
ssdata->sum[layer][pixel] += fit->pdata[layer][ii];
}
}
}
++pixel;
}
}
return 0;
}
// convert the result and store it into gfit
static int sum_stacking_finalize_hook(struct generic_seq_args *args) {
struct sum_stacking_data *ssdata = args->user;
unsigned long max = 0L; // max value of the image's channels
unsigned int i, nbdata;
int layer;
nbdata = args->seq->ry * args->seq->rx * args->seq->nb_layers;
// find the max first
#pragma omp parallel for reduction(max:max)
for (i=0; i < nbdata; ++i)
if (ssdata->sum[0][i] > max)
max = ssdata->sum[0][i];
clearfits(&gfit);
fits *fit = &gfit;
if (new_fit_image(&fit, args->seq->rx, args->seq->ry, args->seq->nb_layers))
return -1;
/* We copy metadata from reference to the final fit */
if (args->seq->type == SEQ_REGULAR) {
int ref = 0;
// FIXME: sum stacking does not have a reference, should we add
// it or use another image here?
if (args->seq->reference_image > 0)
ref = args->seq->reference_image;
if (!seq_open_image(args->seq, ref)) {
import_metadata_from_fitsfile(args->seq->fptr[ref], &gfit);
seq_close_image(args->seq, ref);
}
}
gfit.hi = round_to_WORD(max);
gfit.exposure = ssdata->exposure;
gfit.bitpix = gfit.orig_bitpix = USHORT_IMG;
double ratio = 1.0;
if (max > USHRT_MAX)
ratio = USHRT_MAX_DOUBLE / (double)max;
nbdata = args->seq->ry * args->seq->rx;
for (layer=0; layer<args->seq->nb_layers; ++layer){
unsigned long* from = ssdata->sum[layer];
WORD *to = gfit.pdata[layer];
for (i=0; i < nbdata; ++i) {
if (ratio == 1.0)
*to++ = round_to_WORD(*from++);
else *to++ = round_to_WORD((double)(*from++) * ratio);
}
}
free(ssdata->sum[0]);
return 0;
}
int stack_summing_generic(struct stacking_args *stackargs) {
struct generic_seq_args *args = malloc(sizeof(struct generic_seq_args));
args->seq = stackargs->seq;
args->partial_image = FALSE;
args->filtering_criterion = stackargs->filtering_criterion;
args->filtering_parameter = stackargs->filtering_parameter;
args->nb_filtered_images = stackargs->nb_images_to_stack;
args->prepare_hook = sum_stacking_prepare_hook;
args->image_hook = sum_stacking_image_hook;
args->save_hook = NULL;
args->finalize_hook = sum_stacking_finalize_hook;
args->idle_function = NULL;
args->stop_on_error = TRUE;
args->description = _("Sum stacking");
args->has_output = FALSE;
args->already_in_a_thread = TRUE;
args->parallel = TRUE;
struct sum_stacking_data *ssdata = malloc(sizeof(struct sum_stacking_data));
ssdata->reglayer = stackargs->reglayer;
args->user = ssdata;
//start_in_new_thread(generic_sequence_worker, args);
generic_sequence_worker(args);
return args->retval;
}
|