File indexing completed on 2024-04-28 15:10:10
0001 /* 0002 SPDX-FileCopyrightText: 1993-2011 Emmanuel Bertin -- IAP /CNRS/UPMC 0003 SPDX-FileCopyrightText: 2014 SEP developers 0004 0005 SPDX-License-Identifier: LGPL-3.0-or-later 0006 */ 0007 0008 #include <math.h> 0009 #include <stdio.h> 0010 #include <stdlib.h> 0011 #include <string.h> 0012 #include "sep.h" 0013 #include "sepcore.h" 0014 #include "extract.h" 0015 0016 #ifndef RAND_MAX 0017 #define RAND_MAX 2147483647 0018 #endif 0019 #define NSONMAX 1024 /* max. number per level */ 0020 #define NSONMAX_STR "1024" /* just for error message */ 0021 #define NBRANCH 16 /* starting number per branch */ 0022 0023 0024 0025 int belong(int, objliststruct *, int, objliststruct *); 0026 int *createsubmap(objliststruct *, int, int *, int *, int *, int *); 0027 int gatherup(objliststruct *, objliststruct *); 0028 0029 static objliststruct *objlist=NULL; 0030 static short *son=NULL, *ok=NULL; 0031 0032 /******************************** deblend ************************************/ 0033 /* 0034 Divide a list of isophotal detections in several parts (deblending). 0035 NOTE: Even if the object is not deblended, the output objlist threshold is 0036 recomputed if a variable threshold is used. 0037 0038 This can return two error codes: DEBLEND_OVERFLOW or MEMORY_ALLOC_ERROR 0039 */ 0040 int deblend(objliststruct *objlistin, int l, objliststruct *objlistout, 0041 int deblend_nthresh, double deblend_mincont, int minarea) 0042 { 0043 objstruct *obj; 0044 static objliststruct debobjlist, debobjlist2; 0045 double thresh, thresh0, value0; 0046 int h,i,j,k,m,subx,suby,subh,subw, 0047 xn, 0048 nbm = NBRANCH, 0049 status; 0050 int *submap; 0051 0052 submap = NULL; 0053 status = RETURN_OK; 0054 xn = deblend_nthresh; 0055 0056 /* reset global static objlist for deblending */ 0057 memset(objlist, 0, (size_t)xn*sizeof(objliststruct)); 0058 0059 /* initialize local object lists */ 0060 debobjlist.obj = debobjlist2.obj = NULL; 0061 debobjlist.plist = debobjlist2.plist = NULL; 0062 debobjlist.nobj = debobjlist2.nobj = 0; 0063 debobjlist.npix = debobjlist2.npix = 0; 0064 0065 /* Create the submap for the object. 0066 * The submap is used in lutz(). We create it here because we may call 0067 * lutz multiple times below, and we only want to create it once. 0068 */ 0069 submap = createsubmap(objlistin, l, &subx, &suby, &subw, &subh); 0070 if (!submap) 0071 { 0072 status = MEMORY_ALLOC_ERROR; 0073 goto exit; 0074 } 0075 0076 /* set thresholds of object lists based on object threshold */ 0077 thresh0 = objlistin->obj[l].thresh; 0078 objlistout->thresh = debobjlist2.thresh = thresh0; 0079 0080 /* add input object to global deblending objlist and one local objlist */ 0081 if ((status = addobjdeep(l, objlistin, &objlist[0])) != RETURN_OK) 0082 goto exit; 0083 if ((status = addobjdeep(l, objlistin, &debobjlist2)) != RETURN_OK) 0084 goto exit; 0085 0086 value0 = objlist[0].obj[0].fdflux*deblend_mincont; 0087 ok[0] = (short)1; 0088 for (k=1; k<xn; k++) 0089 { 0090 /*------ Calculate threshold */ 0091 thresh = objlistin->obj[l].fdpeak; 0092 debobjlist.thresh = thresh > 0.0? 0093 thresh0*pow(thresh/thresh0,(double)k/xn) : thresh0; 0094 0095 /*--------- Build tree (bottom->up) */ 0096 if (objlist[k-1].nobj>=NSONMAX) 0097 { 0098 status = DEBLEND_OVERFLOW; 0099 goto exit; 0100 } 0101 0102 for (i=0; i<objlist[k-1].nobj; i++) 0103 { 0104 status = lutz(objlistin->plist, submap, subx, suby, subw, 0105 &objlist[k-1].obj[i], &debobjlist, minarea); 0106 if (status != RETURN_OK) 0107 goto exit; 0108 0109 for (j=h=0; j<debobjlist.nobj; j++) 0110 if (belong(j, &debobjlist, i, &objlist[k-1])) 0111 { 0112 debobjlist.obj[j].thresh = debobjlist.thresh; 0113 if ((status = addobjdeep(j, &debobjlist, &objlist[k])) 0114 != RETURN_OK) 0115 goto exit; 0116 m = objlist[k].nobj - 1; 0117 if (m>=NSONMAX) 0118 { 0119 status = DEBLEND_OVERFLOW; 0120 goto exit; 0121 } 0122 if (h>=nbm-1) 0123 if (!(son = (short *) 0124 realloc(son,xn*NSONMAX*(nbm+=16)*sizeof(short)))) 0125 { 0126 status = MEMORY_ALLOC_ERROR; 0127 goto exit; 0128 } 0129 son[k-1+xn*(i+NSONMAX*(h++))] = (short)m; 0130 ok[k+xn*m] = (short)1; 0131 } 0132 son[k-1+xn*(i+NSONMAX*h)] = (short)-1; 0133 } 0134 } 0135 0136 /*------- cut the right branches (top->down) */ 0137 for (k = xn-2; k>=0; k--) 0138 { 0139 obj = objlist[k+1].obj; 0140 for (i=0; i<objlist[k].nobj; i++) 0141 { 0142 for (m=h=0; (j=(int)son[k+xn*(i+NSONMAX*h)])!=-1; h++) 0143 { 0144 if (obj[j].fdflux - obj[j].thresh * obj[j].fdnpix > value0) 0145 m++; 0146 ok[k+xn*i] &= ok[k+1+xn*j]; 0147 } 0148 if (m>1) 0149 { 0150 for (h=0; (j=(int)son[k+xn*(i+NSONMAX*h)])!=-1; h++) 0151 if (ok[k+1+xn*j] && 0152 obj[j].fdflux - obj[j].thresh * obj[j].fdnpix > value0) 0153 { 0154 objlist[k+1].obj[j].flag |= SEP_OBJ_MERGED; 0155 status = addobjdeep(j, &objlist[k+1], &debobjlist2); 0156 if (status != RETURN_OK) 0157 goto exit; 0158 } 0159 ok[k+xn*i] = (short)0; 0160 } 0161 } 0162 } 0163 0164 if (ok[0]) 0165 status = addobjdeep(0, &debobjlist2, objlistout); 0166 else 0167 status = gatherup(&debobjlist2, objlistout); 0168 0169 exit: 0170 if (status == DEBLEND_OVERFLOW) 0171 put_errdetail("limit of " NSONMAX_STR " sub-objects reached while " 0172 "deblending. Decrease number of deblending thresholds " 0173 "or increase the detection threshold."); 0174 0175 free(submap); 0176 submap = NULL; 0177 free(debobjlist2.obj); 0178 free(debobjlist2.plist); 0179 0180 for (k=0; k<xn; k++) 0181 { 0182 free(objlist[k].obj); 0183 free(objlist[k].plist); 0184 } 0185 0186 free(debobjlist.obj); 0187 free(debobjlist.plist); 0188 0189 return status; 0190 } 0191 0192 0193 /******************************* allocdeblend ******************************/ 0194 /* 0195 Allocate the memory allocated by global pointers in refine.c 0196 */ 0197 int allocdeblend(int deblend_nthresh) 0198 { 0199 int status=RETURN_OK; 0200 QMALLOC(son, short, deblend_nthresh*NSONMAX*NBRANCH, status); 0201 QMALLOC(ok, short, deblend_nthresh*NSONMAX, status); 0202 QMALLOC(objlist, objliststruct, deblend_nthresh, status); 0203 0204 return status; 0205 exit: 0206 freedeblend(); 0207 return status; 0208 } 0209 0210 /******************************* freedeblend *******************************/ 0211 /* 0212 Free the memory allocated by global pointers in refine.c 0213 */ 0214 void freedeblend(void) 0215 { 0216 free(son); 0217 son = NULL; 0218 free(ok); 0219 ok = NULL; 0220 free(objlist); 0221 objlist = NULL; 0222 return; 0223 } 0224 0225 /********************************* gatherup **********************************/ 0226 /* 0227 Collect faint remaining pixels and allocate them to their most probable 0228 progenitor. 0229 */ 0230 int gatherup(objliststruct *objlistin, objliststruct *objlistout) 0231 { 0232 char *bmp; 0233 float *amp, *p, dx,dy, drand, dist, distmin; 0234 objstruct *objin = objlistin->obj, *objout, *objt; 0235 0236 pliststruct *pixelin = objlistin->plist, *pixelout, *pixt,*pixt2; 0237 0238 int i,k,l, *n, iclst, npix, bmwidth, 0239 nobj = objlistin->nobj, xs,ys, x,y, status; 0240 0241 bmp = NULL; 0242 amp = p = NULL; 0243 n = NULL; 0244 status = RETURN_OK; 0245 0246 objlistout->thresh = objlistin->thresh; 0247 0248 QMALLOC(amp, float, nobj, status); 0249 QMALLOC(p, float, nobj, status); 0250 QMALLOC(n, int, nobj, status); 0251 0252 for (i=1; i<nobj; i++) 0253 analyse(i, objlistin, 0, 0.0); 0254 0255 p[0] = 0.0; 0256 bmwidth = objin->xmax - (xs=objin->xmin) + 1; 0257 npix = bmwidth * (objin->ymax - (ys=objin->ymin) + 1); 0258 if (!(bmp = (char *)calloc(1, npix*sizeof(char)))) 0259 { 0260 bmp = NULL; 0261 status = MEMORY_ALLOC_ERROR; 0262 goto exit; 0263 } 0264 0265 for (objt = objin+(i=1); i<nobj; i++, objt++) 0266 { 0267 /*-- Now we have passed the deblending section, reset threshold */ 0268 objt->thresh = objlistin->thresh; 0269 0270 /* ------------ flag pixels which are already allocated */ 0271 for (pixt=pixelin+objin[i].firstpix; pixt>=pixelin; 0272 pixt=pixelin+PLIST(pixt,nextpix)) 0273 bmp[(PLIST(pixt,x)-xs) + (PLIST(pixt,y)-ys)*bmwidth] = '\1'; 0274 0275 status = addobjdeep(i, objlistin, objlistout); 0276 if (status != RETURN_OK) 0277 goto exit; 0278 n[i] = objlistout->nobj - 1; 0279 0280 dist = objt->fdnpix/(2*PI*objt->abcor*objt->a*objt->b); 0281 amp[i] = dist<70.0? objt->thresh * expf(dist) : 4.0*objt->fdpeak; 0282 0283 /* ------------ limitate expansion ! */ 0284 if (amp[i]>4.0*objt->fdpeak) 0285 amp[i] = 4.0*objt->fdpeak; 0286 } 0287 0288 objout = objlistout->obj; /* DO NOT MOVE !!! */ 0289 0290 if (!(pixelout=(pliststruct *)realloc(objlistout->plist, 0291 (objlistout->npix + npix)*plistsize))) 0292 { 0293 status = MEMORY_ALLOC_ERROR; 0294 goto exit; 0295 } 0296 0297 objlistout->plist = pixelout; 0298 k = objlistout->npix; 0299 iclst = 0; /* To avoid gcc -Wall warnings */ 0300 for (pixt=pixelin+objin->firstpix; pixt>=pixelin; 0301 pixt=pixelin+PLIST(pixt,nextpix)) 0302 { 0303 x = PLIST(pixt,x); 0304 y = PLIST(pixt,y); 0305 if (!bmp[(x-xs) + (y-ys)*bmwidth]) 0306 { 0307 pixt2 = pixelout + (l=(k++*plistsize)); 0308 memcpy(pixt2, pixt, (size_t)plistsize); 0309 PLIST(pixt2, nextpix) = -1; 0310 distmin = 1e+31; 0311 for (objt = objin+(i=1); i<nobj; i++, objt++) 0312 { 0313 dx = x - objt->mx; 0314 dy = y - objt->my; 0315 dist=0.5*(objt->cxx*dx*dx+objt->cyy*dy*dy+objt->cxy*dx*dy)/objt->abcor; 0316 p[i] = p[i-1] + (dist<70.0?amp[i]*expf(-dist) : 0.0); 0317 if (dist<distmin) 0318 { 0319 distmin = dist; 0320 iclst = i; 0321 } 0322 } 0323 if (p[nobj-1] > 1.0e-31) 0324 { 0325 drand = p[nobj-1]*rand()/RAND_MAX; 0326 for (i=1; i<nobj && p[i]<drand; i++); 0327 if (i==nobj) 0328 i=iclst; 0329 } 0330 else 0331 i = iclst; 0332 objout[n[i]].lastpix=PLIST(pixelout+objout[n[i]].lastpix,nextpix)=l; 0333 } 0334 } 0335 0336 objlistout->npix = k; 0337 if (!(objlistout->plist = (pliststruct *)realloc(pixelout, 0338 objlistout->npix*plistsize))) 0339 status = MEMORY_ALLOC_ERROR; 0340 0341 exit: 0342 free(bmp); 0343 free(amp); 0344 free(p); 0345 free(n); 0346 0347 return status; 0348 } 0349 0350 /**************** belong (originally in manobjlist.c) ************************/ 0351 /* 0352 * say if an object is "included" in another. Returns 1 if the pixels of the 0353 * first object are included in the pixels of the second object. 0354 */ 0355 0356 int belong(int corenb, objliststruct *coreobjlist, 0357 int shellnb, objliststruct *shellobjlist) 0358 { 0359 objstruct *cobj = &(coreobjlist->obj[corenb]), 0360 *sobj = &(shellobjlist->obj[shellnb]); 0361 pliststruct *cpl = coreobjlist->plist, *spl = shellobjlist->plist, *pixt; 0362 0363 int xc=PLIST(cpl+cobj->firstpix,x), yc=PLIST(cpl+cobj->firstpix,y); 0364 0365 for (pixt = spl+sobj->firstpix; pixt>=spl; pixt = spl+PLIST(pixt,nextpix)) 0366 if ((PLIST(pixt,x) == xc) && (PLIST(pixt,y) == yc)) 0367 return 1; 0368 0369 return 0; 0370 } 0371 0372 0373 /******************************** createsubmap *******************************/ 0374 /* 0375 Create pixel-index submap for deblending. 0376 */ 0377 int *createsubmap(objliststruct *objlistin, int no, 0378 int *subx, int *suby, int *subw, int *subh) 0379 { 0380 objstruct *obj; 0381 pliststruct *pixel, *pixt; 0382 int i, n, xmin,ymin, w, *pix, *pt, *submap; 0383 0384 obj = objlistin->obj+no; 0385 pixel = objlistin->plist; 0386 0387 *subx = xmin = obj->xmin; 0388 *suby = ymin = obj->ymin; 0389 *subw = w = obj->xmax - xmin + 1; 0390 *subh = obj->ymax - ymin + 1; 0391 0392 n = w**subh; 0393 if (!(submap = pix = (int *)malloc(n*sizeof(int)))) 0394 return NULL; 0395 pt = pix; 0396 for (i=n; i--;) 0397 *(pt++) = -1; 0398 0399 for (i=obj->firstpix; i!=-1; i=PLIST(pixt,nextpix)) 0400 { 0401 pixt = pixel+i; 0402 *(pix+(PLIST(pixt,x)-xmin) + (PLIST(pixt,y)-ymin)*w) = i; 0403 } 0404 0405 return submap; 0406 }