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 <math.h>
22 #include <gsl/gsl_matrix.h>
23 #include <string.h>
24
25 #include "core/siril.h"
26 #include "core/proto.h"
27 #include "algos/PSF.h"
28 #include "algos/photometry.h"
29
30 #define hampel_a 1.7
31 #define hampel_b 3.4
32 #define hampel_c 8.5
33 #define sign(x,y) ((y)>=0?fabs(x):-fabs(x))
34 #define epsilon(x) 0.00000001
35 #define maxit 50
36 #define min_sky 5
37
hampel(double x)38 static double hampel(double x) {
39 if (x >= 0) {
40 if (x < hampel_a)
41 return x;
42 if (x < hampel_b)
43 return hampel_a;
44 if (x < hampel_c)
45 return hampel_a * (x - hampel_c) / (hampel_b - hampel_c);
46 } else {
47 if (x > -hampel_a)
48 return x;
49 if (x > -hampel_b)
50 return -hampel_a;
51 if (x > -hampel_c)
52 return hampel_a * (x + hampel_c) / (hampel_b - hampel_c);
53 }
54 return 0.0;
55 }
56
dhampel(double x)57 static double dhampel(double x) {
58 if (x >= 0) {
59 if (x < hampel_a)
60 return 1;
61 if (x < hampel_b)
62 return 0;
63 if (x < hampel_c)
64 return hampel_a / (hampel_b - hampel_c);
65 } else {
66 if (x > -hampel_a)
67 return 1;
68 if (x > -hampel_b)
69 return 0;
70 if (x > -hampel_c)
71 return -hampel_a / (hampel_b - hampel_c);
72 }
73 return 0.0;
74 }
75
qmedD(int n,double * a)76 static double qmedD(int n, double *a)
77 /* Vypocet medianu algoritmem Quick Median (Wirth) */
78 {
79 double w, x;
80 int i, j;
81 int k = ((n & 1) ? (n / 2) : ((n / 2) - 1));
82 int l = 0;
83 int r = n - 1;
84
85 while (l < r) {
86 x = a[k];
87 i = l;
88 j = r;
89 do {
90 while (a[i] < x)
91 i++;
92 while (x < a[j])
93 j--;
94 if (i <= j) {
95 w = a[i];
96 a[i] = a[j];
97 a[j] = w;
98 i++;
99 j--;
100 }
101 } while (i <= j);
102 if (j < k)
103 l = i;
104 if (k < i)
105 r = j;
106 }
107 return a[k];
108 }
109
robustmean(int n,double * x,double * mean,double * stdev)110 static int robustmean(int n, double *x, double *mean, double *stdev)
111 /* Newton's iterations */
112 {
113 int i, it;
114 double a, c, d, dt, r, s, sum1, sum2, sum3, psir;
115 double *xx;
116
117 if (n < 1) {
118 if (mean)
119 *mean = 0.0; /* a few data */
120 if (stdev)
121 *stdev = -1.0;
122 return 1;
123 }
124 if (n == 1) { /* only one point, but correct case */
125 if (mean)
126 *mean = x[0];
127 if (stdev)
128 *stdev = 0.0;
129 return 0;
130 }
131
132 /* initial values:
133 - median is the first approximation of location
134 - MAD/0.6745 is the first approximation of scale */
135 xx = malloc(n * sizeof(double));
136 if (!xx) {
137 PRINT_ALLOC_ERR;
138 return 1;
139 }
140 memcpy(xx, x, n * sizeof(double));
141 a = qmedD(n, xx);
142 for (i = 0; i < n; i++)
143 xx[i] = fabs(x[i] - a);
144 s = qmedD(n, xx) / 0.6745;
145 free(xx);
146
147 /* almost identical points on input */
148 if (fabs(s) < epsilon(s)) {
149 if (mean)
150 *mean = a;
151 if (stdev) {
152 double sum = 0.0;
153 for (i = 0; i < n; i++)
154 sum += (x[i] - a) * (x[i] - a);
155 *stdev = sqrt(sum / n);
156 }
157 return 0;
158 }
159
160 /* corrector's estimation */
161 dt = 0;
162 c = s * s * n * n / (n - 1);
163 for (it = 1; it <= maxit; it++) {
164 sum1 = sum2 = sum3 = 0.0;
165 for (i = 0; i < n; i++) {
166 r = (x[i] - a) / s;
167 psir = hampel(r);
168 sum1 += psir;
169 sum2 += dhampel(r);
170 sum3 += psir * psir;
171 }
172 if (fabs(sum2) < epsilon(sum2))
173 break;
174 d = s * sum1 / sum2;
175 a = a + d;
176 dt = c * sum3 / (sum2 * sum2);
177 if ((it > 2) && ((d * d < 1e-4 * dt) || (fabs(d) < 10.0 * epsilon(d))))
178 break;
179 }
180 if (mean)
181 *mean = a;
182 if (stdev)
183 *stdev = (dt > 0 ? sqrt(dt) : 0);
184 return 0;
185 }
186
phot_alloc()187 static photometry *phot_alloc() {
188 photometry *phot;
189
190 phot = (photometry*) calloc(sizeof(photometry), 1);
191 if (phot == NULL) {
192 PRINT_ALLOC_ERR;
193 }
194 return phot;
195 }
196
getMagnitude(double intensity)197 static double getMagnitude(double intensity) {
198 return -2.5 * log10(intensity);
199 }
200
getCameraGain()201 static double getCameraGain() {
202 if (gfit.type == DATA_FLOAT) {
203 return com.pref.phot_set.gain * USHRT_MAX_DOUBLE;
204 } else {
205 return com.pref.phot_set.gain;
206 }
207 }
208
getInnerRadius()209 static double getInnerRadius() {
210 return com.pref.phot_set.inner;
211 }
212
getOuterRadius()213 static double getOuterRadius() {
214 return com.pref.phot_set.outer;
215 }
216
getMagErr(double intensity,double area,int nsky,double skysig)217 static double getMagErr(double intensity, double area, int nsky, double skysig) {
218 double skyvar, sigsq;
219 double err1, err2, err3;
220 double phpadu;
221
222 skyvar = skysig * skysig; /* variance of the sky brightness */
223 sigsq = skyvar / nsky; /* square of the standard error of the mean sky brightness */
224 phpadu = getCameraGain();
225 err1 = area * skyvar;
226 err2 = intensity / phpadu;
227 err3 = sigsq * area * area;
228
229 return fmin(9.999, 1.0857 * sqrt(err1 + err2 + err3) / intensity);
230 }
231
lo_data()232 static double lo_data() {
233 if (gfit.type == DATA_FLOAT) {
234 return (double) com.pref.phot_set.minval / USHRT_MAX_DOUBLE;
235 } else {
236 return (double) com.pref.phot_set.minval;
237 }
238 }
239
hi_data()240 static double hi_data() {
241 if (gfit.type == DATA_FLOAT) {
242 return (double) com.pref.phot_set.maxval / USHRT_MAX_DOUBLE;
243 } else {
244 return (double) com.pref.phot_set.maxval;
245 }
246 }
247
248 /* Function that compute all photometric data. The result must be freed */
getPhotometryData(gsl_matrix * z,fitted_PSF * psf,gboolean verbose)249 photometry *getPhotometryData(gsl_matrix* z, fitted_PSF *psf, gboolean verbose) {
250 int width = z->size2;
251 int height = z->size1;
252 int n_sky = 0, ret;
253 int x, y, x1, y1, x2, y2;
254 double r1, r2, r, rmin_sq, appRadius;
255 double xc, yc;
256 double apmag = 0.0, mean = 0.0, stdev = 0.0, area = 0.0;
257 double signalIntensity;
258 gboolean valid = TRUE;
259 photometry *phot;
260
261 xc = psf->x0 - 1;
262 yc = psf->y0 - 1;
263
264 if ((xc <= 0.0) || (yc <= 0.0)) return NULL;
265
266 r1 = getInnerRadius();
267 r2 = getOuterRadius();
268 appRadius = psf->fwhmx * 2.0; // in order to be sure to contain star
269 if (appRadius >= r1) {
270 if (verbose) {
271 /* Translator note: radii is plural for radius */
272 siril_log_message(_("Inner and outer radii are too small. Please update values in preferences.\n"));
273 }
274 return NULL;
275 }
276
277 x1 = xc - r2;
278 if (x1 < 1)
279 x1 = 1;
280 x2 = xc + r2;
281 if (x2 > width - 1)
282 x2 = width - 1;
283 y1 = yc - r2;
284 if (y1 < 1)
285 y1 = 1;
286 y2 = yc + r2;
287 if (y2 > height - 1)
288 y2 = height - 1;
289
290 r1 *= r1;
291 r2 *= r2;
292 rmin_sq = (appRadius - 0.5) * (appRadius - 0.5);
293
294 int ndata = (y2 - y1) * (x2 - x1);
295 if (ndata <= 0) {
296 siril_log_color_message(_("An error occurred in your selection. Please make another selection.\n"), "red");
297 return NULL;
298 }
299 double *data = calloc(ndata, sizeof(double));
300 if (!data) {
301 PRINT_ALLOC_ERR;
302 return NULL;
303 }
304
305 for (y = y1; y <= y2; ++y) {
306 int yp = (y - yc) * (y - yc);
307 for (x = x1; x <= x2; ++x) {
308 r = yp + (x - xc) * (x - xc);
309 double pixel = gsl_matrix_get(z, y, x);
310 if (pixel > lo_data() && pixel < hi_data()) {
311 double f = (r < rmin_sq ? 1 : appRadius - sqrt(r) + 0.5);
312 if (f >= 0) {
313 area += f;
314 apmag += pixel * f;
315 }
316 /* annulus */
317 if (r < r2 && r > r1) {
318 data[n_sky] = pixel;
319 n_sky++;
320 }
321 } else {
322 valid = FALSE;
323 }
324 }
325 }
326 if (area < 1) {
327 free(data);
328 return NULL;
329 }
330 if (n_sky < min_sky) {
331 if (verbose) {
332 siril_log_message(_("Warning: There aren't enough pixels"
333 " in the sky annulus. You need to make a larger selection.\n"));
334 }
335 free(data);
336 return NULL;
337 }
338
339 ret = robustmean(n_sky, data, &mean, &stdev);
340 free(data);
341 if (ret > 0) {
342 return NULL;
343 }
344
345 phot = phot_alloc();
346 if (phot) {
347 signalIntensity = apmag - (area * mean);
348
349 phot->mag = getMagnitude(signalIntensity);
350 phot->s_mag = getMagErr(signalIntensity, area, n_sky, stdev);
351 phot->valid = valid;
352 }
353
354 return phot;
355 }
356
initialize_photometric_param()357 void initialize_photometric_param() {
358 com.pref.phot_set.inner = 20;
359 com.pref.phot_set.outer = 30;
360 com.pref.phot_set.gain = 2.3;
361 com.pref.phot_set.minval = 0;
362 com.pref.phot_set.maxval = 60000;
363 }
364