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