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
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
|
/*
* This file is part of the ESO SINFONI Pipeline
* Copyright (C) 2004,2005 European Southern Observatory
*
* This program 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 2 of the License, or
* (at your option) any later version.
*
* This program 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 this program; if not, write to the Free Software
* Foundation, 51 Franklin St, Fifth Floor, Boston, MA 02111-1307 USA
*/
/*----------------------------------------------------------------------------
File name : sinfo_bp_norm.c
Author : J. Schreiber
Created on : May 5, 2003
Description : Different methods for searching for bad pixels
used in the recipe sinfo_rec_mflat
---------------------------------------------------------------------------*/
#ifdef HAVE_CONFIG_H
# include <config.h>
#endif
/*----------------------------------------------------------------------------
Includes
---------------------------------------------------------------------------*/
#include "sinfo_bp_norm.h"
#include "sinfo_image_ops.h"
#include "sinfo_detlin.h"
#include "sinfo_badnorm_ini_by_cpl.h"
#include "sinfo_baddist_ini_by_cpl.h"
#include "sinfo_pro_save.h"
#include "sinfo_functions.h"
#include "sinfo_pro_types.h"
#include "sinfo_hidden.h"
#include "sinfo_error.h"
#include "sinfo_utils_wrappers.h"
/**@{*/
/**
* @addtogroup sinfo_bad_pix_search Bad Pixel Search
*
* TBD
*/
/*----------------------------------------------------------------------------
Defines
---------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------
Function Definitions
---------------------------------------------------------------------------*/
/**
@name sinfo_new_bp_search_norm
@memo finds pixels whose intensity differs from the median more
than a given threshold
@param plugin_id recipe name
@param config input parameters list
@param sof input set of frames
@param procatg PRO.CATG of product.
@return integer (0 if it worked, -1 if it doesn't)
@doc
this function searches for static bad pixels in stacks of flatfield frames,
that means it stacks the flatfield in a data cube and takes the sinfo_median
along the z-axis to be sure to have no cosmics, then the intensity tilt of
each column is removed.
The standard deviation and mean of the pixel values within a defined
rectangular zone is determined in a way that the extreme low and high
values are cut off.
The next step is to indicate each pixel as bad that has a deviation from
the mean greater than a user defined factor times the standard deviation.
This step can be leaped over if wished.
The next step is to calculate the median of the 8 nearest neighbors for
each pixel and to determine the deviation of each pixel value from this
sinfo_median.
If this deviation is greater than a defined factor times the standard
deviation the pixel value is replaced by the sinfo_median. The last step is
repeated to be able to consider also clusters of bad pixels. Then the
resulting image is compared with the input image with the column intensity
tilt removed and each changed pixel is indicated as bad. Finally, a bad
pixel mask is produced that means each good pixel is marked with 1 and each
bad pixel with 0.
*/
int
sinfo_new_bp_search_normal(const char* plugin_id, cpl_parameterlist* config,
cpl_frameset* sof, cpl_frameset* ref_set,
const char* procatg)
{
bad_config * cfg = NULL;
cpl_imagelist * image_list = NULL;
cpl_image ** med = NULL;
cpl_image * medImage = NULL;
cpl_image * medIm = NULL;
cpl_image * colImage = NULL;
cpl_image * compImage = NULL;
cpl_image * maskImage = NULL;
cpl_image * threshIm = NULL;
Stats * stats = NULL;
cpl_parameter *p = NULL;
int no = 0;
float lo_cut = 0.;
float hi_cut = 0.;
int i = 0;
int n = 0;
int half_box_size = 0;
cpl_frameset* raw = NULL;
cpl_table* qclog_tbl = NULL;
char key_value[FILE_NAME_SZ];
/* parse the file names and parameters to the bad_config data structure
cfg */
sinfo_check_rec_status(0);
check_nomsg(raw = cpl_frameset_new());
sinfo_check_rec_status(1);
if (strcmp(procatg, PRO_BP_MAP_NO) == 0) {
cknull(cfg = sinfo_parse_cpl_input_badnorm(config, sof, procatg, &raw),
"could not parse cpl input!");
}
else if (strcmp(procatg, PRO_BP_MAP_DI) == 0) {
cknull(cfg = sinfo_parse_cpl_input_baddist(config, sof, procatg, &raw),
"could not parse cpl input!");
}
else if (strcmp(procatg, PRO_DEFAULT) == 0) {
cknull(cfg = sinfo_parse_cpl_input_badnorm(config, sof, procatg, &raw),
"could not parse cpl input!");
}
else {
sinfo_msg_error("Error: PRO.CATG %s, not supported!", procatg);
goto cleanup;
}
sinfo_check_rec_status(2);
/* take a clean mean of the frames */
sinfo_msg("Takes a clean mean of the frames");
check_nomsg(image_list = cpl_imagelist_new());
sinfo_check_rec_status(3);
for (i = 0; i < cfg->nframes; i++) {
if (sinfo_is_fits_file(cfg->framelist[i]) != 1) {
sinfo_msg_error("Input file %s is not FITS", cfg->framelist[i]);
goto cleanup;
}
check_nomsg(
cpl_imagelist_set(image_list,
cpl_image_load(cfg->framelist[i],
CPL_TYPE_FLOAT, 0, 0),
i));
}
/* finally take the average image of the cube by
rejecting the extreme values */
check_nomsg(no = cpl_imagelist_get_size(image_list));
lo_cut = (floor)(cfg->loReject * no + 0.5);
hi_cut = (floor)(cfg->hiReject * no + 0.5);
cknull(medImage=cpl_imagelist_collapse_minmax_create(image_list,
lo_cut,
hi_cut),
"error in sinfo_average_with_rejection");
/* free memory */
sinfo_free_imagelist(&image_list);
/*----------------------------------------------
* remove the intensity tilt from every column
* and compute the standard deviation on a rectangular zone
*/
cknull(medIm = sinfo_new_thresh_image(medImage, cfg->mincut, cfg->maxcut),
"error sinfo_new_thresh_image");
cknull(colImage = sinfo_new_col_tilt( medIm, cfg->sigmaFactor ),
"sinfo_colTilt failed" );
cknull(stats = sinfo_new_image_stats_on_rectangle(colImage,
cfg->loReject,
cfg->hiReject,
cfg->llx,
cfg->lly,
cfg->urx,
cfg->ury),
" sinfo_get_image_stats_on_vig failed");
if (stats != NULL ) {
sinfo_msg("Clean stdev: %f\n", stats->cleanstdev );
sinfo_msg("Clean mean: %f\n", stats->cleanmean );
}
/*
cknull_nomsg(qclog_tbl = sinfo_qclog_init());
ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC ICTILT STDEV",
stats->cleanstdev,
"Intensity column clean stdev","%g"));
ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC ICTILT MEAN",
stats->cleanmean,
"Intensity column clean mean","%g"));
ck0(sinfo_pro_save_ima(colImage,raw,sof,
(char*) BP_NORM_INT_COL_TILT_CORR_OUT_FILENAME,
PRO_INT_COL_TILT_COR,qclog_tbl,plugin_id,config),
"cannot save ima %s",
BP_NORM_INT_COL_TILT_CORR_OUT_FILENAME);
sinfo_free_table(&qclog_tbl);
*/
/* indicate pixels with great deviations from the clean mean as bad */
if (cfg->threshInd == 1) {
cknull(threshIm = sinfo_new_thresh_image(colImage,
stats->cleanmean-cfg->meanfactor*stats->cleanstdev,
stats->cleanmean+cfg->meanfactor*stats->cleanstdev),
" sinfo_threshImage failed" );
}
if (cfg->threshInd == 0) {
threshIm = colImage;
}
/* AMO here invalid fread? */
med = (cpl_image**) cpl_calloc(cfg->iterations, sizeof(cpl_image*));
/* filter iteratively the images by a sinfo_median filter of the nearest
neighbors under the condition of a deviation greater than a factor
times the standard deviation */
sinfo_msg("Apply sinfo_median filter on pixel nearest neighbors");
if (cfg->methodInd == 1) {
if (cfg->factor > 0) {
cknull(med[0]=sinfo_new_median_image(threshIm,
-cfg->factor*stats->cleanstdev),
" sinfo_medianImage failed (1)" );
for (i = 0; i < cfg->iterations - 1; i++) {
cknull(med[i + 1] = sinfo_new_median_image(med[i],
-cfg->factor* stats->cleanstdev),
"sinfo_medianImage failed (2)");
}
}
else {
cknull(med[0] = sinfo_new_median_image(threshIm, -cfg->factor),
" sinfo_medianImage failed (1)" );
for (i = 0; i < cfg->iterations - 1; i++) {
cknull(med[i+1] = sinfo_new_median_image(med[i], -cfg->factor),
" sinfo_medianImage failed (%d)",i );
}
}
}
else if (cfg->methodInd == 2) {
cknull(med[0] = sinfo_new_abs_dist_image(threshIm, -cfg->factor),
" sinfo_absDistImage failed (1)" );
for (i = 0; i < cfg->iterations - 1; i++) {
cknull(med[i+1] = sinfo_new_abs_dist_image(med[i], -cfg->factor),
" sinfo_absDistImage failed (2)" );
}
}
else if (cfg->methodInd == 3) {
cknull(
med[0] = sinfo_new_mean_image_in_spec(threshIm,
-cfg->factor * stats->cleanstdev),
"sinfo_meanImageInSpec failed (1)");
for (i = 0; i < cfg->iterations - 1; i++) {
cknull(
med[i + 1] = sinfo_new_mean_image_in_spec(med[i],
-cfg->factor * stats->cleanstdev),
" sinfo_meanImageInSpec failed (2)");
}
}
else if (cfg->methodInd == 4) {
half_box_size = (cfg->urx - cfg->llx) / 2;
cknull(med[0] = sinfo_new_local_median_image(threshIm,
-cfg->factor,
cfg->loReject,
cfg->hiReject,
half_box_size),
" sinfo_localMedianImage failed (1)" );
for (i = 0; i < cfg->iterations - 1; i++) {
cknull(med[i+1] = sinfo_new_local_median_image(med[i],
-cfg->factor,
cfg->loReject,
cfg->hiReject,
half_box_size),
" sinfo_localMedianImage failed (2)" );
}
}
else {
sinfo_msg_error (" wrong indicator methodInd !" );
goto cleanup;
}
/* compare the filtered image with the input image */
cknull(compImage = sinfo_new_compare_images(threshIm,
med[cfg->iterations - 1],
medImage),
" sinfo_compareImages failed" );
/* generate the bad pixel mask */
sinfo_msg("Generates bad pixel map");
cknull(maskImage = sinfo_new_promote_image_to_mask( compImage, &n ),
" error in sinfo_promoteImageToMask" );
sinfo_msg("No of bad pixels: %d\n", n );
cknull_nomsg(qclog_tbl = sinfo_qclog_init());
check_nomsg(p = cpl_parameterlist_find(config, "sinfoni.bp.method"));
snprintf(key_value, MAX_NAME_SIZE - 1, "%s", cpl_parameter_get_string(p));
ck0_nomsg(
sinfo_qclog_add_string(qclog_tbl, "QC BP-MAP METHOD",
key_value, "BP search method", "%s"));
ck0_nomsg(
sinfo_qclog_add_int(qclog_tbl, "QC BP-MAP NBADPIX", n,
"No of bad pixels", "%d"));
ck0(
sinfo_pro_save_ima(maskImage, ref_set, sof, cfg->outName,
procatg, qclog_tbl, plugin_id, config),
"cannot save ima %s", cfg->outName);
sinfo_free_table(&qclog_tbl);
sinfo_free_image(&maskImage);
sinfo_free_image(&compImage);
sinfo_free_image_array(&med, cfg->iterations);
/*
if (med != NULL) {
for ( i = 0 ; i < cfg->iterations ; i++ ) {
if(med[i] != NULL) {
cpl_image_delete(med[i]) ;
med[i]=NULL;
}
}
cpl_free(med) ;
med=NULL;
}
*/
if (stats != NULL ) {
cpl_free(stats);
stats = NULL;
}
sinfo_free_image(&medIm);
sinfo_free_image(&medImage);
sinfo_free_image(&colImage);
if (cfg->threshInd == 1) {
sinfo_free_image(&threshIm);
}
sinfo_badnorm_free(&cfg);
sinfo_free_frameset(&raw);
return 0;
cleanup:
sinfo_free_image_array(&med, cfg->iterations);
/*
if (med != NULL) {
for ( i = 0 ; i < cfg->iterations ; i++ ) {
if(med[i] != NULL) {
cpl_image_delete(med[i]) ;
med[i]=NULL;
}
}
cpl_free(med) ;
med=NULL;
}
*/
sinfo_free_image(&compImage);
sinfo_free_image(&maskImage);
sinfo_free_image(&threshIm);
sinfo_free_table(&qclog_tbl);
sinfo_free_image(&threshIm);
if (stats != NULL ) {
cpl_free(stats);
stats = NULL;
}
sinfo_free_image(&medIm);
sinfo_free_image(&medImage);
sinfo_free_image(&colImage);
sinfo_free_imagelist(&image_list);
sinfo_free_frameset(&raw);
sinfo_badnorm_free(&cfg);
return -1;
}
/**@}*/
|