00001
00002
00003
00004
00005
00006
00007 #include <stdio.h>
00008 #include <stdlib.h>
00009 #include <math.h>
00010
00011 #include <cpl.h>
00012
00013 #include "ap.h"
00014 #include "util.h"
00015 #include "imcore.h"
00016 #include "floatmath.h"
00017
00018 #define NITER 6
00019
00020 static void moments_thr(ap_t *, float [], int []);
00021 static void sort_on_zsm_rev(int, plstruct *);
00022 static void update_ov(float [], float, float, float, float);
00023 static void check_term(ap_t *, int *, float [IMNUM][NPAR+1], int [IMNUM][2],
00024 int *);
00025
00026 static float oldthr;
00027 static float curthr;
00028 static float nexthr;
00029 static float lasthr;
00030 static float xbar_start;
00031 static float ybar_start;
00032
00033
00034
00035 extern void overlp(ap_t *ap, float parm[IMNUM][NPAR], int *nbit, float xbar,
00036 float ybar, float total, int npix, float tmax) {
00037 plstruct *pl;
00038 int npl,ipix,ipixo2,npl2,nbitprev,nobj,toomany,i,isnew,k,kk,j,ibitx[IMNUM];
00039 int ibity[IMNUM],ibitl[IMNUM],iwas,iupdate[IMNUM],npl3,iter,lastone,ic,jj;
00040 int ii,conv,ipks[IMNUM][2];
00041 float fconst,offset,tmul,smul,xintmn,itmaxlim,algthr,radmax,xb,yb,radius2;
00042 float results[IMNUM][NPAR+1],distmax,dx,dy,parmnew[IMNUM][NPAR],sumint;
00043 float xlevol,radold,slope,xx,xlevel,radius,xdat[NAREAL+1],xcor[NAREAL+1];
00044 float dlbydr,wt,dist,xeff,polycf[3],ttt,radthr,delb,deli,ratio;
00045 float bitx[IMNUM],bitl[IMNUM],sxx,syy;
00046 ap_t ap2;
00047
00048
00049
00050 pl = ap->plarray;
00051 npl = ap->npl_pix;
00052 ipix = ap->ipnop;
00053 oldthr = ap->thresh;
00054 fconst = ap->fconst;
00055 offset = ap->areal_offset;
00056 xbar_start = xbar;
00057 ybar_start = ybar;
00058
00059
00060
00061 tmul = 1.2589678;
00062 smul = 2.5;
00063 ipixo2 = MAX(2,(ipix + 1)/2);
00064 xintmn = oldthr*ipixo2;
00065 itmaxlim = 0.9*tmax;
00066 lasthr = itmaxlim;
00067 algthr = logf(oldthr);
00068 radmax = sqrtf(((float)npix)/M_PI);
00069
00070
00071
00072
00073
00074
00075
00076
00077 curthr = smul*oldthr;
00078 sort_on_zsm_rev(npl,pl);
00079 while (1) {
00080 npl2 = 0;
00081 while (pl[npl2].zsm > curthr && npl2 < npl-1)
00082 npl2++;
00083 if (npl2 > IDBLIM)
00084 curthr += oldthr;
00085 else
00086 break;
00087 }
00088
00089
00090
00091
00092 if (npl2 < ipix) {
00093 *nbit = 1;
00094 return;
00095 }
00096
00097
00098
00099 ap2.lsiz = ap->lsiz;
00100 ap2.csiz = ap->csiz;
00101 ap2.multiply = 1;
00102 ap2.ipnop = ipixo2;
00103 ap2.areal_offset = offset;
00104 ap2.fconst = fconst;
00105 ap2.mflag = ap->mflag;
00106
00107
00108
00109 *nbit = 0;
00110 nbitprev = 0;
00111 apinit(&ap2);
00112 while (1) {
00113 nexthr = MAX(curthr+oldthr,curthr*tmul);
00114
00115
00116
00117 ap2.thresh = curthr;
00118 apclust(&ap2,npl2,pl);
00119 check_term(&ap2,&nobj,results,ipks,&toomany);
00120 apreinit(&ap2);
00121 if (nobj == 0)
00122 break;
00123
00124
00125
00126
00127 for (i = 0; i < nobj; i++) {
00128
00129
00130
00131
00132
00133
00134 isnew = 1;
00135 xb = results[i][1];
00136 yb = results[i][2];
00137 sxx = MAX(1.0,results[i][4]);
00138 syy = MAX(1.0,results[i][6]);
00139 for (k = 0; k < nbitprev; k++) {
00140 dx = xb - parm[k][1];
00141 dy = yb - parm[k][2];
00142 radius2 = dx*dx/sxx + dy*dy/syy;
00143 if ((ibitx[k] == ipks[i][0] && ibity[k] == ipks[i][1]) ||
00144 radius2 < 1.0) {
00145 isnew = 0;
00146 for (kk = 0; kk < NPAR; kk++)
00147 parmnew[k][kk] = results[i][kk];
00148 ibitl[k] = (int)results[i][NPAR];
00149 break;
00150 }
00151 }
00152
00153
00154
00155
00156
00157
00158
00159 if (isnew && results[i][0] > xintmn) {
00160 if (*nbit >= IMNUM) {
00161 *nbit = IMNUM;
00162 toomany = 1;
00163 break;
00164 }
00165 ibitx[*nbit] = ipks[i][0];
00166 ibity[*nbit] = ipks[i][1];
00167 for (kk = 0; kk < NPAR; kk++)
00168 parm[*nbit][kk] = results[i][kk];
00169 ibitl[*nbit] = (int)results[i][NPAR];
00170 (*nbit)++;
00171 }
00172 }
00173
00174
00175
00176
00177
00178
00179 if (! toomany) {
00180 if (*nbit > nbitprev && nbitprev > 0) {
00181 for (i = 0; i < nbitprev; i++)
00182 iupdate[i] = 0;
00183 for (j = nbitprev; j < *nbit; j++) {
00184 distmax = 0.0;
00185 iwas = 0;
00186 for (i = 0; i < nbitprev; i++) {
00187 if (parmnew[i][0] > 0.0) {
00188 dx = parmnew[i][1] - parm[i][1];
00189 dy = parmnew[i][2] - parm[i][2];
00190 radius2 = dx*dx + dy*dy;
00191 if (radius2 > distmax) {
00192 iwas = i;
00193 distmax = radius2;
00194 }
00195 }
00196 }
00197 iupdate[iwas] = 1;
00198 }
00199 for (i = 0; i < nbitprev; i++)
00200 if (iupdate[i] == 1 && parmnew[i][0] > 0.0)
00201 for (j = 0; j < NPAR; j++)
00202 parm[i][j] = parmnew[i][j];
00203 }
00204
00205
00206
00207 for (i = 0; i <= *nbit; i++)
00208 parmnew[i][0] = -1.0;
00209 nbitprev = *nbit;
00210 }
00211
00212
00213
00214
00215 npl3 = 0;
00216 while (pl[npl3].zsm > nexthr && npl3 < npl2-1)
00217 npl3++;
00218 npl2 = npl3;
00219
00220
00221
00222 if (npl2 == 0 || toomany || nexthr >= itmaxlim)
00223 break;
00224
00225
00226
00227 curthr = nexthr;
00228
00229 }
00230
00231
00232
00233 apclose(&ap2);
00234
00235
00236
00237 if (*nbit == 1)
00238 return;
00239
00240
00241 j = -1;
00242 for (k = 0; k < *nbit; k++) {
00243
00244
00245
00246
00247 if (parm[k][0] > xintmn) {
00248 j++;
00249 if (j != k)
00250 for (i = 0; i < NPAR; i++)
00251 parm[j][i] = parm[k][i];
00252 }
00253 }
00254 *nbit = j + 1;
00255 for (j = 0; j < *nbit; j++) {
00256 bitx[j] = 0.0;
00257 bitl[j] = 0.0;
00258 }
00259
00260
00261
00262
00263 iter = 0;
00264 sumint = 0.0;
00265 lastone = 0;
00266 while (iter < NITER) {
00267 iter++;
00268
00269
00270
00271 for (k = 0; k < *nbit; k++) {
00272 if (parm[k][0] < 0.0)
00273 continue;
00274 xlevol = logf(parm[k][7] + parm[k][3] - bitl[k]);
00275 xlevel = xlevol;
00276 radold = 0.0;
00277 radius = radold;
00278 slope = 1.0;
00279 ic = 0;
00280 for (i = 1; i <= NAREAL; i++) {
00281 jj = NPAR - i;
00282 ii = NAREAL - i;
00283 xx = (float)ii + offset;
00284 if (parm[k][jj] > 0.5) {
00285 if (ii == 0)
00286 xlevel = logf(parm[k][3] - bitl[k] + 0.5);
00287 else
00288 xlevel = logf(powf(2.0,xx) - oldthr + parm[k][3] -
00289 bitl[k] - 0.5);
00290 radius = sqrt(parm[k][jj]/M_PI);
00291 xdat[ic] = xlevel;
00292 xcor[ic++] = radius;
00293 dlbydr = (xlevol - xlevel)/MAX(0.01,radius-radold);
00294 wt = MIN(1.0,MAX((radius-radold)*5.0,0.1));
00295 slope = (1.0 - 0.5*wt)*slope + 0.5*wt*MIN(5.0,dlbydr);
00296 radold = radius;
00297 xlevol = xlevel;
00298 }
00299 }
00300
00301
00302
00303
00304 if (! lastone) {
00305 for (i = 0; i < *nbit; i++) {
00306 if (parm[i][0] >= 0.0 && i != k) {
00307 dx = parm[k][1] - parm[i][1];
00308 dy = parm[k][2] - parm[i][2];
00309 dist = sqrtf(dx*dx + dy*dy);
00310 xeff = xlevel - MAX(0.0,MIN(50.0,slope*(dist-radius)));
00311 bitx[i] += expf(xeff);
00312 }
00313 }
00314
00315
00316
00317
00318 } else {
00319 if (ic > 2) {
00320 polynm(xdat,xcor,ic,polycf,3,0);
00321 ttt = polycf[1] + 2.0*polycf[2]*radius;
00322 } else
00323 ttt = 0.0;
00324 slope = MAX(0.1,MAX(-ttt,slope));
00325 radthr = radius + (xlevel - algthr)/slope;
00326 if (radthr > radmax) {
00327 slope = 1.0;
00328 radthr = radmax;
00329 }
00330
00331
00332
00333 delb = parm[k][8]*(parm[k][3] - bitl[k]);
00334 parm[k][8] = M_PI*radthr*radthr;
00335
00336
00337
00338 parm[k][7] += (parm[k][3] - bitl[k]);
00339
00340
00341
00342 deli = 2.0*M_PI*((parm[k][3] - bitl[k])*(1.0 + slope*radius) -
00343 oldthr*(1.0 + slope*radthr))/(slope*slope);
00344 parm[k][0] += delb + MAX(0.0,deli);
00345 for (i = 0; i < 7; i++)
00346 parm[k][i+9] = -1.0;
00347 if (parm[k][0] > xintmn)
00348 sumint += parm[k][0];
00349 }
00350 }
00351
00352
00353
00354
00355
00356 if (! lastone) {
00357 conv = 1;
00358 for (i = 0; i < *nbit; i++) {
00359 if (parm[i][0] >= 0.0) {
00360 if (abs(bitx[i] - bitl[i]) > 3.0)
00361 conv = 0;
00362 bitl[i] = bitx[i];
00363 bitx[i] = 0;
00364 bitl[i] = MIN(bitl[i],NINT(parm[i][3]-oldthr));
00365 }
00366 }
00367 lastone = (conv || (iter == NITER-1));
00368 } else {
00369 break;
00370 }
00371 }
00372
00373
00374
00375 if (sumint == 0.0) {
00376 *nbit = 1;
00377 return;
00378 } else
00379 ratio = total/sumint;
00380 for (i = 0; i < *nbit; i++)
00381 parm[i][0] = ratio*parm[i][0];
00382 }
00383
00384
00385 static void sort_on_zsm_rev(int npts, plstruct *pts) {
00386 int i,j,ii,jj,ifin;
00387 plstruct tmp;
00388
00389 jj = 4;
00390 while (jj < npts)
00391 jj = 2*jj;
00392 jj = MIN(npts,(3*jj)/4-1);
00393 while (jj > 1) {
00394 jj = jj/2;
00395 ifin = npts - jj;
00396 for (ii = 0; ii < ifin; ii++) {
00397 i = ii;
00398 j = i + jj;
00399 if (pts[i].zsm > pts[j].zsm)
00400 continue;
00401 tmp = pts[j];
00402 do {
00403 pts[j] = pts[i];
00404 j = i;
00405 i = i - jj;
00406 if (i < 0)
00407 break;
00408 } while (pts[i].zsm <= tmp.zsm);
00409 pts[j] = tmp;
00410 }
00411 }
00412 }
00413
00414 static void moments_thr(ap_t *ap, float results[NPAR+1], int ipk[2]) {
00415 int i,np,nnext;
00416 float x,y,xoff,yoff,xsum,ysum,xsumsq,ysumsq,tsum,xysum,t,tmax,twelfth;
00417 float xbar,ybar,sxx,syy,sxy,fconst,offset,xsum_w,ysum_w,wsum,w;
00418 plstruct *plarray;
00419
00420
00421
00422 fconst = ap->fconst;
00423 offset = ap->areal_offset;
00424 plarray = ap->plarray;
00425 np = ap->npl_pix;
00426
00427
00428
00429 xoff = xbar_start;
00430 yoff = ybar_start;
00431 xsum = 0.0;
00432 ysum = 0.0;
00433 xsum_w = 0.0;
00434 ysum_w = 0.0;
00435 wsum = 0.0;
00436 xsumsq = 0.0;
00437 ysumsq = 0.0;
00438 tsum = 0.0;
00439 xysum = 0.0;
00440 tmax = plarray[0].z - curthr;
00441 twelfth = 1.0/12.0;
00442 for (i = 8; i < NPAR; i++)
00443 results[i] = 0.0;
00444
00445
00446
00447 nnext = 0;
00448 for (i = 0; i < np; i++) {
00449 x = (float)plarray[i].x - xoff;
00450 y = (float)plarray[i].y - yoff;
00451 t = plarray[i].z - curthr;
00452 w = plarray[i].zsm;
00453 if (w > nexthr)
00454 nnext++;
00455 xsum += t*x;
00456 ysum += t*y;
00457 tsum += t;
00458 xsum_w += w*t*x;
00459 ysum_w += w*t*y;
00460 wsum += w*t;
00461 xsumsq += (x*x + twelfth)*t;
00462 ysumsq += (y*y + twelfth)*t;
00463 xysum += x*y*t;
00464 update_ov(results+8,t,oldthr,fconst,offset);
00465 if (t > tmax) {
00466 ipk[0] = plarray[i].x;
00467 ipk[1] = plarray[i].y;
00468 tmax = t;
00469 }
00470 }
00471
00472
00473
00474
00475 if (tsum > 0.0) {
00476 results[0] = tsum;
00477 } else {
00478 results[0] = -1.0;
00479 tsum = 1.0;
00480 }
00481 xbar = xsum/tsum;
00482 ybar = ysum/tsum;
00483 sxx = MAX(0.0,(xsumsq/tsum-xbar*xbar));
00484 syy = MAX(0.0,(ysumsq/tsum-ybar*ybar));
00485 sxy = xysum/tsum - xbar*ybar;
00486 wsum = MAX(1.0,wsum);
00487 xbar = xsum_w/wsum;
00488 ybar = ysum_w/wsum;
00489 xbar += xoff;
00490 ybar += yoff;
00491 xbar = MAX(1.0,MIN(xbar,ap->lsiz));
00492 ybar = MAX(1.0,MIN(ybar,ap->csiz));
00493
00494
00495
00496 results[1] = xbar;
00497 results[2] = ybar;
00498 results[3] = curthr;
00499 results[4] = sxx;
00500 results[5] = sxy;
00501 results[6] = syy;
00502 results[7] = tmax;
00503 results[NPAR] = ((nnext > ap->ipnop && nexthr < lasthr) ? 0 : 1);
00504 }
00505
00506 static void update_ov(float iap[NAREAL], float t, float thresh, float fconst,
00507 float offset) {
00508 int nup,i;
00509
00510
00511
00512 if (t <= 0.0)
00513 return;
00514
00515
00516
00517 nup = MAX(1,MIN(NAREAL,(int)(logf(t+thresh)*fconst-offset)+1));
00518 for (i = 0; i < nup; i++)
00519 iap[i] += 1.0;
00520 }
00521
00522 static void check_term(ap_t *ap, int *nobj, float parm[IMNUM][NPAR+1],
00523 int peaks[IMNUM][2], int *toomany) {
00524 int ip,i,ipks[2];
00525 float momresults[NPAR+1];
00526
00527
00528
00529 *nobj = 0;
00530 *toomany = 0;
00531 for (ip = 1; ip <= ap->maxip; ip++) {
00532 if (ap->parent[ip].pnop != -1) {
00533
00534
00535
00536
00537 if ((ap->parent[ip].pnop >= ap->ipnop &&
00538 ap->parent[ip].touch == 0)) {
00539 extract_data(ap,ip);
00540 moments_thr(ap,momresults,ipks);
00541 if (momresults[0] > 0.0) {
00542 if (*nobj == IMNUM-1) {
00543 *toomany = 1;
00544 break;
00545 }
00546 for (i = 0; i <= NPAR; i++)
00547 parm[*nobj][i] = momresults[i];
00548 for (i = 0; i < 2; i++)
00549 peaks[*nobj][i] = ipks[i];
00550 (*nobj)++;
00551 }
00552 }
00553 restack(ap,ip);
00554
00555
00556
00557
00558
00559
00560 }
00561 }
00562 }
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583