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 <stdio.h>
22 #include <assert.h>
23 #include <stdlib.h>
24 #include <string.h>
25 #include <math.h>
26 #include <locale.h>
27 #include <gsl/gsl_matrix.h>
28
29 #include "core/siril.h"
30 #include "core/proto.h"
31 #include "algos/Def_Wavelet.h"
32 #include "gui/utils.h"
33 #include "gui/message_dialog.h"
34 #include "gui/image_display.h"
35 #include "gui/progress_and_log.h"
36 #include "gui/PSF_list.h"
37 #include "algos/PSF.h"
38 #include "algos/star_finder.h"
39 #include "algos/statistics.h"
40 #include "filters/wavelets.h"
41 #include "io/image_format_fits.h"
42 #include "core/OS_utils.h"
43
44 #define WAVELET_SCALE 3
45
guess_resolution(fits * fit)46 static double guess_resolution(fits *fit) {
47 double focal = fit->focal_length;
48 double size = fit->pixel_size_x;
49 double bin;
50
51 if ((focal <= 0.0) || (size <= 0.0)) {
52 focal = com.pref.focal;
53 size = com.pref.pitch;
54 /* if we have no way to guess, we return the
55 * flag -1
56 */
57 if ((focal <= 0.0) || (size <= 0.0))
58 return -1.0;
59 }
60 bin = ((fit->binning_x + fit->binning_y) / 2.0);
61 if (bin <= 0) bin = 1.0;
62
63 double res = RADCONV / focal * size * bin;
64
65 /* test for high value. In this case we increase
66 * the number of detected star in is_star function
67 */
68 if (res > 20.0) return -1.0;
69 /* if res > 1.0 we use default radius value */
70 if (res < 0.1 || res > 1.0) return 1.0;
71 return res;
72 }
73
compute_threshold(fits * fit,double ksigma,int layer,float * norm,double * bg)74 static float compute_threshold(fits *fit, double ksigma, int layer, float *norm, double *bg) {
75 float threshold;
76 imstats *stat;
77
78 assert(layer <= 3);
79
80 stat = statistics(NULL, -1, fit, layer, NULL, STATS_BASIC, FALSE);
81 if (!stat) {
82 siril_log_message(_("Error: statistics computation failed.\n"));
83 *norm = 0;
84 *bg = 0.0;
85 return 0;
86 }
87 threshold = (float)(stat->median + ksigma * stat->bgnoise);
88 *norm = stat->normValue;
89 *bg = stat->median;
90 free_stats(stat);
91
92 return threshold;
93 }
94
is_star(fitted_PSF * result,star_finder_params * sf)95 static gboolean is_star(fitted_PSF *result, star_finder_params *sf ) {
96 /* here this is a bit trick, bu if no resolution is computed
97 * we take a greater factor for star size. This is less optimize
98 * but at least it reproduces almost the old behavior
99 */
100 int factor = sf->no_guess ? 10 : 3;
101 if (isnan(result->fwhmx) || isnan(result->fwhmy))
102 return FALSE;
103 if (isnan(result->x0) || isnan(result->y0))
104 return FALSE;
105 if (isnan(result->mag))
106 return FALSE;
107 if ((result->x0 <= 0.0) || (result->y0 <= 0.0))
108 return FALSE;
109 if (result->sx > factor * sf->adj_radius || result->sy > factor * sf->adj_radius)
110 return FALSE;
111 if (result->fwhmx <= 0.0 || result->fwhmy <= 0.0)
112 return FALSE;
113 if ((result->fwhmy / result->fwhmx) < sf->roundness)
114 return FALSE;
115 return TRUE;
116 }
117
get_structure(star_finder_params * sf)118 static void get_structure(star_finder_params *sf) {
119 static GtkSpinButton *spin_radius = NULL, *spin_sigma = NULL,
120 *spin_roundness = NULL;
121 static GtkToggleButton *toggle_adjust = NULL;
122
123 if (spin_radius == NULL) {
124 spin_radius = GTK_SPIN_BUTTON(lookup_widget("spinstarfinder_radius"));
125 spin_sigma = GTK_SPIN_BUTTON(lookup_widget("spinstarfinder_threshold"));
126 spin_roundness = GTK_SPIN_BUTTON(lookup_widget("spinstarfinder_round"));
127 toggle_adjust = GTK_TOGGLE_BUTTON(lookup_widget("toggle_radius_adjust"));
128 }
129 sf->radius = (int) gtk_spin_button_get_value(spin_radius);
130 sf->sigma = gtk_spin_button_get_value(spin_sigma);
131 sf->adjust = gtk_toggle_button_get_active(toggle_adjust);
132 sf->roundness = gtk_spin_button_get_value(spin_roundness);
133 }
134
init_peaker_GUI()135 void init_peaker_GUI() {
136 /* TODO someday: read values from conf file and set them in the GUI.
137 * Until then, storing values in com.starfinder_conf instead of getting
138 * them in the GUI while running the peaker.
139 * see also init_peaker_default below */
140 get_structure(&com.starfinder_conf);
141 }
142
init_peaker_default()143 void init_peaker_default() {
144 /* values taken from siril3.glade */
145 com.starfinder_conf.radius = 10;
146 com.starfinder_conf.adjust = TRUE;
147 com.starfinder_conf.sigma = 1.0;
148 com.starfinder_conf.roundness = 0.5;
149 com.starfinder_conf.no_guess = FALSE;
150 }
151
on_toggle_radius_adjust_toggled(GtkToggleButton * togglebutton,gpointer user_data)152 void on_toggle_radius_adjust_toggled(GtkToggleButton *togglebutton, gpointer user_data) {
153 com.starfinder_conf.adjust = gtk_toggle_button_get_active(togglebutton);
154 }
155
on_spin_sf_radius_changed(GtkSpinButton * spinbutton,gpointer user_data)156 void on_spin_sf_radius_changed(GtkSpinButton *spinbutton, gpointer user_data) {
157 com.starfinder_conf.radius = (int)gtk_spin_button_get_value(spinbutton);
158 }
159
on_spin_sf_threshold_changed(GtkSpinButton * spinbutton,gpointer user_data)160 void on_spin_sf_threshold_changed(GtkSpinButton *spinbutton, gpointer user_data) {
161 com.starfinder_conf.sigma = gtk_spin_button_get_value(spinbutton);
162 }
163
on_spin_sf_roundness_changed(GtkSpinButton * spinbutton,gpointer user_data)164 void on_spin_sf_roundness_changed(GtkSpinButton *spinbutton, gpointer user_data) {
165 com.starfinder_conf.roundness = gtk_spin_button_get_value(spinbutton);
166 }
167
update_peaker_GUI()168 void update_peaker_GUI() {
169 static GtkSpinButton *spin_radius = NULL, *spin_sigma = NULL,
170 *spin_roundness = NULL;
171 static GtkToggleButton *toggle_adjust = NULL;
172
173 if (spin_radius == NULL) {
174 spin_radius = GTK_SPIN_BUTTON(lookup_widget("spinstarfinder_radius"));
175 spin_sigma = GTK_SPIN_BUTTON(lookup_widget("spinstarfinder_threshold"));
176 spin_roundness = GTK_SPIN_BUTTON(lookup_widget("spinstarfinder_round"));
177 toggle_adjust = GTK_TOGGLE_BUTTON(lookup_widget("toggle_radius_adjust"));
178 }
179 gtk_spin_button_set_value(spin_radius, (double) com.starfinder_conf.radius);
180 gtk_toggle_button_set_active(toggle_adjust, com.starfinder_conf.adjust);
181 gtk_spin_button_set_value(spin_sigma, com.starfinder_conf.sigma);
182 gtk_spin_button_set_value(spin_roundness, com.starfinder_conf.roundness);
183 }
184
confirm_peaker_GUI()185 void confirm_peaker_GUI() {
186 static GtkSpinButton *spin_radius = NULL, *spin_sigma = NULL,
187 *spin_roundness = NULL;
188
189 if (spin_radius == NULL) {
190 spin_radius = GTK_SPIN_BUTTON(lookup_widget("spinstarfinder_radius"));
191 spin_sigma = GTK_SPIN_BUTTON(lookup_widget("spinstarfinder_threshold"));
192 spin_roundness = GTK_SPIN_BUTTON(lookup_widget("spinstarfinder_round"));
193 }
194 gtk_spin_button_update(spin_radius);
195 gtk_spin_button_update(spin_sigma);
196 gtk_spin_button_update(spin_roundness);
197 }
198
199
200 /*
201 This is an implementation of a simple peak detector algorithm which
202 identifies any pixel that is greater than any of its eight neighbors.
203
204 Original algorithm come from:
205 Copyleft (L) 1998 Kenneth J. Mighell (Kitt Peak National Observatory)
206 */
207 static int minimize_candidates(fits *image, star_finder_params *sf, double bg, pointi *candidates, int nb_candidates, int layer, fitted_PSF ***retval);
208
peaker(fits * fit,int layer,star_finder_params * sf,int * nb_stars,rectangle * area,gboolean showtime)209 fitted_PSF **peaker(fits *fit, int layer, star_finder_params *sf, int *nb_stars, rectangle *area, gboolean showtime) {
210 int nx = fit->rx;
211 int ny = fit->ry;
212 int areaX0 = 0;
213 int areaY0 = 0;
214 int areaX1 = nx;
215 int areaY1 = ny;
216 int nbstars = 0;
217 double bg;
218 float threshold, norm;
219 float **wave_image;
220 fits wave_fit = { 0 };
221 pointi *candidates;
222
223 struct timeval t_start, t_end;
224
225 assert(nx > 0 && ny > 0);
226
227 siril_log_color_message(_("Findstar: processing...\n"), "green");
228 gettimeofday(&t_start, NULL);
229
230 /* running statistics on the input image is best as it caches them */
231 threshold = compute_threshold(fit, sf->sigma * 5.0, layer, &norm, &bg);
232 if (norm == 0.0f)
233 return NULL;
234
235 siril_debug_print("Threshold: %f (background: %f, norm: %f)\n", threshold, bg, norm);
236
237 /* if fit is ushort, we need to get it in float [0, 65535] to run
238 * statistics only once: stats give the threshold used in the
239 * wavelets-filtered image for candidate detection and the background
240 * used in PSF fitting from real image data which can be still ushort.
241 */
242 if (extract_fits(fit, &wave_fit, layer, TRUE)) {
243 siril_log_message(_("Failed to copy the image for processing\n"));
244 return NULL;
245 }
246 get_wavelet_layers(&wave_fit, WAVELET_SCALE, 2, TO_PAVE_BSPLINE, layer);
247
248 /* Build 2D representation of wavelet image upside-down */
249 wave_image = malloc(ny * sizeof(float *));
250 if (!wave_image) {
251 PRINT_ALLOC_ERR;
252 clearfits(&wave_fit);
253 return NULL;
254 }
255 for (int k = 0; k < ny; k++) {
256 wave_image[ny - k - 1] = wave_fit.fdata + k * nx;
257 }
258
259 if (area) {
260 areaX0 = area->x;
261 areaY0 = area->y;
262 areaX1 = area->w + areaX0;
263 areaY1 = area->h + areaY0;
264 }
265
266 candidates = malloc(MAX_STARS * sizeof(pointi));
267 if (!candidates) {
268 PRINT_ALLOC_ERR;
269 return NULL;
270 }
271
272 double res = guess_resolution(fit);
273
274 if (res < 0) {
275 sf->no_guess = TRUE;
276 res = 1.0;
277 }
278
279 sf->adj_radius = sf->adjust ? sf->radius / res : sf->radius;
280
281 /* Search for candidate stars in the filtered image */
282 for (int y = sf->adj_radius + areaY0; y < areaY1 - sf->adj_radius; y++) {
283 for (int x = sf->adj_radius + areaX0; x < areaX1 - sf->adj_radius; x++) {
284 float pixel = wave_image[y][x];
285 if (pixel > threshold && pixel < norm) {
286 gboolean bingo = TRUE;
287 float neighbor;
288 for (int yy = y - 1; yy <= y + 1; yy++) {
289 for (int xx = x - 1; xx <= x + 1; xx++) {
290 if (xx == x && yy == y)
291 continue;
292 neighbor = wave_image[yy][xx];
293 if (neighbor > pixel) {
294 bingo = FALSE;
295 break;
296 } else if (neighbor == pixel) {
297 if ((xx <= x && yy <= y) || (xx > x && yy < y)) {
298 bingo = FALSE;
299 break;
300 }
301 }
302 }
303 }
304 if (bingo && nbstars < MAX_STARS) {
305 candidates[nbstars].x = x;
306 candidates[nbstars].y = y;
307 nbstars++;
308 if (nbstars == MAX_STARS) break;
309 }
310 }
311 }
312 if (nbstars == MAX_STARS) break;
313 }
314 clearfits(&wave_fit);
315 siril_debug_print("Candidates for stars: %d\n", nbstars);
316
317 /* Check if candidates are stars by minimizing a PSF on each */
318 fitted_PSF **results;
319 nbstars = minimize_candidates(fit, sf, bg, candidates, nbstars, layer, &results);
320 if (nbstars == 0)
321 results = NULL;
322 sort_stars(results, nbstars);
323 free(wave_image);
324 free(candidates);
325
326 gettimeofday(&t_end, NULL);
327 if (showtime)
328 show_time(t_start, t_end);
329 if (nb_stars)
330 *nb_stars = nbstars;
331 return results;
332 }
333
334 /* returns number of stars found, result is in parameters */
minimize_candidates(fits * image,star_finder_params * sf,double bg,pointi * candidates,int nb_candidates,int layer,fitted_PSF *** retval)335 static int minimize_candidates(fits *image, star_finder_params *sf, double bg, pointi *candidates, int nb_candidates, int layer, fitted_PSF ***retval) {
336 int radius = sf->adj_radius;
337 int nx = image->rx;
338 int ny = image->ry;
339 gsl_matrix *z = gsl_matrix_alloc(radius * 2, radius * 2);
340 WORD **image_ushort = NULL;
341 float **image_float = NULL;
342 gint nbstars = 0;
343
344 if (image->type == DATA_USHORT) {
345 image_ushort = malloc(ny * sizeof(WORD *));
346 for (int k = 0; k < ny; k++)
347 image_ushort[ny - k - 1] = image->pdata[layer] + k * nx;
348 }
349 else if (image->type == DATA_FLOAT) {
350 image_float = malloc(ny * sizeof(float *));
351 for (int k = 0; k < ny; k++)
352 image_float[ny - k - 1] = image->fpdata[layer] + k * nx;
353 }
354 else return 0;
355
356 fitted_PSF **results = malloc((nb_candidates + 1) * sizeof(fitted_PSF *));
357 if (!results) {
358 PRINT_ALLOC_ERR;
359 return 0;
360 }
361
362 for (int candidate = 0; candidate < nb_candidates; candidate++) {
363 int x = candidates[candidate].x, y = candidates[candidate].y;
364 int ii, jj, i, j;
365 /* FILL z */
366 if (image->type == DATA_USHORT) {
367 for (jj = 0, j = y - radius; j < y + radius; j++, jj++) {
368 for (ii = 0, i = x - radius; i < x + radius;
369 i++, ii++) {
370 gsl_matrix_set(z, ii, jj, (double)image_ushort[j][i]);
371 }
372 }
373 } else {
374 for (jj = 0, j = y - radius; j < y + radius; j++, jj++) {
375 for (ii = 0, i = x - radius; i < x + radius;
376 i++, ii++) {
377 gsl_matrix_set(z, ii, jj, (double)image_float[j][i]);
378 }
379 }
380 }
381
382 fitted_PSF *cur_star = psf_global_minimisation(z, bg, FALSE, FALSE, FALSE);
383 if (cur_star) {
384 if (is_star(cur_star, sf)) {
385 //fwhm_to_arcsec_if_needed(image, cur_star); // should we do this here?
386 cur_star->layer = layer;
387 int result_index = g_atomic_int_add(&nbstars, 1);
388 cur_star->xpos = (x - radius) + cur_star->x0 - 1.0;
389 cur_star->ypos = (y - radius) + cur_star->y0 - 1.0;
390 results[result_index] = cur_star;
391 //fprintf(stdout, "%03d: %11f %11f %f\n",
392 // result_index, cur_star->xpos, cur_star->ypos, cur_star->mag);
393 }
394 else free(cur_star);
395 }
396 }
397
398 results[nbstars] = NULL;
399 if (retval)
400 *retval = results;
401 gsl_matrix_free(z);
402 if (image_ushort) free(image_ushort);
403 if (image_float) free(image_float);
404 return nbstars;
405 }
406
407 /* Function to add star one by one, from the selection rectangle, the
408 * minimization is run and the star is detected and added to the list of stars.
409 *
410 * IF A STAR IS FOUND and not already present in com.stars, the return value is
411 * the new star and index is set to the index of the new star in com.stars.
412 * IF NO NEW STAR WAS FOUND, either because it was already in the list, or a
413 * star failed to be detected in the selection, or any other error, the return
414 * value is NULL and index is set to -1.
415 */
add_star(fits * fit,int layer,int * index)416 fitted_PSF *add_star(fits *fit, int layer, int *index) {
417 int i = 0;
418 gboolean already_found = FALSE;
419
420 *index = -1;
421 fitted_PSF * result = psf_get_minimisation(&gfit, layer, &com.selection, FALSE, TRUE, TRUE);
422 if (!result)
423 return NULL;
424 /* We do not check if it's matching with the "is_star()" criteria.
425 * Indeed, in this case the user can add manually stars missed by star_finder */
426
427 if (com.stars) {
428 // check if the star was already detected/peaked
429 while (com.stars[i]) {
430 if (fabs(result->x0 + com.selection.x - com.stars[i]->xpos) < 0.9
431 && fabs(
432 com.selection.y + com.selection.h - result->y0
433 - com.stars[i]->ypos) < 0.9)
434 already_found = TRUE;
435 i++;
436 }
437 } else {
438 com.stars = malloc((MAX_STARS + 1) * sizeof(fitted_PSF*));
439 if (!com.stars) {
440 PRINT_ALLOC_ERR;
441 return NULL;
442 }
443 com.star_is_seqdata = FALSE;
444 }
445
446 if (already_found) {
447 free(result);
448 result = NULL;
449 char *msg = siril_log_message(_("This star has already been picked !\n"));
450 siril_message_dialog( GTK_MESSAGE_INFO, _("Peaker"), msg);
451 } else {
452 if (i < MAX_STARS) {
453 result->xpos = result->x0 + com.selection.x;
454 result->ypos = com.selection.y + com.selection.h - result->y0;
455 com.stars[i] = result;
456 com.stars[i + 1] = NULL;
457 *index = i;
458 } else {
459 free(result);
460 result = NULL;
461 }
462 }
463 return result;
464 }
465
get_size_star_tab()466 int get_size_star_tab() {
467 int i = 0;
468 while (com.stars[i])
469 i++;
470 return i;
471 }
472
473 /* Remove a star from com.stars, at index index. The star is freed. */
remove_star(int index)474 int remove_star(int index) {
475 if (index < 0 || !com.stars || !com.stars[index])
476 return 1;
477
478 int N = get_size_star_tab() + 1;
479
480 free(com.stars[index]);
481 memmove(&com.stars[index], &com.stars[index + 1],
482 (N - index - 1) * sizeof(*com.stars));
483 redraw(com.cvport, REMAP_NONE);
484 return 0;
485 }
486
compare_stars(const void * star1,const void * star2)487 int compare_stars(const void* star1, const void* star2) {
488 fitted_PSF *s1 = *(fitted_PSF**) star1;
489 fitted_PSF *s2 = *(fitted_PSF**) star2;
490
491 if (s1->mag < s2->mag)
492 return -1;
493 else if (s1->mag > s2->mag)
494 return 1;
495 else
496 return 0;
497 }
498
sort_stars(fitted_PSF ** stars,int total)499 void sort_stars(fitted_PSF **stars, int total) {
500 if (*(&stars))
501 qsort(*(&stars), total, sizeof(fitted_PSF*), compare_stars);
502 }
503
free_fitted_stars(fitted_PSF ** stars)504 void free_fitted_stars(fitted_PSF **stars) {
505 int i = 0;
506 while (stars && stars[i])
507 free(stars[i++]);
508 free(stars);
509 }
510
FWHM_average(fitted_PSF ** stars,int nb,float * FWHMx,float * FWHMy,char ** units)511 void FWHM_average(fitted_PSF **stars, int nb, float *FWHMx, float *FWHMy, char **units) {
512 if (stars && stars[0]) {
513 int i = 0;
514
515 *FWHMx = 0.0f;
516 *FWHMy = 0.0f;
517 *units = stars[0]->units;
518 while (i < nb) {
519 *FWHMx += (float) stars[i]->fwhmx;
520 *FWHMy += (float) stars[i]->fwhmy;
521 i++;
522 }
523 *FWHMx = *FWHMx / (float)i;
524 *FWHMy = *FWHMy / (float)i;
525 }
526 }
527