1 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2
3 This file is part of SEP
4
5 SPDX-FileCopyrightText: 1993-2011 Emmanuel Bertin -- IAP /CNRS/UPMC
6 SPDX-FileCopyrightText: 2014 SEP developers
7
8 SPDX-License-Identifier: LGPL-3.0-or-later
9 */
10
11 /* Note: was scan.c in SExtractor. */
12
13 #include <math.h>
14 #include <stdio.h>
15 #include <stdlib.h>
16 #include <string.h>
17 #include "sep.h"
18 #include "sepcore.h"
19 #include "extract.h"
20
21 #define DETECT_MAXAREA 0 /* replaces prefs.ext_maxarea */
22 #define WTHRESH_CONVFAC 1e-4 /* Factor to apply to weights when */
23 /* thresholding filtered weight-maps */
24
25 /* globals */
26 int plistexist_cdvalue, plistexist_thresh, plistexist_var;
27 int plistoff_value, plistoff_cdvalue, plistoff_thresh, plistoff_var;
28 int plistsize;
29 size_t extract_pixstack = 1000000;
30
31 /* get and set pixstack */
sep_set_extract_pixstack(size_t val)32 void sep_set_extract_pixstack(size_t val)
33 {
34 extract_pixstack = val;
35 }
36
sep_get_extract_pixstack()37 size_t sep_get_extract_pixstack()
38 {
39 return extract_pixstack;
40 }
41
42 int sortit(infostruct *info, objliststruct *objlist, int minarea,
43 objliststruct *finalobjlist,
44 int deblend_nthresh, double deblend_mincont, double gain);
45 void plistinit(int hasconv, int hasvar);
46 void clean(objliststruct *objlist, double clean_param, int *survives);
47 int convert_to_catalog(objliststruct *objlist, int *survives,
48 sep_catalog *cat, int w, int include_pixels);
49
50 int arraybuffer_init(arraybuffer *buf, void *arr, int dtype, int w, int h,
51 int bufw, int bufh);
52 void arraybuffer_readline(arraybuffer *buf);
53 void arraybuffer_free(arraybuffer *buf);
54
55 /********************* array buffer functions ********************************/
56
57 /* initialize buffer */
58 /* bufw must be less than or equal to w */
arraybuffer_init(arraybuffer * buf,void * arr,int dtype,int w,int h,int bufw,int bufh)59 int arraybuffer_init(arraybuffer *buf, void *arr, int dtype, int w, int h,
60 int bufw, int bufh)
61 {
62 int status, yl;
63 status = RETURN_OK;
64
65 /* data info */
66 buf->dptr = arr;
67 buf->dw = w;
68 buf->dh = h;
69
70 /* buffer array info */
71 buf->bptr = NULL;
72 QMALLOC(buf->bptr, PIXTYPE, bufw*bufh, status);
73 buf->bw = bufw;
74 buf->bh = bufh;
75
76 /* pointers to within buffer */
77 buf->midline = buf->bptr + bufw*(bufh/2); /* ptr to middle buffer line */
78 buf->lastline = buf->bptr + bufw*(bufh-1); /* ptr to last buffer line */
79
80 status = get_array_converter(dtype, &(buf->readline), &(buf->elsize));
81 if (status != RETURN_OK)
82 goto exit;
83
84 /* initialize yoff */
85 buf->yoff = -bufh;
86
87 /* read in lines until the first data line is one line above midline */
88 for (yl=0; yl < bufh - bufh/2 - 1; yl++)
89 arraybuffer_readline(buf);
90
91 return status;
92
93 exit:
94 free(buf->bptr);
95 buf->bptr = NULL;
96 return status;
97 }
98
99 /* read a line into the buffer at the top, shifting all lines down one */
arraybuffer_readline(arraybuffer * buf)100 void arraybuffer_readline(arraybuffer *buf)
101 {
102 PIXTYPE *line;
103 int y;
104
105 /* shift all lines down one */
106 for (line = buf->bptr; line < buf->lastline; line += buf->bw)
107 memcpy(line, line + buf->bw, sizeof(PIXTYPE) * buf->bw);
108
109 /* which image line now corresponds to the last line in buffer? */
110 buf->yoff++;
111 y = buf->yoff + buf->bh - 1;
112
113 if (y < buf->dh)
114 buf->readline(buf->dptr + buf->elsize * buf->dw * y, buf->dw,
115 buf->lastline);
116
117 return;
118 }
119
arraybuffer_free(arraybuffer * buf)120 void arraybuffer_free(arraybuffer *buf)
121 {
122 free(buf->bptr);
123 buf->bptr = NULL;
124 }
125
126 /* apply_mask_line: Apply the mask to the image and noise buffers.
127 *
128 * If convolution is off, masked values should simply be not
129 * detected. For this, would be sufficient to either set data to zero or
130 * set noise (if present) to infinity.
131 *
132 * If convolution is on, strictly speaking, a masked (unknown) pixel
133 * should "poison" the convolved value whenever it is present in the
134 * convolution kernel (e.g., NaN behavior). However, practically we'd
135 * rather use a "best guess" for the value. Without doing
136 * interpolation from neighbors, 0 is the best guess (assuming image
137 * is background subtracted).
138 *
139 * For the purpose of the full matched filter, we should set noise = infinity.
140 *
141 * So, this routine sets masked pixels to zero in the image buffer and
142 * infinity in the noise buffer (if present). It affects the first
143 */
apply_mask_line(arraybuffer * mbuf,arraybuffer * imbuf,arraybuffer * nbuf)144 void apply_mask_line(arraybuffer *mbuf, arraybuffer *imbuf, arraybuffer *nbuf)
145 {
146 int i;
147
148 for (i=0; i<mbuf->bw; i++)
149 {
150 if (mbuf->lastline[i] > 0.0)
151 {
152 imbuf->lastline[i] = 0.0;
153 if (nbuf)
154 nbuf->lastline[i] = BIG;
155 }
156 }
157 }
158
159 /****************************** extract **************************************/
sep_extract(sep_image * image,float thresh,int thresh_type,int minarea,float * conv,int convw,int convh,int filter_type,int deblend_nthresh,double deblend_cont,int clean_flag,double clean_param,sep_catalog ** catalog)160 int sep_extract(sep_image *image, float thresh, int thresh_type,
161 int minarea, float *conv, int convw, int convh,
162 int filter_type, int deblend_nthresh, double deblend_cont,
163 int clean_flag, double clean_param,
164 sep_catalog **catalog)
165 {
166 arraybuffer dbuf, nbuf, mbuf;
167 infostruct curpixinfo, initinfo, freeinfo;
168 objliststruct objlist;
169 char newmarker;
170 size_t mem_pixstack;
171 int nposize, oldnposize;
172 int w, h;
173 int co, i, luflag, pstop, xl, xl2, yl, cn;
174 int stacksize, convn, status;
175 int bufh;
176 int isvarthresh, isvarnoise;
177 short trunflag;
178 PIXTYPE relthresh, cdnewsymbol, pixvar, pixsig;
179 float sum;
180 pixstatus cs, ps;
181
182 infostruct *info, *store;
183 objliststruct *finalobjlist;
184 pliststruct *pixel, *pixt;
185 char *marker;
186 PIXTYPE *scan, *cdscan, *wscan, *dummyscan;
187 PIXTYPE *sigscan, *workscan;
188 float *convnorm;
189 int *start, *end, *survives;
190 pixstatus *psstack;
191 char errtext[512];
192 sep_catalog *cat;
193
194 status = RETURN_OK;
195 pixel = NULL;
196 convnorm = NULL;
197 scan = wscan = cdscan = dummyscan = NULL;
198 sigscan = workscan = NULL;
199 info = NULL;
200 store = NULL;
201 marker = NULL;
202 psstack = NULL;
203 start = end = NULL;
204 finalobjlist = NULL;
205 survives = NULL;
206 cat = NULL;
207 convn = 0;
208 sum = 0.0;
209 w = image->w;
210 h = image->h;
211 isvarthresh = 0;
212 relthresh = 0.0;
213 pixvar = 0.0;
214 pixsig = 0.0;
215 isvarnoise = 0;
216
217 mem_pixstack = sep_get_extract_pixstack();
218
219 /* seed the random number generator consistently on each call to get
220 * consistent results. rand() is used in deblending. */
221 srand(1);
222
223 /* Noise characteristics of the image: None, scalar or variable? */
224 if (image->noise_type == SEP_NOISE_NONE) { } /* nothing to do */
225 else if (image->noise == NULL) {
226 /* noise is constant; we can set pixel noise now. */
227 if (image->noise_type == SEP_NOISE_STDDEV) {
228 pixsig = image->noiseval;
229 pixvar = pixsig * pixsig;
230 }
231 else if (image->noise_type == SEP_NOISE_VAR) {
232 pixvar = image->noiseval;
233 pixsig = sqrt(pixvar);
234 }
235 else {
236 return UNKNOWN_NOISE_TYPE;
237 }
238 }
239 else {
240 /* noise is variable; we deal with setting pixvar and pixsig at each
241 * pixel. */
242 isvarnoise = 1;
243 }
244
245 /* Deal with relative thresholding. (For an absolute threshold
246 * nothing needs to be done, as `thresh` should already contain the constant
247 * threshold, and `isvarthresh` is already 0.) */
248 if (thresh_type == SEP_THRESH_REL) {
249
250 /* The image must have noise information. */
251 if (image->noise_type == SEP_NOISE_NONE) return RELTHRESH_NO_NOISE;
252
253 isvarthresh = isvarnoise; /* threshold is variable if noise is */
254 if (isvarthresh) {
255 relthresh = thresh; /* used below to set `thresh` for each pixel. */
256 }
257 else {
258 /* thresh is constant; convert relative threshold to absolute */
259 thresh *= pixsig;
260 }
261 }
262
263 /* this is input `thresh` regardless of thresh_type. */
264 objlist.thresh = thresh;
265
266 /*Allocate memory for buffers */
267 stacksize = w+1;
268 QMALLOC(info, infostruct, stacksize, status);
269 QCALLOC(store, infostruct, stacksize, status);
270 QMALLOC(marker, char, stacksize, status);
271 QMALLOC(dummyscan, PIXTYPE, stacksize, status);
272 QMALLOC(psstack, pixstatus, stacksize, status);
273 QCALLOC(start, int, stacksize, status);
274 QMALLOC(end, int, stacksize, status);
275 if ((status = lutzalloc(w, h)) != RETURN_OK)
276 goto exit;
277 if ((status = allocdeblend(deblend_nthresh)) != RETURN_OK)
278 goto exit;
279
280 /* Initialize buffers for input array(s).
281 * The buffer size depends on whether or not convolution is active.
282 * If not convolving, the buffer size is just a single line. If convolving,
283 * the buffer height equals the height of the convolution kernel.
284 */
285 bufh = conv ? convh : 1;
286 status = arraybuffer_init(&dbuf, image->data, image->dtype, w, h, stacksize,
287 bufh);
288 if (status != RETURN_OK) goto exit;
289 if (isvarnoise) {
290 status = arraybuffer_init(&nbuf, image->noise, image->ndtype, w, h,
291 stacksize, bufh);
292 if (status != RETURN_OK) goto exit;
293 }
294 if (image->mask) {
295 status = arraybuffer_init(&mbuf, image->mask, image->mdtype, w, h,
296 stacksize, bufh);
297 if (status != RETURN_OK) goto exit;
298 }
299
300 /* `scan` (or `wscan`) is always a pointer to the current line being
301 * processed. It might be the only line in the buffer, or it might be the
302 * middle line. */
303 scan = dbuf.midline;
304 if (isvarnoise)
305 wscan = nbuf.midline;
306
307 /* More initializations */
308 initinfo.pixnb = 0;
309 initinfo.flag = 0;
310 initinfo.firstpix = initinfo.lastpix = -1;
311
312 for (xl=0; xl<stacksize; xl++)
313 {
314 marker[xl] = 0;
315 dummyscan[xl] = -BIG;
316 }
317
318 co = pstop = 0;
319 objlist.nobj = 1;
320 curpixinfo.pixnb = 1;
321
322 /* Init finalobjlist */
323 QMALLOC(finalobjlist, objliststruct, 1, status);
324 finalobjlist->obj = NULL;
325 finalobjlist->plist = NULL;
326 finalobjlist->nobj = finalobjlist->npix = 0;
327
328
329 /* Allocate memory for the pixel list */
330 plistinit((conv != NULL), (image->noise_type != SEP_NOISE_NONE));
331 if (!(pixel = objlist.plist = malloc(nposize=mem_pixstack*plistsize)))
332 {
333 status = MEMORY_ALLOC_ERROR;
334 goto exit;
335 }
336
337 /*----- at the beginning, "free" object fills the whole pixel list */
338 freeinfo.firstpix = 0;
339 freeinfo.lastpix = nposize-plistsize;
340 pixt = pixel;
341 for (i=plistsize; i<nposize; i += plistsize, pixt += plistsize)
342 PLIST(pixt, nextpix) = i;
343 PLIST(pixt, nextpix) = -1;
344
345 /* can only use a matched filter when convolving and when there is a noise
346 * array */
347 if (!(conv && isvarnoise))
348 filter_type = SEP_FILTER_CONV;
349
350 if (conv)
351 {
352 /* allocate memory for convolved buffers */
353 QMALLOC(cdscan, PIXTYPE, stacksize, status);
354 if (filter_type == SEP_FILTER_MATCHED)
355 {
356 QMALLOC(sigscan, PIXTYPE, stacksize, status);
357 QMALLOC(workscan, PIXTYPE, stacksize, status);
358 }
359
360 /* normalize the filter */
361 convn = convw * convh;
362 QMALLOC(convnorm, PIXTYPE, convn, status);
363 for (i=0; i<convn; i++)
364 sum += fabs(conv[i]);
365 for (i=0; i<convn; i++)
366 convnorm[i] = conv[i] / sum;
367 }
368
369 /*----- MAIN LOOP ------ */
370 for (yl=0; yl<=h; yl++)
371 {
372
373 ps = COMPLETE;
374 cs = NONOBJECT;
375
376 /* Need an empty line for Lutz' algorithm to end gracely */
377 if (yl==h)
378 {
379 if (conv)
380 {
381 free(cdscan);
382 cdscan = NULL;
383 if (filter_type == SEP_FILTER_MATCHED)
384 {
385 for (xl=0; xl<stacksize; xl++)
386 {
387 sigscan[xl] = -BIG;
388 }
389 }
390 }
391 cdscan = dummyscan;
392 }
393
394 else
395 {
396 arraybuffer_readline(&dbuf);
397 if (isvarnoise)
398 arraybuffer_readline(&nbuf);
399 if (image->mask)
400 {
401 arraybuffer_readline(&mbuf);
402 apply_mask_line(&mbuf, &dbuf, (isvarnoise? &nbuf: NULL));
403 }
404
405 /* filter the lines */
406 if (conv)
407 {
408 status = convolve(&dbuf, yl, convnorm, convw, convh, cdscan);
409 if (status != RETURN_OK)
410 goto exit;
411
412 if (filter_type == SEP_FILTER_MATCHED)
413 {
414 status = matched_filter(&dbuf, &nbuf, yl, convnorm, convw,
415 convh, workscan, sigscan,
416 image->noise_type);
417
418 if (status != RETURN_OK)
419 goto exit;
420 }
421 }
422 else
423 {
424 cdscan = scan;
425 }
426 }
427
428 trunflag = (yl==0 || yl==h-1)? SEP_OBJ_TRUNC: 0;
429
430 for (xl=0; xl<=w; xl++)
431 {
432 if (xl == w)
433 cdnewsymbol = -BIG;
434 else
435 cdnewsymbol = cdscan[xl];
436
437 newmarker = marker[xl]; /* marker at this pixel */
438 marker[xl] = 0;
439
440 curpixinfo.flag = trunflag;
441
442 /* set pixel variance/noise based on noise array */
443 if (isvarthresh) {
444 if (xl == w || yl == h) {
445 pixsig = pixvar = 0.0;
446 }
447 else if (image->noise_type == SEP_NOISE_VAR) {
448 pixvar = wscan[xl];
449 pixsig = sqrt(pixvar);
450 }
451 else if (image->noise_type == SEP_NOISE_STDDEV) {
452 pixsig = wscan[xl];
453 pixvar = pixsig * pixsig;
454 }
455 else {
456 status = UNKNOWN_NOISE_TYPE;
457 goto exit;
458 }
459
460 /* set `thresh` (This is needed later, even
461 * if filter_type is SEP_FILTER_MATCHED */
462 thresh = relthresh * pixsig;
463 }
464
465 /* luflag: is pixel above thresh (Y/N)? */
466 if (filter_type == SEP_FILTER_MATCHED)
467 luflag = ((xl != w) && (sigscan[xl] > relthresh))? 1: 0;
468 else
469 luflag = cdnewsymbol > thresh? 1: 0;
470
471 if (luflag)
472 {
473 /* flag the current object if we're near the image bounds */
474 if (xl==0 || xl==w-1)
475 curpixinfo.flag |= SEP_OBJ_TRUNC;
476
477 /* point pixt to first free pixel in pixel list */
478 /* and increment the "first free pixel" */
479 pixt = pixel + (cn=freeinfo.firstpix);
480 freeinfo.firstpix = PLIST(pixt, nextpix);
481 curpixinfo.lastpix = curpixinfo.firstpix = cn;
482
483 /* set values for the new pixel */
484 PLIST(pixt, nextpix) = -1;
485 PLIST(pixt, x) = xl;
486 PLIST(pixt, y) = yl;
487 PLIST(pixt, value) = scan[xl];
488 if (PLISTEXIST(cdvalue))
489 PLISTPIX(pixt, cdvalue) = cdnewsymbol;
490 if (PLISTEXIST(var))
491 PLISTPIX(pixt, var) = pixvar;
492 if (PLISTEXIST(thresh))
493 PLISTPIX(pixt, thresh) = thresh;
494
495 /* Check if we have run out of free pixels in objlist.plist */
496 if (freeinfo.firstpix==freeinfo.lastpix)
497 {
498 status = PIXSTACK_FULL;
499 sprintf(errtext,
500 "The limit of %d active object pixels over the "
501 "detection threshold was reached. Check that "
502 "the image is background subtracted and the "
503 "detection threshold is not too low. If you "
504 "need to increase the limit, use "
505 "set_extract_pixstack.",
506 (int)mem_pixstack);
507 put_errdetail(errtext);
508 goto exit;
509
510 /* The code in the rest of this block increases the
511 * stack size as needed. Currently, it is never
512 * executed. This is because it isn't clear that
513 * this is a good idea: most times when the stack
514 * overflows it is due to user error: too-low
515 * threshold or image not background subtracted. */
516
517 /* increase the stack size */
518 oldnposize = nposize;
519 mem_pixstack = (int)(mem_pixstack * 2);
520 nposize = mem_pixstack * plistsize;
521 pixel = (pliststruct *)realloc(pixel, nposize);
522 objlist.plist = pixel;
523 if (!pixel)
524 {
525 status = MEMORY_ALLOC_ERROR;
526 goto exit;
527 }
528
529 /* set next free pixel to the start of the new block
530 * and link up all the pixels in the new block */
531 PLIST(pixel+freeinfo.firstpix, nextpix) = oldnposize;
532 pixt = pixel + oldnposize;
533 for (i=oldnposize + plistsize; i<nposize;
534 i += plistsize, pixt += plistsize)
535 PLIST(pixt, nextpix) = i;
536 PLIST(pixt, nextpix) = -1;
537
538 /* last free pixel is now at the end of the new block */
539 freeinfo.lastpix = nposize - plistsize;
540 }
541 /*------------------------------------------------------------*/
542
543 /* if the current status on this line is not already OBJECT... */
544 /* start segment */
545 if (cs != OBJECT)
546 {
547 cs = OBJECT;
548 if (ps == OBJECT)
549 {
550 if (start[co] == UNKNOWN)
551 {
552 marker[xl] = 'S';
553 start[co] = xl;
554 }
555 else
556 marker[xl] = 's';
557 }
558 else
559 {
560 psstack[pstop++] = ps;
561 marker[xl] = 'S';
562 start[++co] = xl;
563 ps = COMPLETE;
564 info[co] = initinfo;
565 }
566 }
567
568 } /* closes if pixel above threshold */
569
570 /* process new marker ---------------------------------------------*/
571 /* newmarker is marker[ ] at this pixel position before we got to
572 it. We'll only enter this if marker[ ] was set on a previous
573 loop iteration. */
574 if (newmarker)
575 {
576 if (newmarker == 'S')
577 {
578 psstack[pstop++] = ps;
579 if (cs == NONOBJECT)
580 {
581 psstack[pstop++] = COMPLETE;
582 info[++co] = store[xl];
583 start[co] = UNKNOWN;
584 }
585 else
586 update(&info[co], &store[xl], pixel);
587 ps = OBJECT;
588 }
589
590 else if (newmarker == 's')
591 {
592 if ((cs == OBJECT) && (ps == COMPLETE))
593 {
594 pstop--;
595 xl2 = start[co];
596 update (&info[co-1],&info[co], pixel);
597 if (start[--co] == UNKNOWN)
598 start[co] = xl2;
599 else
600 marker[xl2] = 's';
601 }
602 ps = OBJECT;
603 }
604
605 else if (newmarker == 'f')
606 ps = INCOMPLETE;
607
608 else if (newmarker == 'F')
609 {
610 ps = psstack[--pstop];
611 if ((cs == NONOBJECT) && (ps == COMPLETE))
612 {
613 if (start[co] == UNKNOWN)
614 {
615 if ((int)info[co].pixnb >= minarea)
616 {
617 /* update threshold before object is processed */
618 objlist.thresh = thresh;
619
620 status = sortit(&info[co], &objlist, minarea,
621 finalobjlist,
622 deblend_nthresh,deblend_cont,
623 image->gain);
624 if (status != RETURN_OK)
625 goto exit;
626 }
627
628 /* free the chain-list */
629 PLIST(pixel+info[co].lastpix, nextpix) =
630 freeinfo.firstpix;
631 freeinfo.firstpix = info[co].firstpix;
632 }
633 else
634 {
635 marker[end[co]] = 'F';
636 store[start[co]] = info[co];
637 }
638 co--;
639 ps = psstack[--pstop];
640 }
641 }
642 }
643 /* end of if (newmarker) ------------------------------------------*/
644
645 /* update the info or end segment */
646 if (luflag)
647 {
648 update(&info[co], &curpixinfo, pixel);
649 }
650 else if (cs == OBJECT)
651 {
652 cs = NONOBJECT;
653 if (ps != COMPLETE)
654 {
655 marker[xl] = 'f';
656 end[co] = xl;
657 }
658 else
659 {
660 ps = psstack[--pstop];
661 marker[xl] = 'F';
662 store[start[co]] = info[co];
663 co--;
664 }
665 }
666
667 } /*------------ End of the loop over the x's -----------------------*/
668
669 } /*---------------- End of the loop over the y's -----------------------*/
670
671 /* convert `finalobjlist` to an array of `sepobj` structs */
672 /* if cleaning, see which objects "survive" cleaning. */
673 if (clean_flag)
674 {
675 /* Calculate mthresh for all objects in the list (needed for cleaning) */
676 for (i=0; i<finalobjlist->nobj; i++)
677 {
678 status = analysemthresh(i, finalobjlist, minarea, thresh);
679 if (status != RETURN_OK)
680 goto exit;
681 }
682
683 QMALLOC(survives, int, finalobjlist->nobj, status);
684 clean(finalobjlist, clean_param, survives);
685 }
686
687 /* convert to output catalog */
688 QCALLOC(cat, sep_catalog, 1, status);
689 status = convert_to_catalog(finalobjlist, survives, cat, w, 1);
690 if (status != RETURN_OK) goto exit;
691
692 exit:
693 free(finalobjlist->obj);
694 free(finalobjlist->plist);
695 free(finalobjlist);
696 freedeblend();
697 free(pixel);
698 lutzfree();
699 free(info);
700 free(store);
701 free(marker);
702 free(dummyscan);
703 free(psstack);
704 free(start);
705 free(end);
706 free(survives);
707 arraybuffer_free(&dbuf);
708 if (image->noise)
709 arraybuffer_free(&nbuf);
710 if (image->mask)
711 arraybuffer_free(&mbuf);
712 if (conv)
713 free(convnorm);
714 if (filter_type == SEP_FILTER_MATCHED)
715 {
716 free(sigscan);
717 free(workscan);
718 }
719
720 if (status != RETURN_OK)
721 {
722 /* free cdscan if we didn't do it on the last `yl` line */
723 if ((cdscan != NULL) && (cdscan != dummyscan))
724 free(cdscan);
725 /* clean up catalog if it was allocated */
726 sep_catalog_free(cat);
727 cat = NULL;
728 }
729
730 *catalog = cat;
731 return status;
732 }
733
734
735 /********************************* sortit ************************************/
736 /*
737 build the object structure.
738 */
sortit(infostruct * info,objliststruct * objlist,int minarea,objliststruct * finalobjlist,int deblend_nthresh,double deblend_mincont,double gain)739 int sortit(infostruct *info, objliststruct *objlist, int minarea,
740 objliststruct *finalobjlist,
741 int deblend_nthresh, double deblend_mincont, double gain)
742 {
743 objliststruct objlistout, *objlist2;
744 static objstruct obj;
745 int i, status;
746
747 status=RETURN_OK;
748 objlistout.obj = NULL;
749 objlistout.plist = NULL;
750 objlistout.nobj = objlistout.npix = 0;
751
752 /*----- Allocate memory to store object data */
753 objlist->obj = &obj;
754 objlist->nobj = 1;
755
756 memset(&obj, 0, (size_t)sizeof(objstruct));
757 objlist->npix = info->pixnb;
758 obj.firstpix = info->firstpix;
759 obj.lastpix = info->lastpix;
760 obj.flag = info->flag;
761 obj.thresh = objlist->thresh;
762
763 preanalyse(0, objlist);
764
765 status = deblend(objlist, 0, &objlistout, deblend_nthresh, deblend_mincont,
766 minarea);
767 if (status)
768 {
769 /* formerly, this wasn't a fatal error, so a flag was set for
770 * the object and we continued. I'm leaving the flag-setting
771 * here in case we want to change this to a non-fatal error in
772 * the future, but currently the flag setting is irrelevant. */
773 objlist2 = objlist;
774 for (i=0; i<objlist2->nobj; i++)
775 objlist2->obj[i].flag |= SEP_OBJ_DOVERFLOW;
776 goto exit;
777 }
778 else
779 objlist2 = &objlistout;
780
781 /* Analyze the deblended objects and add to the final list */
782 for (i=0; i<objlist2->nobj; i++)
783 {
784 analyse(i, objlist2, 1, gain);
785
786 /* this does nothing if DETECT_MAXAREA is 0 (and it currently is) */
787 if (DETECT_MAXAREA && objlist2->obj[i].fdnpix > DETECT_MAXAREA)
788 continue;
789
790 /* add the object to the final list */
791 status = addobjdeep(i, objlist2, finalobjlist);
792 if (status != RETURN_OK)
793 goto exit;
794 }
795
796 exit:
797 free(objlistout.plist);
798 free(objlistout.obj);
799 return status;
800 }
801
802
803 /********** addobjdeep (originally in manobjlist.c) **************************/
804 /*
805 Add object number `objnb` from list `objl1` to list `objl2`.
806 Unlike `addobjshallow` this also copies plist pixels to the second list.
807 */
808
addobjdeep(int objnb,objliststruct * objl1,objliststruct * objl2)809 int addobjdeep(int objnb, objliststruct *objl1, objliststruct *objl2)
810 {
811 objstruct *objl2obj;
812 pliststruct *plist1 = objl1->plist, *plist2 = objl2->plist;
813 int fp, i, j, npx, objnb2;
814
815 fp = objl2->npix; /* 2nd list's plist size in pixels */
816 j = fp*plistsize; /* 2nd list's plist size in bytes */
817 objnb2 = objl2->nobj; /* # of objects currently in 2nd list*/
818
819 /* Allocate space in `objl2` for the new object */
820 if (objnb2)
821 objl2obj = (objstruct *)realloc(objl2->obj,
822 (++objl2->nobj)*sizeof(objstruct));
823 else
824 objl2obj = (objstruct *)malloc((++objl2->nobj)*sizeof(objstruct));
825
826 if (!objl2obj)
827 goto earlyexit;
828 objl2->obj = objl2obj;
829
830 /* Allocate space for the new object's pixels in 2nd list's plist */
831 npx = objl1->obj[objnb].fdnpix;
832 if (fp)
833 plist2 = (pliststruct *)realloc(plist2, (objl2->npix+=npx)*plistsize);
834 else
835 plist2 = (pliststruct *)malloc((objl2->npix=npx)*plistsize);
836
837 if (!plist2)
838 goto earlyexit;
839 objl2->plist = plist2;
840
841 /* copy the plist */
842 plist2 += j;
843 for(i=objl1->obj[objnb].firstpix; i!=-1; i=PLIST(plist1+i,nextpix))
844 {
845 memcpy(plist2, plist1+i, (size_t)plistsize);
846 PLIST(plist2,nextpix) = (j+=plistsize);
847 plist2 += plistsize;
848 }
849 PLIST(plist2-=plistsize, nextpix) = -1;
850
851 /* copy the object itself */
852 objl2->obj[objnb2] = objl1->obj[objnb];
853 objl2->obj[objnb2].firstpix = fp*plistsize;
854 objl2->obj[objnb2].lastpix = j-plistsize;
855
856 return RETURN_OK;
857
858 /* if early exit, reset 2nd list */
859 earlyexit:
860 objl2->nobj--;
861 objl2->npix = fp;
862 return MEMORY_ALLOC_ERROR;
863 }
864
865
866
867 /****************************** plistinit ************************************
868 * (originally init_plist() in sextractor)
869 PURPOSE initialize a pixel-list and its components.
870 ***/
plistinit(int hasconv,int hasvar)871 void plistinit(int hasconv, int hasvar)
872 {
873 pbliststruct *pbdum = NULL;
874
875 plistsize = sizeof(pbliststruct);
876 plistoff_value = (char *)&pbdum->value - (char *)pbdum;
877
878 if (hasconv)
879 {
880 plistexist_cdvalue = 1;
881 plistoff_cdvalue = plistsize;
882 plistsize += sizeof(PIXTYPE);
883 }
884 else
885 {
886 plistexist_cdvalue = 0;
887 plistoff_cdvalue = plistoff_value;
888 }
889
890 if (hasvar)
891 {
892 plistexist_var = 1;
893 plistoff_var = plistsize;
894 plistsize += sizeof(PIXTYPE);
895
896 plistexist_thresh = 1;
897 plistoff_thresh = plistsize;
898 plistsize += sizeof(PIXTYPE);
899 }
900 else
901 {
902 plistexist_var = 0;
903 plistexist_thresh = 0;
904 }
905
906 return;
907
908 }
909
910
911 /************************** clean an objliststruct ***************************/
912 /*
913 Fill a list with whether each object in the list survived the cleaning
914 (assumes that mthresh has already been calculated for all objects in the list)
915 */
916
clean(objliststruct * objlist,double clean_param,int * survives)917 void clean(objliststruct *objlist, double clean_param, int *survives)
918 {
919 objstruct *obj1, *obj2;
920 int i,j;
921 double amp,ampin,alpha,alphain, unitarea,unitareain,beta,val;
922 float dx,dy,rlim;
923
924 beta = clean_param;
925
926 /* initialize to all surviving */
927 for (i=0; i<objlist->nobj; i++)
928 survives[i] = 1;
929
930 obj1 = objlist->obj;
931 for (i=0; i<objlist->nobj; i++, obj1++)
932 {
933 if (!survives[i])
934 continue;
935
936 /* parameters for test object */
937 unitareain = PI*obj1->a*obj1->b;
938 ampin = obj1->fdflux/(2*unitareain*obj1->abcor);
939 alphain = (pow(ampin/obj1->thresh, 1.0/beta)-1)*unitareain/obj1->fdnpix;
940
941 /* loop over remaining objects in list*/
942 obj2 = obj1 + 1;
943 for (j=i+1; j<objlist->nobj; j++, obj2++)
944 {
945 if (!survives[j])
946 continue;
947
948 dx = obj1->mx - obj2->mx;
949 dy = obj1->my - obj2->my;
950 rlim = obj1->a + obj2->a;
951 rlim *= rlim;
952 if (dx*dx + dy*dy > rlim*CLEAN_ZONE*CLEAN_ZONE)
953 continue;
954
955 /* if obj1 is bigger, see if it eats obj2 */
956 if (obj2->fdflux < obj1->fdflux)
957 {
958 val = 1 + alphain*(obj1->cxx*dx*dx + obj1->cyy*dy*dy +
959 obj1->cxy*dx*dy);
960 if (val>1.0 && ((float)(val<1e10?ampin*pow(val,-beta):0.0) >
961 obj2->mthresh))
962 survives[j] = 0; /* the test object eats this one */
963 }
964
965 /* if obj2 is bigger, see if it eats obj1 */
966 else
967 {
968 unitarea = PI*obj2->a*obj2->b;
969 amp = obj2->fdflux/(2*unitarea*obj2->abcor);
970 alpha = (pow(amp/obj2->thresh, 1.0/beta) - 1) *
971 unitarea/obj2->fdnpix;
972 val = 1 + alpha*(obj2->cxx*dx*dx + obj2->cyy*dy*dy +
973 obj2->cxy*dx*dy);
974 if (val>1.0 && ((float)(val<1e10?amp*pow(val,-beta):0.0) >
975 obj1->mthresh))
976 survives[i] = 0; /* this object eats the test object */
977 }
978
979 } /* inner loop over objlist (obj2) */
980 } /* outer loop of objlist (obj1) */
981 }
982
983
984 /*****************************************************************************/
985 /* sep_catalog manipulations */
986
free_catalog_fields(sep_catalog * catalog)987 void free_catalog_fields(sep_catalog *catalog)
988 {
989 free(catalog->thresh);
990 free(catalog->npix);
991 free(catalog->tnpix);
992 free(catalog->xmin);
993 free(catalog->xmax);
994 free(catalog->ymin);
995 free(catalog->ymax);
996 free(catalog->x);
997 free(catalog->y);
998 free(catalog->x2);
999 free(catalog->y2);
1000 free(catalog->xy);
1001 free(catalog->errx2);
1002 free(catalog->erry2);
1003 free(catalog->errxy);
1004 free(catalog->a);
1005 free(catalog->b);
1006 free(catalog->theta);
1007 free(catalog->cxx);
1008 free(catalog->cyy);
1009 free(catalog->cxy);
1010 free(catalog->cflux);
1011 free(catalog->flux);
1012 free(catalog->cpeak);
1013 free(catalog->peak);
1014 free(catalog->xcpeak);
1015 free(catalog->ycpeak);
1016 free(catalog->xpeak);
1017 free(catalog->ypeak);
1018 free(catalog->flag);
1019
1020 free(catalog->pix);
1021 free(catalog->objectspix);
1022
1023 memset(catalog, 0, sizeof(sep_catalog));
1024 }
1025
1026
1027 /* convert_to_catalog()
1028 *
1029 * Convert the final object list to an output catalog.
1030 *
1031 * `survives`: array of 0 or 1 indicating whether or not to output each object
1032 * (ignored if NULL)
1033 * `cat`: catalog object to be filled.
1034 * `w`: width of image (used to calculate linear indicies).
1035 */
convert_to_catalog(objliststruct * objlist,int * survives,sep_catalog * cat,int w,int include_pixels)1036 int convert_to_catalog(objliststruct *objlist, int *survives,
1037 sep_catalog *cat, int w, int include_pixels)
1038 {
1039 int i, j, k;
1040 int totnpix;
1041 int nobj = 0;
1042 int status = RETURN_OK;
1043 objstruct *obj;
1044 pliststruct *pixt, *pixel;
1045
1046 /* Set struct to zero in case the caller didn't.
1047 * This is important if there is a memory error and we have to call
1048 * free_catalog_fields() to recover at exit */
1049 memset(cat, 0, sizeof(sep_catalog));
1050
1051 /* Count number of surviving objects so that we can allocate the
1052 appropriate amount of space in the output catalog. */
1053 if (survives)
1054 for (i=0; i<objlist->nobj; i++) nobj += survives[i];
1055 else
1056 nobj = objlist->nobj;
1057
1058 /* allocate catalog fields */
1059 cat->nobj = nobj;
1060 QMALLOC(cat->thresh, float, nobj, status);
1061 QMALLOC(cat->npix, int, nobj, status);
1062 QMALLOC(cat->tnpix, int, nobj, status);
1063 QMALLOC(cat->xmin, int, nobj, status);
1064 QMALLOC(cat->xmax, int, nobj, status);
1065 QMALLOC(cat->ymin, int, nobj, status);
1066 QMALLOC(cat->ymax, int, nobj, status);
1067 QMALLOC(cat->x, double, nobj, status);
1068 QMALLOC(cat->y, double, nobj, status);
1069 QMALLOC(cat->x2, double, nobj, status);
1070 QMALLOC(cat->y2, double, nobj, status);
1071 QMALLOC(cat->xy, double, nobj, status);
1072 QMALLOC(cat->errx2, double, nobj, status);
1073 QMALLOC(cat->erry2, double, nobj, status);
1074 QMALLOC(cat->errxy, double, nobj, status);
1075 QMALLOC(cat->a, float, nobj, status);
1076 QMALLOC(cat->b, float, nobj, status);
1077 QMALLOC(cat->theta, float, nobj, status);
1078 QMALLOC(cat->cxx, float, nobj, status);
1079 QMALLOC(cat->cyy, float, nobj, status);
1080 QMALLOC(cat->cxy, float, nobj, status);
1081 QMALLOC(cat->cflux, float, nobj, status);
1082 QMALLOC(cat->flux, float, nobj, status);
1083 QMALLOC(cat->cpeak, float, nobj, status);
1084 QMALLOC(cat->peak, float, nobj, status);
1085 QMALLOC(cat->xcpeak, int, nobj, status);
1086 QMALLOC(cat->ycpeak, int, nobj, status);
1087 QMALLOC(cat->xpeak, int, nobj, status);
1088 QMALLOC(cat->ypeak, int, nobj, status);
1089 QMALLOC(cat->cflux, float, nobj, status);
1090 QMALLOC(cat->flux, float, nobj, status);
1091 QMALLOC(cat->flag, short, nobj, status);
1092
1093 /* fill output arrays */
1094 j = 0; /* running index in output array */
1095 for (i=0; i<objlist->nobj; i++)
1096 {
1097 if ((survives == NULL) || survives[i])
1098 {
1099 obj = objlist->obj + i;
1100 cat->thresh[j] = obj->thresh;
1101 cat->npix[j] = obj->fdnpix;
1102 cat->tnpix[j] = obj->dnpix;
1103 cat->xmin[j] = obj->xmin;
1104 cat->xmax[j] = obj->xmax;
1105 cat->ymin[j] = obj->ymin;
1106 cat->ymax[j] = obj->ymax;
1107 cat->x[j] = obj->mx;
1108 cat->y[j] = obj->my;
1109 cat->x2[j] = obj->mx2;
1110 cat->y2[j] = obj->my2;
1111 cat->xy[j] = obj->mxy;
1112 cat->errx2[j] = obj->errx2;
1113 cat->erry2[j] = obj->erry2;
1114 cat->errxy[j] = obj->errxy;
1115
1116 cat->a[j] = obj->a;
1117 cat->b[j] = obj->b;
1118 cat->theta[j] = obj->theta;
1119
1120 cat->cxx[j] = obj->cxx;
1121 cat->cyy[j] = obj->cyy;
1122 cat->cxy[j] = obj->cxy;
1123
1124 cat->cflux[j] = obj->fdflux; /* these change names */
1125 cat->flux[j] = obj->dflux;
1126 cat->cpeak[j] = obj->fdpeak;
1127 cat->peak[j] = obj->dpeak;
1128
1129 cat->xpeak[j] = obj->xpeak;
1130 cat->ypeak[j] = obj->ypeak;
1131 cat->xcpeak[j] = obj->xcpeak;
1132 cat->ycpeak[j] = obj->ycpeak;
1133
1134 cat->flag[j] = obj->flag;
1135
1136 j++;
1137 }
1138 }
1139
1140 if (include_pixels)
1141 {
1142 /* count the total number of pixels */
1143 totnpix = 0;
1144 for (i=0; i<cat->nobj; i++) totnpix += cat->npix[i];
1145
1146 /* allocate buffer for all objects' pixels */
1147 QMALLOC(cat->objectspix, int, totnpix, status);
1148
1149 /* allocate array of pointers into the above buffer */
1150 QMALLOC(cat->pix, int*, nobj, status);
1151
1152 pixel = objlist->plist;
1153
1154 /* for each object, fill buffer and direct object's to it */
1155 k = 0; /* current position in `objectspix` buffer */
1156 j = 0; /* output object index */
1157 for (i=0; i<objlist->nobj; i++)
1158 {
1159 obj = objlist->obj + i; /* input object */
1160 if ((survives == NULL) || survives[i])
1161 {
1162 /* point this object's pixel list into the buffer. */
1163 cat->pix[j] = cat->objectspix + k;
1164
1165 /* fill the actual pixel values */
1166 for (pixt=pixel+obj->firstpix; pixt>=pixel;
1167 pixt=pixel+PLIST(pixt,nextpix), k++)
1168 {
1169 cat->objectspix[k] = PLIST(pixt,x) + w*PLIST(pixt,y);
1170 }
1171 j++;
1172 }
1173 }
1174 }
1175
1176 exit:
1177 if (status != RETURN_OK)
1178 free_catalog_fields(cat);
1179
1180 return status;
1181 }
1182
1183
sep_catalog_free(sep_catalog * catalog)1184 void sep_catalog_free(sep_catalog *catalog)
1185 {
1186 if (catalog != NULL)
1187 free_catalog_fields(catalog);
1188 free(catalog);
1189 }
1190