imcore_radii.c

00001 /*
00002 
00003 $Id: imcore_radii.c,v 1.1 2005/09/13 13:25:30 jim Exp $
00004 
00005 */
00006 
00007 #include <stdio.h>
00008 #include <math.h>
00009 
00010 #include "imcore.h"
00011 #include "imcore_radii.h"
00012 #include "floatmath.h"
00013 #include "util.h"
00014 #include "ap.h"
00015 
00016 static float fraction (float x, float y, float r_out);
00017 static void dchole (double a[IMNUM+1][IMNUM+1], double b[IMNUM+1], int n);
00018 
00019 extern float imcore_exprad(float thresh, float peak, float areal0, 
00020                            float rcores[], int naper) {
00021     float pk,r_t,rad;
00022 
00023     /* Work out the radius... */
00024 
00025     pk = MAX(1.5*thresh,peak);
00026     r_t = sqrtf(areal0/M_PI);
00027     rad = 5.0*r_t/logf(pk/thresh);
00028     rad = MAX(r_t,MIN(5.0*r_t,MIN(rad,rcores[naper-1])));
00029     return(rad);
00030 }
00031     
00032 extern float imcore_kronrad(float areal0, float rcores[], float cflux[], 
00033                             int naper) {
00034     int i;
00035     float r_t,rad,wt;
00036 
00037     /* Work out the radius... */
00038 
00039     r_t = sqrtf(areal0/M_PI);
00040     rad = 0.5*rcores[0]*cflux[0];
00041     for (i = 1; i < naper; i++) {
00042         wt = MAX(0.0,cflux[i]-cflux[i-1]);
00043         rad += 0.5*(rcores[i] + rcores[i-1])*wt;
00044     }
00045     rad /= cflux[MIN(naper-1,7)];
00046     rad = MAX(r_t,MIN(5.0*r_t,MIN(2.0*rad,rcores[naper-1])));    
00047     return(rad);
00048 }
00049 
00050 extern float imcore_petrad (float areal0, float rcores[], float cflux[], 
00051                             int naper) {
00052     int j;
00053     float eta,r_t,etaold,r1,r2,r3,r4,r5,r_petr;
00054 
00055     /* Work out petrosian radius */
00056 
00057     r_t = sqrtf(areal0/M_PI);
00058     eta = 1.0;
00059     etaold = eta;
00060     j = 1;
00061     while (eta > 0.2 && j < naper) {
00062         etaold = eta;
00063         r1 = rcores[j]*rcores[j]/(rcores[j-1]*rcores[j-1]) - 1.0;
00064         r2 = cflux[j]/cflux[j-1] - 1.0;
00065         eta = r2/r1;
00066         j++;
00067     }
00068     if (j == naper) {
00069         r_petr = rcores[naper-1];
00070     } else {
00071         r1 = rcores[j]*rcores[j];
00072         r2 = rcores[j-1]*rcores[j-1];
00073         r3 = rcores[j-2]*rcores[j-2];
00074         r4 = (etaold - 0.2)/(etaold - eta);
00075         r5 = (0.2 - eta)/(etaold - eta);
00076         r_petr = r4*sqrt(0.5*(r1 + r2)) + r5*sqrt(0.5*(r2 + r3));
00077     }
00078     r_petr = MAX(r_t,MIN(5.0*r_t,MIN(2.0*r_petr,rcores[naper-1])));
00079     return(r_petr);
00080 }
00081 
00082 extern void imcore_flux(ap_t *ap, float parm[IMNUM][NPAR], int nbit, 
00083                         float apertures[], float cflux[]) {
00084     double aa[IMNUM+1][IMNUM+1],bb[IMNUM+1];
00085     float *map,parrad[IMNUM],cn[IMNUM],xi,yi,d,a2j,a2i,di,dj,argi,argj;
00086     float xmin,xmax,ymin,ymax,t,xj,yj,tj,xk,yk,tk;
00087     unsigned char *mflag;
00088     long nx,ny;
00089     int i,ix1,ix2,iy1,iy2,ii,kk,j,k;
00090 
00091     /* Set up some local variables */
00092  
00093     map = ap->indata;
00094     mflag = ap->mflag;
00095     nx = ap->lsiz;
00096     ny = ap->csiz;
00097 
00098     /* Set up some convenience variables */
00099 
00100     for (i = 0; i < nbit; i++) {
00101         parrad[i] = apertures[i] + 0.5;
00102         cn[i] = 1.0/(M_PI*apertures[i]*apertures[i]);
00103     }
00104  
00105     /* Set up covariance matrix - analytic special case for cores */
00106 
00107     for(i = 0; i < nbit; i++) {
00108         aa[i][i] = cn[i];              /* overlaps totally area=pi*r**2 */
00109         if (nbit > 1) {
00110             xi = parm[i][1];
00111             yi = parm[i][2];
00112             for(j = i+1; j < nbit; j++) {
00113                 d = sqrtf((xi-parm[j][1])*(xi-parm[j][1])
00114                     + (yi-parm[j][2])*(yi-parm[j][2]));
00115                 if(d >= apertures[i]+apertures[j]) {
00116                     aa[j][i] = 0.0;
00117                     aa[i][j] = aa[j][i];
00118                 } else {
00119                     a2i = apertures[i]*apertures[i];
00120                     a2j = apertures[j]*apertures[j];
00121                     di = 0.5*(d + (a2i - a2j)/d);
00122                     dj = 0.5*(d - (a2i - a2j)/d);
00123                     argi = di/apertures[i];
00124                     argj = dj/apertures[j];
00125                     if (di < 0.0) {
00126                         argi = 1.0;
00127                         argj = 0.0;
00128                     }
00129                     if (dj < 0.0) {
00130                         argi = 0.0;
00131                         argj = 1.0;
00132                     }
00133                     aa[j][i] = cn[i]*cn[j]*
00134                         (a2i*(acosf(argi) - argi*(sqrtf(1.0 - argi*argi))) +
00135                          a2j*(acosf(argj) - argj*(sqrtf(1.0 - argj*argj))));
00136                     aa[i][j] = aa[j][i];
00137                 }
00138             }
00139         }
00140     }
00141 
00142     /* Clear accumulators */
00143 
00144     for (i = 0; i < nbit; i++)
00145         bb[i] = 0.0;
00146 
00147     /* generate image-blend outer boundaries */
00148 
00149     xmin = 1.0e6;
00150     xmax = -1.0e6;
00151     ymin = 1.0e6;
00152     ymax = -1.0e6;
00153     for(i = 0; i < nbit; i++) {
00154         xi = parm[i][1];
00155         yi = parm[i][2];
00156         xmin = MIN(xmin, xi - parrad[i]);
00157         xmax = MAX(xmax, xi + parrad[i]);
00158         ymin = MIN(ymin, yi - parrad[i]);
00159         ymax = MAX(ymax, yi + parrad[i]);
00160     }
00161     ix1 = MAX(0,(int)xmin-1);
00162     ix2 = MIN(nx-1,(int)xmax);
00163     iy1 = MAX(0,(int)ymin-1);
00164     iy2 = MIN(ny-1,(int)ymax);
00165 
00166 
00167     /* Now go through pixel region */
00168 
00169     for(ii = iy1; ii <= iy2; ii++) {
00170         kk = ii*nx;
00171         for(i = ix1; i <= ix2; i++) {
00172             if (mflag[kk+i] == MF_CLEANPIX || mflag[kk+i] == MF_OBJPIX) {
00173                 t = map[kk+i];   
00174                 for(j = 0; j < nbit; j++) {
00175                     xj = i - parm[j][1] + 1.0;
00176                     yj = ii - parm[j][2] + 1.0;
00177                     bb[j] += fraction(xj,yj,apertures[j])*t;
00178                 }
00179             } else {
00180                 for (j = 0; j < nbit; j++) {
00181                     xj = i - parm[j][1] + 1.0;
00182                     yj = ii - parm[j][2] + 1.0;
00183                     tj = fraction(xj,yj,apertures[j]);
00184                     aa[j][j] -= tj*tj*cn[j]*cn[j];
00185                     for (k = j + 1; k < nbit; k++) {
00186                         xk = i - parm[k][1] + 1.0;
00187                         yk = ii - parm[k][2] + 1.0;
00188                         tk = fraction(xk,yk,apertures[k]);
00189                         aa[j][k] -= tk*tj*cn[j]*cn[k];
00190                         aa[k][j] = aa[j][k];
00191                     }
00192                 }
00193             } 
00194         }
00195     }
00196 
00197     /* Supress matrix instabilities */
00198 
00199     if (nbit > 1) {
00200         for (k = 0; k < nbit; k++) {
00201             for (j = k + 1; j < nbit; j++) {
00202                 aa[j][k] = 0.99*MIN(aa[k][k],aa[j][j]);
00203                 aa[k][j] = aa[j][k];
00204             }
00205         }
00206     }
00207 
00208     /* Trivial solution for single object */
00209 
00210     if (nbit == 1) {
00211         cflux[0] = bb[0];
00212         return;
00213     }
00214 
00215     /* solve for profile intensities */
00216 
00217     dchole(aa,bb,nbit);
00218     for(i = 0; i < nbit; i++) 
00219         cflux[i] = cn[i]*bb[i];
00220 }
00221 
00222 
00223 /* returns fraction of pixel bounded by 0 -  r_out
00224  * x,y coordinates relative to centre
00225  * Uses linear approximation ok if pixel located >>1 away from centre */
00226 
00227 static float fraction (float x, float y, float r_out) {
00228     float r,t,x_a,x_b,frac,tanao2,cosa,tanp2a,sqrt2o2;
00229 
00230     r = sqrtf(x*x + y*y);
00231     sqrt2o2 = 0.5*M_SQRT2;
00232 
00233     /* is it worth bothering? */
00234 
00235     if(r > r_out+sqrt2o2) 
00236         return(0.0);
00237 
00238     /* is it trivially all in? */
00239 
00240     if(r < r_out-sqrt2o2) 
00241         return(1.0);
00242 
00243     /* bugger - have to do some work then ... ok first ...
00244      * use 8-fold symmetry to convert to 0-45 degree range */
00245 
00246     x = fabsf(x);
00247     y = fabsf(y);
00248     if(y > x) {
00249         t = x;
00250         x = y;
00251         y = t;
00252     }
00253 
00254     /* If the angles are too close to cardinal points, then fudge something */
00255 
00256     if (x > 0.0 && y > 0.0) {
00257         tanao2 = 0.5*y/x;
00258         tanp2a = x/y;
00259         cosa = x/sqrt(x*x + y*y);
00260     } else {
00261         tanao2 = 0.00005;
00262         tanp2a = 10000.0;
00263         cosa = 1.0;
00264     }
00265 
00266     /* only outer radius - compute linear intersections top and bot of pixel */
00267 
00268     x_a = x - tanao2 + (r_out - r)/cosa;
00269     if(x_a < x+0.5) {
00270 
00271         /* intersects */
00272 
00273         x_b = x + tanao2 + (r_out - r)/cosa;
00274 
00275         /* three cases to consider */
00276 
00277         if(x_a < x-0.5)
00278             frac = 0.5*MAX(0.0,x_b-(x-0.5))*MAX(0.0,x_b-(x-0.5))*tanp2a;
00279         else {
00280             if(x_b > x+0.5)
00281                 frac = 1.0 - 0.5*(x+0.5-x_a)*(x+0.5-x_a)*tanp2a;
00282             else
00283                 frac = 0.5-(x-x_a)+0.5*(x_b-x_a);
00284         }
00285     } else  /* missed entirely */
00286         frac = 1.0;
00287 
00288     return(frac);
00289 }
00290 
00291 /* CHOLEsky decomposition of +ve definite symmetric matrix to solve Ax = b */
00292 
00293 static void dchole (double a[IMNUM+1][IMNUM+1], double b[IMNUM+1], int n) {
00294   double sum, l[IMNUM+1][IMNUM+1], y[IMNUM+1];
00295   double aveigv, offset;
00296   int i, j, k;
00297 
00298 restart:
00299     l[0][0] = sqrt(a[0][0]);
00300 
00301     for(k = 1; k < n; k++) {
00302         for(j = 0; j <= k-1; j++) {
00303             sum = a[j][k];
00304             if(j != 0) 
00305                 for(i = 0; i <= j-1; i++) 
00306                     sum -= l[i][k]*l[i][j];
00307             l[j][k] = sum/l[j][j];
00308         }
00309         sum = a[k][k];
00310         for(i = 0; i <= k-1; i++) 
00311             sum -= l[i][k]*l[i][k];
00312         if(sum <= 0.0) {
00313 /*          fprintf(stderr, "dchole: warning: matrix ill-conditioned\n"); */
00314             aveigv = a[0][0];
00315             for(i = 1; i < n; i++) 
00316                 aveigv += a[i][i];
00317             /* max eigenvalue < trace */
00318             offset = 0.1*aveigv/((double) n);
00319             for(i = 0; i < n; i++) 
00320                 a[i][i] += offset;
00321 /*          fprintf(stderr, "dchole: Offset added to diagonal = %f\n", offset); */
00322             goto restart;
00323         }
00324         l[k][k] = sqrt(sum);
00325     }
00326 
00327     /* solve Ly = b */
00328 
00329     y[0] = b[0]/l[0][0];
00330     for(i = 1; i < n; i++) {
00331         sum = b[i];
00332         for(k = 0; k <= i-1; k++) 
00333             sum -= l[k][i]*y[k];
00334         y[i] = sum/l[i][i];
00335     }
00336 
00337     /* solve L(T)x = y */
00338 
00339     b[n-1] = y[n-1]/l[n-1][n-1];
00340     for(i = n-2; i >= 0; i--) {
00341         sum = y[i];
00342         for(k = i+1; k < n; k++) 
00343             sum -= l[i][k]*b[k];
00344         b[i] = sum/l[i][i];
00345     }
00346 }
00347 
00348 /*
00349 
00350 $Log: imcore_radii.c,v $
00351 Revision 1.1  2005/09/13 13:25:30  jim
00352 Initial entry after modifications to make cpl compliant
00353 
00354 
00355 */

Generated on Sat Apr 6 04:03:07 2013 for VIRCAM Pipeline by  doxygen 1.5.1