1 /*
2 scan.c
3
4 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5 *
6 * Part of: SExtractor
7 *
8 * Author: E.BERTIN (IAP)
9 *
10 * Contents: functions for extraction of connected pixels from
11 * a pixmap.
12 *
13 * Last modify: 29/11/2005
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 "back.h"
31 #include "check.h"
32 #include "clean.h"
33 #include "extract.h"
34 #include "filter.h"
35 #include "image.h"
36 #include "plist.h"
37
38 /****************************** scanimage ************************************
39 PROTO void scanimage(picstruct *field, picstruct *dfield, picstruct *ffield,
40 picstruct *wfield, picstruct *dwfield)
41 PURPOSE Scan of the large pixmap(s). Main loop and heart of the program.
42 INPUT Measurement field pointer,
43 Detection field pointer,
44 Flag field pointer,
45 Measurement weight-map field pointer,
46 Detection weight-map field pointer,
47 OUTPUT -.
48 NOTES -.
49 AUTHOR E. Bertin (IAP)
50 VERSION 29/11/2005
51 ***/
52
53 PIXTYPE *dumscan;
54
scanimage(picstruct * field,picstruct * dfield,picstruct ** pffield,int nffield,picstruct * wfield,picstruct * dwfield)55 void scanimage(picstruct *field, picstruct *dfield, picstruct **pffield,
56 int nffield, picstruct *wfield, picstruct *dwfield)
57
58 {
59 static infostruct curpixinfo, *info, *store,
60 initinfo, freeinfo, *victim;
61 picstruct *ffield;
62 checkstruct *check;
63 objliststruct objlist;
64 objstruct *cleanobj;
65 pliststruct *pixel, *pixt;
66 picstruct *cfield, *cdwfield;
67
68 char *marker, newmarker, *blankpad, *bpt,*bpt0;
69 int co, i,j, flag, luflag,pstop, xl,xl2,yl, cn,
70 nposize, stacksize, w, h, blankh, maxpixnb,
71 varthreshflag;
72 short trunflag;
73 PIXTYPE thresh, relthresh, cdnewsymbol, cdvar,
74 *scan,*dscan,*cdscan,*dwscan,*cdwscan,*cdwscanp,
75 *scant, *wscan,*wscanp;
76 FLAGTYPE *pfscan[MAXFLAG];
77 status cs, ps, *psstack;
78 int *start, *end, ymax;
79
80 /* Avoid gcc -Wall warnings */
81 scan = dscan = cdscan = cdwscan = cdwscanp = wscan = wscanp = NULL;
82 victim = NULL; /* Avoid gcc -Wall warnings */
83 blankh = 0; /* Avoid gcc -Wall warnings */
84 /*----- Beginning of the main loop: Initialisations */
85 thecat.ntotal = thecat.ndetect = 0;
86
87 /* cfield is the detection field in any case */
88 cfield = dfield? dfield:field;
89
90 /* cdwfield is the detection weight-field if available */
91 cdwfield = dwfield? dwfield:(prefs.dweight_flag?wfield:NULL);
92
93 /* If WEIGHTing and no absolute thresholding, activate threshold scaling */
94 varthreshflag = (cdwfield && prefs.thresh_type[0]!=THRESH_ABSOLUTE);
95 relthresh = varthreshflag ? prefs.dthresh[0] : 0.0;/* To avoid gcc warnings*/
96 w = cfield->width;
97 h = cfield->height;
98 objlist.dthresh = cfield->dthresh;
99 objlist.thresh = cfield->thresh;
100 cfield->yblank = 1;
101 field->y = field->stripy = 0;
102 field->ymin = field->stripylim = 0;
103 field->stripysclim = 0;
104 if (dfield)
105 {
106 dfield->y = dfield->stripy = 0;
107 dfield->ymin = dfield->stripylim = 0;
108 dfield->stripysclim = 0;
109 }
110 if (nffield)
111 for (i=0; i<nffield; i++)
112 {
113 ffield = pffield[i];
114 ffield->y = ffield->stripy = 0;
115 ffield->ymin = ffield->stripylim = 0;
116 ffield->stripysclim = 0;
117 }
118 if (wfield)
119 {
120 wfield->y = wfield->stripy = 0;
121 wfield->ymin = wfield->stripylim = 0;
122 wfield->stripysclim = 0;
123 }
124 if (dwfield)
125 {
126 dwfield->y = dwfield->stripy = 0;
127 dwfield->ymin = dwfield->stripylim = 0;
128 dwfield->stripysclim = 0;
129 }
130
131 /*Allocate memory for buffers */
132 stacksize = w+1;
133 QMALLOC(info, infostruct, stacksize);
134 QMALLOC(store, infostruct, stacksize);
135 QMALLOC(marker, char, stacksize);
136 QMALLOC(dumscan, PIXTYPE, stacksize);
137 QMALLOC(psstack, status, stacksize);
138 QMALLOC(start, int, stacksize);
139 QMALLOC(end, int, stacksize);
140 blankpad = bpt = NULL;
141 lutzalloc(w,h);
142 allocparcelout();
143 if (prefs.filter_flag)
144 {
145 QMALLOC(cdscan, PIXTYPE, stacksize);
146 if (cdwfield)
147 {
148 QCALLOC(cdwscan, PIXTYPE, stacksize);
149 if (PLISTEXIST(wflag))
150 {
151 QMALLOC(cdwscanp, PIXTYPE, stacksize);
152 }
153 }
154 /*-- One needs a buffer to protect filtering if source-blanking applies */
155 if (prefs.blank_flag)
156 {
157 blankh = thefilter->convh/2+1;
158 QMALLOC(blankpad, char, w*blankh);
159 cfield->yblank -= blankh;
160 if (dfield)
161 field->yblank = cfield->yblank;
162 bpt = blankpad;
163 }
164 }
165
166 /* Some initializations */
167
168 thresh = objlist.dthresh;
169 initinfo.pixnb = 0;
170 initinfo.flag = 0;
171 initinfo.firstpix = initinfo.lastpix = -1;
172
173 for (xl=0; xl<stacksize; xl++)
174 {
175 marker[xl] = 0 ;
176 dumscan[xl] = -BIG ;
177 }
178
179 co = pstop = 0;
180 objlist.nobj = 1;
181 curpixinfo.pixnb = 1;
182
183 /* Init cleaning procedure */
184 initclean();
185
186 /*----- Allocate memory for the pixel list */
187 init_plist();
188
189 if (!(pixel = objlist.plist = malloc(nposize=prefs.mem_pixstack*plistsize)))
190 error(EXIT_FAILURE, "Not enough memory to store the pixel stack:\n",
191 " Try to decrease MEMORY_PIXSTACK");
192
193 /*----- at the beginning, "free" object fills the whole pixel list */
194 freeinfo.firstpix = 0;
195 freeinfo.lastpix = nposize-plistsize;
196 pixt = pixel;
197 for (i=plistsize; i<nposize; i += plistsize, pixt += plistsize)
198 PLIST(pixt, nextpix) = i;
199 PLIST(pixt, nextpix) = -1;
200
201 /*----- Here we go */
202 for (yl=0; yl<=h;)
203 {
204 ps = COMPLETE;
205 cs = NONOBJECT;
206 if (yl==h)
207 {
208 /*---- Need an empty line for Lutz' algorithm to end gracely */
209 if (prefs.filter_flag)
210 {
211 free(cdscan);
212 if (cdwfield)
213 {
214 free(cdwscan);
215 if (PLISTEXIST(wflag))
216 free(cdwscanp);
217 }
218 }
219 cdwscan = cdwscanp = cdscan = dumscan;
220 }
221 else
222 {
223 if (nffield)
224 for (i=0; i<nffield; i++)
225 {
226 ffield = pffield[i];
227 pfscan[i] = (ffield->stripy==ffield->stripysclim)?
228 (FLAGTYPE *)loadstrip(ffield, (picstruct *)NULL)
229 : &ffield->fstrip[ffield->stripy*ffield->width];
230 }
231 if (wfield)
232 {
233 /*------ Copy the previous weight line to track bad pixel limits */
234 wscan = (wfield->stripy==wfield->stripysclim)?
235 (PIXTYPE *)loadstrip(wfield, (picstruct *)NULL)
236 : &wfield->strip[wfield->stripy*wfield->width];
237 if (yl && PLISTEXIST(wflag))
238 wscanp = &wfield->strip[((yl-1)%wfield->stripheight)*wfield->width];
239 }
240 else
241 wscan = NULL;
242 scan = (field->stripy==field->stripysclim)?
243 (PIXTYPE *)loadstrip(field, wfield)
244 : &field->strip[field->stripy*field->width];
245 if (dwfield)
246 dwscan = (dwfield->stripy==dwfield->stripysclim)?
247 (PIXTYPE *)loadstrip(dwfield,
248 dfield?(picstruct *)NULL:dwfield)
249 : &dwfield->strip[dwfield->stripy*dwfield->width];
250 else
251 dwscan = wscan;
252 if (dfield)
253 dscan = (dfield->stripy==dfield->stripysclim)?
254 (PIXTYPE *)loadstrip(dfield, dwfield)
255 : &dfield->strip[dfield->stripy*dfield->width];
256 else
257 dscan = scan;
258
259 if (PLISTEXIST(wflag) && cdwfield)
260 /*------ Copy the previous filtered weight line to track bad pixel limits */
261 memcpy(cdwscanp, cdwscan, wfield->width*sizeof(PIXTYPE));
262 if (prefs.filter_flag)
263 {
264 filter(cfield, cdscan);
265 if (cdwfield)
266 filter(cdwfield, cdwscan);
267 }
268 else
269 {
270 cdscan = dscan;
271 cdwscan = dwscan;
272 }
273
274 if ((check=prefs.check[CHECK_FILTERED]))
275 writecheck(check, cdscan, w);
276 }
277
278 trunflag = (yl==0 || yl==h-1)? OBJ_TRUNC:0;
279
280 for (xl=0; xl<=w; xl++)
281 {
282 if (xl == w)
283 cdnewsymbol = -BIG;
284 else
285 cdnewsymbol = cdscan[xl];
286
287 newmarker = marker[xl];
288 marker[xl] = 0;
289
290 curpixinfo.flag = trunflag;
291 if (varthreshflag)
292 thresh = relthresh*sqrt(cdvar = ((xl==w || yl==h)? 0.0:cdwscan[xl]));
293 luflag = cdnewsymbol > thresh?1:0;
294
295 if (luflag)
296 {
297 if (xl==0 || xl==w-1)
298 curpixinfo.flag |= OBJ_TRUNC;
299 pixt = pixel + (cn=freeinfo.firstpix);
300 freeinfo.firstpix = PLIST(pixt, nextpix);
301
302 /*------- Running out of pixels, the largest object becomes a "victim" ------*/
303
304 if (freeinfo.firstpix==freeinfo.lastpix)
305 {
306 sprintf(gstr, "%d,%d", xl+1, yl+1);
307 warning("Pixel stack overflow at position ", gstr);
308 maxpixnb = 0;
309 for (i=0; i<=w; i++)
310 if (store[i].pixnb>maxpixnb)
311 if (marker[i]=='S' || (newmarker=='S' && i==xl))
312 {
313 flag = 0;
314 if (i<xl)
315 for (j=0; j<=co; j++)
316 flag |= (start[j]==i);
317 if (!flag)
318 maxpixnb = (victim = &store[i])->pixnb;
319 }
320 for (j=1; j<=co; j++)
321 if (info[j].pixnb>maxpixnb)
322 maxpixnb = (victim = &info[j])->pixnb;
323
324 if (!maxpixnb)
325 error(EXIT_FAILURE, "*Fatal Error*: something is badly bugged in ",
326 "scanimage()!");
327 if (maxpixnb <= 1)
328 error(EXIT_FAILURE, "Pixel stack overflow in ", "scanimage()");
329 freeinfo.firstpix = PLIST(pixel+victim->firstpix, nextpix);
330 PLIST(pixel+victim->lastpix, nextpix) = freeinfo.lastpix;
331 PLIST(pixel+(victim->lastpix=victim->firstpix), nextpix) = -1;
332 victim->pixnb = 1;
333 victim->flag |= OBJ_OVERFLOW;
334 }
335
336 /*---------------------------------------------------------------------------*/
337 curpixinfo.lastpix = curpixinfo.firstpix = cn;
338 PLIST(pixt, nextpix) = -1;
339 PLIST(pixt, x) = xl;
340 PLIST(pixt, y) = yl;
341 PLIST(pixt, value) = scan[xl];
342 if (PLISTEXIST(dvalue))
343 PLISTPIX(pixt, dvalue) = dscan[xl];
344 if (PLISTEXIST(cdvalue))
345 PLISTPIX(pixt, cdvalue) = cdnewsymbol;
346 if (PLISTEXIST(flag))
347 for (i=0; i<nffield; i++)
348 PLISTFLAG(pixt, flag[i]) = pfscan[i][xl];
349 if (PLISTEXIST(wflag) && wscan)
350 {
351 PLISTFLAG(pixt, wflag) = 0;
352 if (xl>0)
353 {
354 if (wscan[xl-1] >= BIG)
355 PLISTFLAG(pixt, wflag) |= OBJ_WEIGHTZERO;
356 if (cdwscan[xl-1] >= BIG)
357 PLISTFLAG(pixt, wflag) = OBJ_DWEIGHTZERO;
358 }
359 PLISTFLAG(pixt, wflag) = 0;
360 if (xl<w-1)
361 {
362 if (wscan[xl+1] >= BIG)
363 PLISTFLAG(pixt, wflag) |= OBJ_WEIGHTZERO;
364 if (cdwscan[xl+1] >= BIG)
365 PLISTFLAG(pixt, wflag) = OBJ_DWEIGHTZERO;
366 }
367 if (yl>0)
368 {
369 if (wscanp[xl] >= BIG)
370 PLISTFLAG(pixt, wflag) |= OBJ_WEIGHTZERO;
371 if (cdwscanp[xl] >= BIG)
372 PLISTFLAG(pixt, wflag) = OBJ_DWEIGHTZERO;
373 }
374 }
375 if (PLISTEXIST(dthresh))
376 PLISTPIX(pixt, dthresh) = thresh;
377 if (PLISTEXIST(var))
378 PLISTPIX(pixt, var) = wscan[xl];
379
380 if (cs != OBJECT)
381 /*------------------------------- Start Segment -----------------------------*/
382
383 {
384 cs = OBJECT;
385 if (ps == OBJECT)
386 {
387 if (start[co] == UNKNOWN)
388 {
389 marker[xl] = 'S';
390 start[co] = xl;
391 }
392 else
393 marker[xl] = 's';
394 }
395 else
396 {
397 psstack[pstop++] = ps;
398 marker[xl] = 'S';
399 start[++co] = xl;
400 ps = COMPLETE;
401 info[co] = initinfo;
402 }
403 }
404
405 /*---------------------------------------------------------------------------*/
406 }
407
408 if (newmarker)
409
410 /*---------------------------- Process New Marker ---------------------------*/
411
412 {
413 if (newmarker == 'S')
414 {
415 psstack[pstop++] = ps;
416 if (cs == NONOBJECT)
417 {
418 psstack[pstop++] = COMPLETE;
419 info[++co] = store[xl];
420 start[co] = UNKNOWN;
421 }
422 else
423 update (&info[co],&store[xl], pixel);
424 ps = OBJECT;
425 }
426 else if (newmarker == 's')
427 {
428 if ((cs == OBJECT) && (ps == COMPLETE))
429 {
430 pstop--;
431 xl2 = start[co];
432 update (&info[co-1],&info[co], pixel);
433 if (start[--co] == UNKNOWN)
434 start[co] = xl2;
435 else
436 marker[xl2] = 's';
437 }
438 ps = OBJECT;
439 }
440 else if (newmarker == 'f')
441 ps = INCOMPLETE;
442 else if (newmarker == 'F')
443 {
444 ps = psstack[--pstop];
445 if ((cs == NONOBJECT) && (ps == COMPLETE))
446 {
447 if (start[co] == UNKNOWN)
448 {
449 if ((int)info[co].pixnb >= prefs.ext_minarea)
450 {
451 sortit(field, dfield, wfield, cdwfield, &info[co], &objlist,
452 cdwscan, wscan);
453 }
454 /* ------------------------------------ free the chain-list */
455
456 PLIST(pixel+info[co].lastpix, nextpix) = freeinfo.firstpix;
457 freeinfo.firstpix = info[co].firstpix;
458 }
459 else
460 {
461 marker[end[co]] = 'F';
462 store[start[co]] = info[co];
463 }
464 co--;
465 ps = psstack[--pstop];
466 }
467 }
468 }
469 /*---------------------------------------------------------------------------*/
470
471 if (luflag)
472 update (&info[co],&curpixinfo, pixel);
473 else
474 {
475 if (cs == OBJECT)
476 /*-------------------------------- End Segment ------------------------------*/
477 {
478 cs = NONOBJECT;
479 if (ps != COMPLETE)
480 {
481 marker[xl] = 'f';
482 end[co] = xl;
483 }
484 else
485 {
486 ps = psstack[--pstop];
487 marker[xl] = 'F';
488 store[start[co]] = info[co];
489 co--;
490 }
491 }
492 }
493
494 if (prefs.blank_flag && xl<w)
495 {
496 if (prefs.filter_flag)
497 *(bpt++) = (luflag)?1:0;
498 else if (luflag)
499 dscan[xl] = -BIG;
500 if (dfield && luflag)
501 scan[xl] = -BIG;
502 }
503 /*--------------------- End of the loop over the x's -----------------------*/
504 }
505
506 /* Detected pixel removal at the end of each line */
507 if (prefs.blank_flag && yl<h)
508 {
509 if (prefs.filter_flag)
510 {
511 bpt = bpt0 = blankpad + w*((yl+1)%blankh);
512 if (cfield->yblank >= 0)
513 {
514 scant = &PIX(cfield, 0, cfield->yblank);
515 for (i=w; i--; scant++)
516 if (*(bpt++))
517 *scant = -BIG;
518 if (dfield)
519 {
520 bpt = bpt0;
521 scant = &PIX(field, 0, cfield->yblank);
522 for (i=w; i--; scant++)
523 if (*(bpt++))
524 *scant = -BIG;
525 }
526 bpt = bpt0;
527 }
528 }
529 cfield->yblank++;
530 if (dfield)
531 field->yblank = cfield->yblank;
532 }
533
534 /*-- Prepare markers for the next line */
535 yl++;
536 field->stripy = (field->y=yl)%field->stripheight;
537 if (dfield)
538 dfield->stripy = (dfield->y=yl)%dfield->stripheight;
539 if (nffield)
540 for (i=0; i<nffield; i++)
541 {
542 ffield = pffield[i];
543 ffield->stripy = (ffield->y=yl)%ffield->stripheight;
544 }
545 if (wfield)
546 wfield->stripy = (wfield->y=yl)%wfield->stripheight;
547 if (dwfield)
548 dwfield->stripy = (dwfield->y=yl)%dwfield->stripheight;
549
550 /*-- Remove objects close to the ymin limit if ymin is ready to increase */
551 if (cfield->stripy==cfield->stripysclim)
552 {
553 cleanobj = cleanobjlist->obj+cleanobjlist->nobj-1;
554 for (i=cleanobjlist->nobj; i--; cleanobj--)
555 {
556 if (cleanobj->ycmin <= cfield->ymin)
557 {
558 /*-------- Warn if there is a possibility for any aperture to be truncated */
559 if ((ymax=cleanobj->ycmax) > cfield->ymax)
560 {
561 sprintf(gstr, "Object at position %.0f,%.0f ",
562 cleanobj->mx+1, cleanobj->my+1);
563 QWARNING(gstr, "may have some apertures truncated:\n"
564 " You might want to increase MEMORY_BUFSIZE");
565 }
566 else if (ymax>cfield->yblank && prefs.blank_flag)
567 {
568 sprintf(gstr, "Object at position %.0f,%.0f ",
569 cleanobj->mx+1, cleanobj->my+1);
570 QWARNING(gstr, "may have some unBLANKed neighbours:\n"
571 " You might want to increase MEMORY_PIXSTACK");
572 }
573 endobject(field, dfield, wfield, cdwfield, i, cleanobjlist);
574 subcleanobj(i);
575 cleanobj = cleanobjlist->obj+i; /* realloc in subcleanobj() */
576 }
577 }
578 }
579
580 if (!((yl+1)%16))
581 NPRINTF(OUTPUT, "\33[1M> Line:%5d "
582 "Objects: %8d detected / %8d sextracted\n\33[1A",
583 yl+1, thecat.ndetect, thecat.ntotal);
584 /*--------------------- End of the loop over the y's -----------------------*/
585 }
586
587 /* Removal or the remaining pixels */
588 if (prefs.blank_flag && prefs.filter_flag && (cfield->yblank >= 0))
589 for (j=blankh-1; j--; yl++)
590 {
591 bpt = bpt0 = blankpad + w*(yl%blankh);
592 scant = &PIX(cfield, 0, cfield->yblank);
593 for (i=w; i--; scant++)
594 if (*(bpt++))
595 *scant = -BIG;
596 if (dfield)
597 {
598 bpt = bpt0;
599 scant = &PIX(field, 0, cfield->yblank);
600 for (i=w; i--; scant++)
601 if (*(bpt++))
602 *scant = -BIG;
603 }
604 cfield->yblank++;
605 if (dfield)
606 field->yblank = cfield->yblank;
607 }
608
609 /* Now that all "detected" pixels have been removed, analyse detections */
610 for (j=cleanobjlist->nobj; j--;)
611 {
612 endobject(field, dfield, wfield, cdwfield, 0, cleanobjlist);
613 subcleanobj(0);
614 }
615
616 endclean();
617
618 /*Free memory */
619 freeparcelout();
620 free(pixel);
621 lutzfree();
622 free(info);
623 free(store);
624 free(marker);
625 free(dumscan);
626 free(psstack);
627 free(start);
628 free(end);
629 if (prefs.blank_flag && prefs.filter_flag)
630 free(blankpad);
631
632 return;
633 }
634
635
636 /********************************* update ************************************/
637 /*
638 update object's properties each time one of its pixels is scanned by lutz()
639 */
update(infostruct * infoptr1,infostruct * infoptr2,pliststruct * pixel)640 void update(infostruct *infoptr1, infostruct *infoptr2, pliststruct *pixel)
641
642 {
643 infoptr1->pixnb += infoptr2->pixnb;
644 infoptr1->flag |= infoptr2->flag;
645 if (infoptr1->firstpix == -1)
646 {
647 infoptr1->firstpix = infoptr2->firstpix;
648 infoptr1->lastpix = infoptr2->lastpix;
649 }
650 else if (infoptr2->lastpix != -1)
651 {
652 PLIST(pixel+infoptr1->lastpix, nextpix) = infoptr2->firstpix;
653 infoptr1->lastpix = infoptr2->lastpix;
654 }
655
656 return;
657 }
658
659 /********************************* sortit ************************************/
660 /*
661 build the object structure.
662 */
sortit(picstruct * field,picstruct * dfield,picstruct * wfield,picstruct * dwfield,infostruct * info,objliststruct * objlist,PIXTYPE * cdwscan,PIXTYPE * wscan)663 void sortit(picstruct *field, picstruct *dfield, picstruct *wfield,
664 picstruct *dwfield, infostruct *info, objliststruct *objlist,
665 PIXTYPE *cdwscan, PIXTYPE *wscan)
666
667 {
668 picstruct *cfield;
669 objliststruct objlistout, *objlist2;
670 static objstruct obj;
671 objstruct *cobj;
672 pliststruct *pixel;
673 int i,j,n;
674
675 cfield = dfield? dfield: field;
676
677 pixel = objlist->plist;
678 objlistout.obj = NULL;
679 objlistout.plist = NULL;
680 objlistout.nobj = objlistout.npix = 0;
681
682 /*----- Allocate memory to store object data */
683
684 objlist->obj = &obj;
685 objlist->nobj = 1;
686
687 memset(&obj, 0, (size_t)sizeof(objstruct));
688 objlist->npix = info->pixnb;
689 obj.firstpix = info->firstpix;
690 obj.lastpix = info->lastpix;
691 obj.flag = info->flag;
692 obj.dthresh = objlist->dthresh;
693 obj.thresh = objlist->thresh;
694
695 preanalyse(0, objlist, ANALYSE_FAST);
696
697 /*----- Check if the current strip contains the lower isophote... */
698 if ((int)obj.ymin < cfield->ymin)
699 obj.flag |= OBJ_ISO_PB;
700
701 if (!(obj.flag & OBJ_OVERFLOW) && (createsubmap(objlist, 0) == RETURN_OK))
702 {
703 if (parcelout(objlist, &objlistout) == RETURN_OK)
704 objlist2 = &objlistout;
705 else
706 {
707 objlist2 = objlist;
708 for (i=0; i<objlist2->nobj; i++)
709 objlist2->obj[i].flag |= OBJ_DOVERFLOW;
710 sprintf(gstr, "%.0f,%.0f", obj.mx+1, obj.my+1);
711 warning("Deblending overflow for detection at ", gstr);
712 }
713 free(obj.submap);
714 }
715 else
716 objlist2 = objlist;
717
718 for (i=0; i<objlist2->nobj; i++)
719 {
720 preanalyse(i, objlist2, ANALYSE_FULL|ANALYSE_ROBUST);
721 analyse(field, dfield, i, objlist2);
722 cobj = objlist2->obj + i;
723 if (prefs.blank_flag)
724 {
725 if (createblank(objlist2,i) != RETURN_OK)
726 {
727 /*------ Not enough mem. for the BLANK vignet: flag the object now */
728 cobj->flag |= OBJ_OVERFLOW;
729 cobj->blank = cobj->dblank = NULL;
730 sprintf(gstr, "%.0f,%.0f", cobj->mx+1, cobj->my+1);
731 warning("Memory overflow during masking for detection at ", gstr);
732 }
733 }
734
735 if ((n=cleanobjlist->nobj) >= prefs.clean_stacksize)
736 {
737 objstruct *cleanobj;
738 int ymin, ymax, victim=0;
739
740 ymin = 2000000000; /* No image is expected to be that tall ! */
741 cleanobj = cleanobjlist->obj;
742 for (j=0; j<n; j++, cleanobj++)
743 if (cleanobj->ycmax < ymin)
744 {
745 victim = j;
746 ymin = cleanobj->ycmax;
747 }
748
749 /*---- Warn if there is a possibility for any aperture to be truncated */
750 if (field->ymax < field->height)
751 {
752 cleanobj = &cleanobjlist->obj[victim];
753 if ((ymax=cleanobj->ycmax) > field->ymax)
754 {
755 sprintf(gstr, "Object at position %.0f,%.0f ",
756 cleanobj->mx+1, cleanobj->my+1);
757 QWARNING(gstr, "may have some apertures truncated:\n"
758 " You might want to increase MEMORY_OBJSTACK");
759 }
760 else if (ymax>field->yblank && prefs.blank_flag)
761 {
762 sprintf(gstr, "Object at position %.0f,%.0f ",
763 cleanobj->mx+1, cleanobj->my+1);
764 QWARNING(gstr, "may have some unBLANKed neighbours\n"
765 " You might want to increase MEMORY_OBJSTACK");
766 }
767 }
768
769 endobject(field, dfield, wfield, dwfield, victim, cleanobjlist);
770 subcleanobj(victim);
771 }
772
773 /* Only add the object if it is not swallowed by cleaning */
774 if (!prefs.clean_flag || clean(field, dfield, i, objlist2))
775 addcleanobj(cobj);
776 }
777
778 free(objlistout.plist);
779 free(objlistout.obj);
780
781 return;
782 }
783
784
785 /******************************** preanalyse *********************************
786 PROTO void preanalyse(int no, objliststruct *objlist, int analyse_type)
787 PURPOSE Compute basic image parameters from the pixel-list for each detection.
788 INPUT objlist number,
789 objlist pointer,
790 analysis switch flag.
791 OUTPUT -.
792 NOTES -.
793 AUTHOR E. Bertin (IAP & Leiden & ESO)
794 VERSION 28/11/2003
795 ***/
preanalyse(int no,objliststruct * objlist,int analyse_type)796 void preanalyse(int no, objliststruct *objlist, int analyse_type)
797
798 {
799 objstruct *obj = &objlist->obj[no];
800 pliststruct *pixel = objlist->plist, *pixt;
801 PIXTYPE peak, cpeak, val, cval, minthresh, thresht;
802 double thresh,thresh2, t1t2,darea,
803 mx,my, mx2,my2,mxy, rv, tv,
804 xm,ym, xm2,ym2,xym,
805 temp,temp2, theta,pmx2,pmy2;
806 int x, y, xmin,xmax, ymin,ymax,area2, fdnpix, dnpix;
807
808
809 /*----- initialize stacks and bounds */
810 thresh = obj->dthresh;
811 if (PLISTEXIST(dthresh))
812 minthresh = BIG;
813 else
814 minthresh = 0.0;
815 fdnpix = dnpix = 0;
816 rv = 0.0;
817 peak = cpeak = -BIG;
818 ymin = xmin = 2*MAXPICSIZE; /* to be really sure!! */
819 ymax = xmax = 0;
820
821 /*----- integrate results */
822 for (pixt=pixel+obj->firstpix; pixt>=pixel; pixt=pixel+PLIST(pixt,nextpix))
823 {
824 x = PLIST(pixt, x);
825 y = PLIST(pixt, y);
826 val=PLISTPIX(pixt, dvalue);
827 if (cpeak < (cval=PLISTPIX(pixt, cdvalue)))
828 cpeak = cval;
829 if (PLISTEXIST(dthresh) && (thresht=PLISTPIX(pixt, dthresh))<minthresh)
830 minthresh = thresht;
831 if (peak < val)
832 peak = val;
833 rv += cval;
834 if (xmin > x)
835 xmin = x;
836 if (xmax < x)
837 xmax = x;
838 if (ymin > y)
839 ymin = y;
840 if (ymax < y)
841 ymax = y;
842 fdnpix++;
843 }
844
845 if (PLISTEXIST(dthresh))
846 obj->dthresh = thresh = minthresh;
847
848 /* copy some data to "obj" structure */
849
850 obj->fdnpix = (LONG)fdnpix;
851 obj->fdflux = (float)rv;
852 obj->fdpeak = cpeak;
853 obj->dpeak = peak;
854 obj->xmin = xmin;
855 obj->xmax = xmax;
856 obj->ymin = ymin;
857 obj->ymax = ymax;
858
859 if (analyse_type & ANALYSE_FULL)
860 {
861 mx = my = tv = 0.0;
862 mx2 = my2 = mxy = 0.0;
863 thresh2 = (thresh + peak)/2.0;
864 area2 = 0;
865 for (pixt=pixel+obj->firstpix; pixt>=pixel; pixt=pixel+PLIST(pixt,nextpix))
866 {
867 x = PLIST(pixt,x)-xmin; /* avoid roundoff errors on big images */
868 y = PLIST(pixt,y)-ymin; /* avoid roundoff errors on big images */
869 cval = PLISTPIX(pixt, cdvalue);
870 tv += (val = PLISTPIX(pixt, dvalue));
871 if (val>thresh)
872 dnpix++;
873 if (val > thresh2)
874 area2++;
875 mx += cval * x;
876 my += cval * y;
877 mx2 += cval * x*x;
878 my2 += cval * y*y;
879 mxy += cval * x*y;
880 }
881
882 /*----- compute object's properties */
883 xm = mx / rv; /* mean x */
884 ym = my / rv; /* mean y */
885
886 /*-- In case of blending, use previous barycenters */
887 if ((analyse_type&ANALYSE_ROBUST) && (obj->flag&OBJ_MERGED))
888 {
889 double xn,yn;
890
891 xn = obj->mx-xmin;
892 yn = obj->my-ymin;
893 xm2 = mx2 / rv + xn*xn - 2*xm*xn;
894 ym2 = my2 / rv + yn*yn - 2*ym*yn;
895 xym = mxy / rv + xn*yn - xm*yn - xn*ym;
896 xm = xn;
897 ym = yn;
898 }
899 else
900 {
901 xm2 = mx2 / rv - xm * xm; /* variance of x */
902 ym2 = my2 / rv - ym * ym; /* variance of y */
903 xym = mxy / rv - xm * ym; /* covariance */
904 }
905
906 /* Handle fully correlated x/y (which cause a singularity...) */
907 if ((temp2=xm2*ym2-xym*xym)<0.00694)
908 {
909 xm2 += 0.0833333;
910 ym2 += 0.0833333;
911 temp2 = xm2*ym2-xym*xym;
912 obj->singuflag = 1;
913 }
914 else
915 obj->singuflag = 0;
916
917 if ((fabs(temp=xm2-ym2)) > 0.0)
918 theta = atan2(2.0 * xym,temp) / 2.0;
919 else
920 theta = PI/4.0;
921
922 temp = sqrt(0.25*temp*temp+xym*xym);
923 pmy2 = pmx2 = 0.5*(xm2+ym2);
924 pmx2+=temp;
925 pmy2-=temp;
926
927 obj->dnpix = (obj->flag & OBJ_OVERFLOW)? obj->fdnpix:(LONG)dnpix;
928 obj->dflux = tv;
929 obj->mx = xm+xmin; /* add back xmin */
930 obj->my = ym+ymin; /* add back ymin */
931 obj->mx2 = xm2;
932 obj->my2 = ym2;
933 obj->mxy = xym;
934 obj->a = (float)sqrt(pmx2);
935 obj->b = (float)sqrt(pmy2);
936 obj->theta = theta*180.0/PI;
937
938 obj->cxx = (float)(ym2/temp2);
939 obj->cyy = (float)(xm2/temp2);
940 obj->cxy = (float)(-2*xym/temp2);
941
942 darea = (double)area2 - dnpix;
943 t1t2 = thresh/thresh2;
944 if (t1t2>0.0 && !prefs.dweight_flag)
945 {
946 obj->abcor = (darea<0.0?darea:-1.0)/(2*PI*log(t1t2<1.0?t1t2:0.99)
947 *obj->a*obj->b);
948 if (obj->abcor>1.0)
949 obj->abcor = 1.0;
950 }
951 else
952 obj->abcor = 1.0;
953 }
954
955 return;
956 }
957
958