00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031
00032
00033
00034
00035
00036 #include <math.h>
00037 #include <cpl.h>
00038
00039 #include "irplib_spectrum.h"
00040
00041
00042
00043
00044
00045 #define SPECTRUM_HW 16
00046 #define MIN_THRESH_FACT 0.9
00047 #define MAX_THRESH_FACT 1.1
00048 #define SPEC_SHADOW_FACT 30.0
00049 #define SPEC_MAXWIDTH 48
00050
00051
00052
00053
00054
00055 static int select_valid_spectra(cpl_image *, cpl_apertures *, int,
00056 spec_shadows, int, int *, int **) ;
00057 static int valid_spectrum(cpl_image *, cpl_apertures *, int, spec_shadows, int,
00058 int) ;
00059
00060
00064
00065
00068
00083
00084 int irplib_spectrum_find_brightest(
00085 cpl_image * in,
00086 int offset,
00087 spec_shadows shadows,
00088 double min_bright,
00089 int orient,
00090 double * pos)
00091 {
00092 cpl_image * loc_ima ;
00093 cpl_image * filt_image ;
00094 cpl_matrix * kernel ;
00095 cpl_image * collapsed ;
00096 float * pcollapsed ;
00097 cpl_vector * line ;
00098 double * pline ;
00099 cpl_vector * line_filt ;
00100 double threshold ;
00101 double median, stdev, max, mean ;
00102 cpl_mask * mask ;
00103 cpl_image * labels ;
00104 int nlabels ;
00105 cpl_apertures * aperts ;
00106 int n_valid_specs ;
00107 int * valid_specs ;
00108 double brightness, brightest ;
00109 int i ;
00110
00111
00112 if (in == NULL) return -1 ;
00113 if (orient!=0 && orient!=1) return -1 ;
00114
00115
00116 if (orient == 1) {
00117 loc_ima = cpl_image_duplicate(in) ;
00118 cpl_image_flip(loc_ima, 1) ;
00119 } else {
00120 loc_ima = cpl_image_duplicate(in) ;
00121 }
00122
00123
00124 kernel = cpl_matrix_new(3, 3) ;
00125 cpl_matrix_fill(kernel, 1.0) ;
00126 if ((filt_image = cpl_image_filter_median(loc_ima, kernel)) == NULL) {
00127 cpl_matrix_delete(kernel) ;
00128 cpl_image_delete(loc_ima) ;
00129 cpl_msg_error(cpl_func, "cannot filter the image") ;
00130 return -1 ;
00131 }
00132 cpl_image_delete(loc_ima) ;
00133 cpl_matrix_delete(kernel) ;
00134
00135
00136 if ((collapsed = cpl_image_collapse_median_create(filt_image, 1, 0,
00137 0)) == NULL) {
00138 cpl_msg_error(cpl_func, "collapsing image: aborting spectrum detection");
00139 cpl_image_delete(filt_image) ;
00140 return -1 ;
00141 }
00142 cpl_image_delete(filt_image) ;
00143
00144
00145 line = cpl_vector_new_from_image_column(collapsed, 1) ;
00146 cpl_image_delete(collapsed) ;
00147 line_filt = cpl_vector_filter_median_create(line, SPECTRUM_HW) ;
00148 cpl_vector_subtract(line, line_filt) ;
00149 cpl_vector_delete(line_filt) ;
00150
00151
00152 median = cpl_vector_get_median(line) ;
00153 stdev = cpl_vector_get_stdev(line) ;
00154 max = cpl_vector_get_max(line) ;
00155 mean = cpl_vector_get_mean(line) ;
00156
00157
00158 threshold = median + stdev ;
00159 if (threshold > MIN_THRESH_FACT * max) threshold = MIN_THRESH_FACT * max ;
00160 if (threshold < MAX_THRESH_FACT * mean) threshold = MAX_THRESH_FACT * mean;
00161
00162
00163 collapsed = cpl_image_new(1, cpl_vector_get_size(line), CPL_TYPE_FLOAT) ;
00164 pcollapsed = cpl_image_get_data_float(collapsed) ;
00165 pline = cpl_vector_get_data(line) ;
00166 for (i=0 ; i<cpl_vector_get_size(line) ; i++)
00167 pcollapsed[i] = (float)pline[i] ;
00168 cpl_vector_delete(line) ;
00169
00170
00171 if ((mask = cpl_mask_threshold_image_create(collapsed, threshold,
00172 CPL_PIXEL_MAXVAL)) == NULL) {
00173 cpl_msg_error(cpl_func, "cannot binarise") ;
00174 cpl_image_delete(collapsed) ;
00175 return -1 ;
00176 }
00177 if (cpl_mask_count(mask) < 1) {
00178 cpl_msg_error(cpl_func, "not enough signal to detect spectra") ;
00179 cpl_image_delete(collapsed) ;
00180 cpl_mask_delete(mask) ;
00181 return -1 ;
00182 }
00183
00184 if ((labels = cpl_image_labelise_mask_create(mask, &nlabels))==NULL) {
00185 cpl_msg_error(cpl_func, "cannot labelise") ;
00186 cpl_image_delete(collapsed) ;
00187 cpl_mask_delete(mask) ;
00188 return -1 ;
00189 }
00190 cpl_mask_delete(mask) ;
00191
00192
00193 if ((aperts = cpl_apertures_new_from_image(collapsed, labels)) == NULL) {
00194 cpl_msg_error(cpl_func, "cannot compute apertures") ;
00195 cpl_image_delete(collapsed) ;
00196 cpl_image_delete(labels) ;
00197 return -1 ;
00198 }
00199 cpl_image_delete(labels) ;
00200
00201
00202 if (select_valid_spectra(collapsed, aperts, offset, shadows, SPEC_MAXWIDTH,
00203 &n_valid_specs, &valid_specs) == -1) {
00204 cpl_msg_debug(cpl_func, "cannot select valid spectra") ;
00205 cpl_image_delete(collapsed) ;
00206 cpl_apertures_delete(aperts) ;
00207 return -1 ;
00208 }
00209 cpl_image_delete(collapsed) ;
00210 if (n_valid_specs < 1) {
00211 cpl_msg_error(cpl_func, "no valid spectrum detected") ;
00212 cpl_free(valid_specs) ;
00213 cpl_apertures_delete(aperts) ;
00214 return -1 ;
00215 }
00216
00217
00218 *pos = cpl_apertures_get_centroid_y(aperts, valid_specs[0]+1) ;
00219 brightest = valid_specs[0] ;
00220 brightness = cpl_apertures_get_flux(aperts, valid_specs[0]+1) ;
00221 for (i=0 ; i<n_valid_specs ; i++) {
00222 if (cpl_apertures_get_flux(aperts, valid_specs[i]+1) > brightness) {
00223 *pos = cpl_apertures_get_centroid_y(aperts, valid_specs[i]+1) ;
00224 brightest = valid_specs[i] ;
00225 brightness = cpl_apertures_get_flux(aperts, valid_specs[i]+1) ;
00226 }
00227 }
00228 cpl_apertures_delete(aperts) ;
00229 cpl_free(valid_specs) ;
00230
00231
00232 if (brightness < min_bright) {
00233 cpl_msg_error(cpl_func, "brightness %f too low <%f", brightness,
00234 min_bright) ;
00235 return -1 ;
00236 }
00237
00238
00239 return 0 ;
00240 }
00241
00244
00256
00257 static int select_valid_spectra(
00258 cpl_image * in,
00259 cpl_apertures * aperts,
00260 int offset,
00261 spec_shadows shadows,
00262 int max_spec_width,
00263 int * n_valid_specs,
00264 int ** valid_specs)
00265 {
00266 int nb_aperts ;
00267 int i, j ;
00268
00269
00270 *valid_specs = NULL ;
00271 nb_aperts = cpl_apertures_get_size(aperts) ;
00272 *n_valid_specs = 0 ;
00273
00274
00275 if (nb_aperts < 1) return -1 ;
00276
00277
00278 j = 0 ;
00279 for (i=0 ; i<nb_aperts ; i++)
00280 if (valid_spectrum(in, aperts, offset, shadows, max_spec_width,
00281 i+1)) (*n_valid_specs)++ ;
00282
00283
00284 if (*n_valid_specs) {
00285 *valid_specs = cpl_calloc(*n_valid_specs, sizeof(int)) ;
00286 j = 0 ;
00287 for (i=0 ; i<nb_aperts ; i++)
00288 if (valid_spectrum(in, aperts, offset, shadows, max_spec_width,
00289 i+1)) {
00290 (*valid_specs)[j] = i ;
00291 j++ ;
00292 }
00293 } else return -1 ;
00294
00295 return 0 ;
00296 }
00297
00298
00309
00310 static int valid_spectrum(
00311 cpl_image * in,
00312 cpl_apertures * aperts,
00313 int offset,
00314 spec_shadows shadows,
00315 int max_spec_width,
00316 int objnum)
00317 {
00318 int objwidth ;
00319 double valover, valunder, valcenter ;
00320
00321
00322 objwidth = cpl_apertures_get_top(aperts, objnum) -
00323 cpl_apertures_get_bottom(aperts, objnum) + 1 ;
00324 if (objwidth > max_spec_width) {
00325 cpl_msg_error(cpl_func, "object is too wide") ;
00326 return 0 ;
00327 }
00328
00329
00330 if (cpl_apertures_get_npix(aperts, objnum) < 2) return 0 ;
00331
00332
00333 if (shadows == NO_SHADOW) return 1 ;
00334
00335
00336 valcenter = cpl_apertures_get_median(aperts, objnum) ;
00337
00338
00339 if (cpl_apertures_get_bottom(aperts, objnum) - offset < 1) valunder = 0.0 ;
00340 else valunder = cpl_image_get_median_window(in, 1,
00341 cpl_apertures_get_bottom(aperts, objnum) - offset, 1,
00342 cpl_apertures_get_top(aperts, objnum) - offset) ;
00343
00344 if (cpl_apertures_get_top(aperts, objnum) + offset > 1024) valover = 0.0 ;
00345 else valover = cpl_image_get_median_window(in, 1,
00346 cpl_apertures_get_bottom(aperts, objnum) + offset, 1,
00347 cpl_apertures_get_top(aperts, objnum) + offset) ;
00348
00349 switch (shadows) {
00350 case TWO_SHADOWS:
00351 if ((valunder < -fabs(valcenter/SPEC_SHADOW_FACT)) &&
00352 (valover < -fabs(valcenter/SPEC_SHADOW_FACT)) &&
00353 (valunder/valover > 0.5) &&
00354 (valunder/valover < 2.0)) return 1 ;
00355 else return 0 ;
00356
00357 case ONE_SHADOW:
00358 if ((valunder < -fabs(valcenter/SPEC_SHADOW_FACT)) ||
00359 (valover < -fabs(valcenter/SPEC_SHADOW_FACT))) return 1 ;
00360 else return 0 ;
00361
00362 case NO_SHADOW:
00363 return 1 ;
00364
00365 default:
00366 cpl_msg_error(cpl_func, "unknown spec_detect_mode") ;
00367 break ;
00368 }
00369
00370 return 0 ;
00371 }