public software.sextractor

[/] [trunk/] [src/] [scan.c] - Blame information for rev 219

Go to most recent revision | Details | Compare with Previous | View Log

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