1 /*
2 filter.c
3
4 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5 *
6 * Part of: SExtractor
7 *
8 * Author: E.BERTIN (IAP)
9 *
10 * Contents: functions dealing with on-line filtering of the image
11 * (for detection).
12 *
13 * Last modify: 26/11/2003
14 *
15 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
16 */
17
18 #ifdef HAVE_CONFIG_H
19 #include "config.h"
20 #endif
21
22 #include <math.h>
23 #include <stdio.h>
24 #include <stdlib.h>
25 #include <string.h>
26
27 #include "define.h"
28 #include "globals.h"
29 #include "prefs.h"
30 #include "fits/fitscat.h"
31 #include "bpro.h"
32 #include "filter.h"
33 #include "image.h"
34
35 filterstruct *thefilter;
36
37 /******************************** convolve ***********************************/
38 /*
39 Convolve a scan line with an array.
40 */
convolve(picstruct * field,PIXTYPE * mscan)41 void convolve(picstruct *field, PIXTYPE *mscan)
42
43 {
44 int mw,mw2,m0,me,m,mx,dmx, y0,dy, sw,sh;
45 float *mask;
46 PIXTYPE *mscane, *s,*s0, *d,*de, mval;
47
48 sw = field->width;
49 sh = field->stripheight;
50 mw = thefilter->convw;
51 mw2 = mw/2;
52 mscane = mscan+sw;
53 y0 = field->y - (thefilter->convh/2);
54 if ((dy = field->ymin-y0) > 0)
55 {
56 m0 = mw*dy;
57 y0 = field->ymin;
58 }
59 else
60 m0 = 0;
61
62 if ((dy = field->ymax - y0) < thefilter->convh)
63 me = mw*dy;
64 else
65 me = mw*thefilter->convh;
66
67 memset(mscan, 0, sw*sizeof(PIXTYPE));
68 s0 = NULL; /* To avoid gcc -Wall warnings */
69 mask = thefilter->conv+m0;
70 for (m = m0, mx = 0; m<me; m++, mx++)
71 {
72 if (mx==mw)
73 mx = 0;
74 if (!mx)
75 s0 = field->strip+sw*((y0++)%sh);
76
77 if ((dmx = mx-mw2)>=0)
78 {
79 s = s0 + dmx;
80 d = mscan;
81 de = mscane - dmx;
82 }
83 else
84 {
85 s = s0;
86 d = mscan - dmx;
87 de = mscane;
88 }
89
90 mval = *(mask++);
91 while (d<de)
92 *(d++) += mval**(s++);
93 }
94
95 return;
96 }
97
98
99 /********************************* getconv **********************************/
100 /*
101 Read the convolution mask from a file.
102 */
getconv(char * filename)103 int getconv(char *filename)
104
105 {
106 FILE *file;
107 char str[MAXCHAR], *sstr, *null = NULL;
108 double sum, var, pix;
109 int i,j, n, normflag;
110
111
112 /* Open the file which may contain a convolution mask */
113
114 if (!(file = fopen(filename, "r")))
115 error(EXIT_FAILURE, "*Error*: cannot open ", filename);
116
117 /* Check it is a convolution mask. Otherwise, exit */
118 fgets(str,MAXCHAR,file);
119 if (strncmp(str,"CONV",4))
120 {
121 fclose(file);
122 return RETURN_ERROR;
123 }
124
125 if (strstr(str,"NORM"))
126 normflag = strstr(str,"NONORM")?0:1;
127 else
128 {
129 warning("No normalization info found in convolution file (old format?)\n",
130 "> => I will assume you want the mask to be normalized...");
131 normflag = 1;
132 }
133 /* Allocate memory for storing mask elements */
134 QMALLOC(thefilter->conv, float, MAXMASK);
135
136 for (i=0, n=0; fgets(str,MAXCHAR,file);)
137 {
138 j=0;
139 sstr = strtok(str, " \t\n");
140 if (sstr && sstr[0]!=(char)'#')
141 do
142 {
143 j++;
144 thefilter->conv[i++] = (float)atof(sstr);
145 if (i>MAXMASK)
146 error(EXIT_FAILURE, "*Error*: Convolution mask too large in ",
147 filename);
148 } while ((sstr = strtok(null, " \t\n")));
149 if (j>n)
150 n = j;
151 }
152
153 fclose(file);
154
155 if ((thefilter->convw = n)<1)
156 error(EXIT_FAILURE, "*Error*: unappropriate convolution mask width in ",
157 filename);
158 if (i%n)
159 error(EXIT_FAILURE, "*Error*: unappropriate convolution mask line in ",
160 filename);
161
162 QREALLOC(thefilter->conv, float, i);
163
164 thefilter->convh = i/n;
165 var = 0.0;
166 for (j=0, sum=0.0; j<i; j++)
167 {
168 sum += fabs(pix = thefilter->conv[j]);
169 var += pix*pix;
170 }
171
172 thefilter->varnorm = (float)(var=sqrt(var));
173
174 if (normflag)
175 {
176 if (sum == 0.0)
177 {
178 warning("Zero-sum filter: ", "Normalization switched to variance-mode");
179 sum = var;
180 }
181 for (j=0; j<i; j++)
182 thefilter->conv[j] /= sum;
183 }
184
185 thefilter->nconv = thefilter->convw*thefilter->convh;
186
187 return RETURN_OK;
188 }
189
190
191 /********************************* getfilter *********************************/
192 /*
193 Read either a convolution mask or an ANN file.
194 */
getfilter(char * filename)195 void getfilter(char *filename)
196
197 {
198 QCALLOC(thefilter, filterstruct, 1);
199 if (getconv(filename) != RETURN_OK && getneurfilter(filename) != RETURN_OK)
200 error(EXIT_FAILURE, "*Error*: not a suitable filter in ", filename);
201
202 return;
203 }
204
205
206 /******************************** getneurfilter ******************************/
207 /*
208 Read an ANN RETINA-filter file.
209 */
getneurfilter(char * filename)210 int getneurfilter(char *filename)
211
212 {
213 #define FILTEST(x) \
214 if (x != RETURN_OK) \
215 error(EXIT_FAILURE, "*Error*: incorrect filter header in ", filename)
216
217 catstruct *fcat;
218 tabstruct *ftab;
219 int ival;
220
221 /* We first map the catalog */
222 if (!(fcat = read_cat(filename)))
223 return RETURN_ERROR;
224 /* Test if the requested table is present */
225 if (!(ftab = name_to_tab(fcat, "BP-ANN", 0)))
226 error(EXIT_FAILURE, "*Error*: no BP-ANN info found in ", filename);
227 FILTEST(fitsread(ftab->headbuf, "BPTYPE ", gstr,H_STRING,T_STRING));
228 if (strcmp(gstr, "RETINA"))
229 error(EXIT_FAILURE, "*Error*: not a retina-filter in ", filename);
230 FILTEST(fitsread(ftab->headbuf, "RENAXIS ", &ival ,H_INT, T_LONG));
231 if (ival != 2)
232 error(EXIT_FAILURE, "*Error*: not a 2D retina in ", filename);
233 FILTEST(fitsread(ftab->headbuf, "RENAXIS1", &thefilter->convw ,H_INT,T_LONG));
234 FILTEST(fitsread(ftab->headbuf, "RENAXIS2", &thefilter->convh ,H_INT,T_LONG));
235 thefilter->nconv = thefilter->convw*thefilter->convh;
236 QMALLOC(thefilter->conv, float, thefilter->nconv);
237 thefilter->bpann = loadtab_bpann(ftab, filename);
238
239 close_cat(fcat);
240 free_cat(&fcat,1);
241
242 return RETURN_OK;
243 }
244
245
246 /********************************* endfilter *********************************/
247 /*
248 Terminate filtering procedures.
249 */
endfilter()250 void endfilter()
251
252 {
253 QFREE(thefilter->conv);
254 if (thefilter->bpann)
255 {
256 free_bpann(thefilter->bpann);
257 thefilter->bpann = NULL;
258 }
259
260 QFREE(thefilter);
261
262 return;
263 }
264
265
266 /********************************** filter ***********************************/
267 /*
268 Switch to the appropriate filtering routine.
269 */
filter(picstruct * field,PIXTYPE * mscan)270 void filter(picstruct *field, PIXTYPE *mscan)
271
272 {
273 if (thefilter->bpann)
274 neurfilter(field, mscan);
275 else
276 convolve(field, mscan);
277
278 return;
279 }
280
281
282 /******************************** neurfilter *********************************/
283 /*
284 Filter a scan line using an artificial retina.
285 */
neurfilter(picstruct * field,PIXTYPE * mscan)286 void neurfilter(picstruct *field, PIXTYPE *mscan)
287
288 {
289 PIXTYPE cval;
290 float *pix, resp, sig, threshlow, threshhigh;
291 int i,x, tflag;
292
293 sig = field->backsig;
294 if ((tflag = (prefs.nfilter_thresh>0)))
295 {
296 threshlow = prefs.filter_thresh[0]*field->backsig;
297 threshhigh = (prefs.nfilter_thresh<2)?
298 BIG : prefs.filter_thresh[1]*field->backsig;
299 }
300 else
301 threshlow = threshhigh = 0.0; /* To avoid gcc -Wall warnings */
302
303 for (x=0;x<field->width;x++)
304 {
305 if (tflag)
306 {
307 cval = PIX(field, x, field->y);
308 if (cval<threshlow || cval>threshhigh)
309 {
310 *(mscan++) = cval;
311 continue;
312 }
313 }
314 /*-- Copy the surrounding image area to the retina */
315 copyimage(field, thefilter->conv, thefilter->convw, thefilter->convh,
316 x, field->y);
317 pix = thefilter->conv;
318 /*-- Apply a transform of the intensity scale */
319 for (i=thefilter->nconv; i--; pix++)
320 *pix = *pix>0.0?log(1+*pix/sig):-log(1-*pix/sig);
321 play_bpann(thefilter->bpann, thefilter->conv, &resp);
322 if (resp>70.0)
323 resp = 70.0;
324 else if (resp<-70.0)
325 resp = -70.0;
326 /*-- Recover the linear intensity scale */
327 *(mscan++) = resp>0.0?sig*(exp(resp)-1):sig*(exp(-resp)-1);
328 }
329
330 return;
331 }
332
333
334 /******************************* convolve_image ******************************/
335 /*
336 Convolve a vignet with an array.
337 */
convolve_image(picstruct * field,double * vig1,double * vig2,int width,int height)338 void convolve_image(picstruct *field, double *vig1,
339 double *vig2, int width, int height)
340
341 {
342 int mw,mw2,m0,me,m,mx,dmx, y, y0,dy;
343 float *mask;
344 double *mscane, *s,*s0, *vig3, *d,*de, mval, val;
345
346 /* If no filtering, just return a copy of the input data */
347 if (!thefilter || thefilter->bpann)
348 {
349 memcpy(vig2, vig1, width*height*sizeof(double));
350 return;
351 }
352 s0 = NULL; /* To avoid gcc -Wall warnings */
353 mw = thefilter->convw;
354 mw2 = mw/2;
355 memset(vig2, 0, width*height*sizeof(double));
356 vig3 = vig2;
357 for (y=0; y<height; y++, vig3+=width)
358 {
359 mscane = vig3+width;
360 y0 = y - (thefilter->convh/2);
361 if (y0 < 0)
362 {
363 m0 = -mw*y0;
364 y0 = 0;
365 }
366 else
367 m0 = 0;
368 if ((dy = height - y0) < thefilter->convh)
369 me = mw*dy;
370 else
371 me = mw*thefilter->convh;
372 mask = thefilter->conv+m0;
373 for (m = m0, mx = 0; m<me; m++, mx++)
374 {
375 if (mx==mw)
376 mx = 0;
377 if (!mx)
378 s0 = vig1+width*(y0++);
379 if ((dmx = mx-mw2)>=0)
380 {
381 s = s0 + dmx;
382 d = vig3;
383 de = mscane - dmx;
384 }
385 else
386 {
387 s = s0;
388 d = vig3 - dmx;
389 de = mscane;
390 }
391
392 mval = *(mask++);
393 while (d<de)
394 *(d++) += ((val = *(s++))>-BIG? mval*val:0.0);
395 }
396 }
397
398 return;
399 }
400
401
402