imcore_phopt.c

00001 /*
00002 
00003 $Id: imcore_phopt.c,v 1.1 2005/09/13 13:25:29 jim Exp $
00004 
00005 */
00006 
00007 #include <stdio.h>
00008 #include <stdlib.h>
00009 #include <math.h>
00010 
00011 #include <cpl.h>
00012 
00013 #include "floatmath.h"
00014 #include "util.h"
00015 #include "imcore.h"
00016 #include "ap.h"
00017 
00018 /* Function Prototypes */
00019 
00020 static void dchole (double a[IMNUM+1][IMNUM+1], double b[IMNUM+1], int n);
00021 static float fraction (float x, float y, float r_out);
00022 
00023 /* does multiple profile fitting to determine intensities */
00024 
00025 extern void phopt(ap_t *ap, float parm[IMNUM][NPAR], int nbit, int naper, 
00026                   float apertures[], float cflux[], float badpix[], 
00027                   int nrcore) {
00028     double aa[IMNUM+1][IMNUM+1],bb[IMNUM+1];
00029     float d,arg,*map,rcirc;
00030     float cn,parrad,xmin,xmax,xi,yi,ymin,ymax;
00031     float t,xj,yj,cnsq,tj,xk,yk,tk;
00032     int i,ii,j,kk,ix1,ix2,iy1,iy2,nx,ny,k,iaper;
00033     unsigned char *mflag;
00034 
00035     /* Set up some local variables */
00036 
00037     map = ap->indata;
00038     mflag = ap->mflag;
00039     nx = ap->lsiz;
00040     ny = ap->csiz;   
00041 
00042     /* Loop for each of the apertures */
00043 
00044     for (iaper = 0; iaper < naper; iaper++) {
00045         rcirc = apertures[iaper];
00046         parrad = rcirc + 0.5;
00047         cn = 1.0/(M_PI*rcirc*rcirc);       /* profile normalising constant */
00048         cnsq = cn*cn;
00049     
00050         /* set up covariance matrix - analytic special case for cores */
00051 
00052         for(i = 0; i < nbit; i++) {
00053             aa[i][i] = cn;                 /* overlaps totally area=pi*r**2 */
00054             if(nbit > 1) {
00055                 xi = parm[i][1];
00056                 yi = parm[i][2];
00057                 for(j = i+1; j < nbit; j++) {
00058                     d = sqrtf((xi-parm[j][1])*(xi-parm[j][1]) 
00059                         + (yi-parm[j][2])*(yi-parm[j][2]));
00060                     if(d >= 2.0*rcirc) {
00061                         aa[j][i] = 0.0;
00062                         aa[i][j] = aa[j][i];
00063                     } else {
00064                         arg = d/(2.0*rcirc);
00065                         aa[j][i] = cnsq*2.0*rcirc*rcirc*
00066                             (acosf(arg)-arg*(sqrtf(1.0-arg*arg)));
00067                         aa[i][j] = aa[j][i];
00068                     }
00069                 }
00070             }
00071         }
00072 
00073         /* clear accumulators */
00074 
00075         for(i = 0; i < nbit; i++) 
00076             bb[i] = 0.0;
00077 
00078         /* generate image-blend outer boundaries */
00079 
00080         xmin = 1.0e6;
00081         xmax = -1.0e6;
00082         ymin = 1.0e6;
00083         ymax = -1.0e6;
00084         for(i = 0; i < nbit; i++) {
00085             xi = parm[i][1];
00086             yi = parm[i][2];
00087             xmin = MIN(xmin, xi);
00088             xmax = MAX(xmax, xi);
00089             ymin = MIN(ymin, yi);
00090             ymax = MAX(ymax, yi);
00091         }
00092         ix1 = MAX(0,(int)(xmin-parrad)-1);
00093         ix2 = MIN(nx-1,(int)(xmax+parrad));
00094         iy1 = MAX(0,(int)(ymin-parrad)-1);
00095         iy2 = MIN(ny-1,(int)(ymax+parrad));
00096 
00097         /* now go through pixel region */
00098 
00099         for(ii = iy1; ii <= iy2; ii++) {
00100             kk = ii*nx;
00101             for(i = ix1; i <= ix2; i++) {
00102                 if (mflag[kk+i] == MF_CLEANPIX || mflag[kk+i] == MF_OBJPIX) {
00103                     t = map[kk+i];   
00104                     for(j = 0; j < nbit; j++) {
00105                         xj = i - parm[j][1] + 1.0;
00106                         yj = ii - parm[j][2] + 1.0;
00107                         bb[j] += fraction(xj,yj,rcirc)*t;
00108                     }
00109                 } else {
00110                     for (j = 0; j < nbit; j++) {
00111                         xj = i - parm[j][1] + 1.0;
00112                         yj = ii - parm[j][2] + 1.0;
00113                         tj = fraction(xj,yj,rcirc);
00114                         aa[j][j] -= tj*tj*cnsq;
00115                         for (k = j + 1; k < nbit; k++) {
00116                             xk = i - parm[k][1] + 1.0;
00117                             yk = ii - parm[k][2] + 1.0;
00118                             tk = fraction(xk,yk,rcirc);
00119                             aa[j][k] -= tk*tj*cnsq;
00120                             aa[k][j] = aa[j][k];
00121                         }
00122                         if (iaper == nrcore)
00123                             badpix[j] += tj;
00124                     }
00125                 } 
00126             }
00127         }
00128 
00129         /* Trivial solution for single object */
00130 
00131         if (nbit == 1) {
00132             cflux[iaper] = bb[0];
00133 
00134         /* solve for profile intensities */
00135 
00136         } else {
00137             for (i = 0; i < nbit; i++) 
00138                 aa[i][i] = MAX(aa[i][i],cnsq);
00139             dchole(aa,bb,nbit);
00140             for(i = 0; i < nbit; i++) 
00141                 cflux[i*naper+iaper] = cn*bb[i];
00142         }
00143     }
00144 }
00145 
00146 /* CHOLEsky decomposition of +ve definite symmetric matrix to solve Ax = b */
00147 
00148 static void dchole (double a[IMNUM+1][IMNUM+1], double b[IMNUM+1], int n) {
00149   double sum, l[IMNUM+1][IMNUM+1], y[IMNUM+1];
00150   double aveigv, offset;
00151   int i, j, k;
00152 
00153 restart:
00154     l[0][0] = sqrt(a[0][0]);
00155 
00156     for(k = 1; k < n; k++) {
00157         for(j = 0; j <= k-1; j++) {
00158             sum = a[j][k];
00159             if(j != 0) 
00160                 for(i = 0; i <= j-1; i++) 
00161                     sum -= l[i][k]*l[i][j];
00162             l[j][k] = sum/l[j][j];
00163         }
00164         sum = a[k][k];
00165         for(i = 0; i <= k-1; i++) 
00166             sum -= l[i][k]*l[i][k];
00167         if(sum <= 0.0) {
00168 /*          fprintf(stderr, "dchole: warning: matrix ill-conditioned\n"); */
00169             aveigv = a[0][0];
00170             for(i = 1; i < n; i++) 
00171                 aveigv += a[i][i];
00172             /* max eigenvalue < trace */
00173             offset = 0.1*aveigv/((double) n);
00174             for(i = 0; i < n; i++) 
00175                 a[i][i] += offset;
00176 /*          fprintf(stderr, "dchole: Offset added to diagonal = %f\n", offset); */
00177             goto restart;
00178         }
00179         l[k][k] = sqrt(sum);
00180     }
00181 
00182     /* solve Ly = b */
00183 
00184     y[0] = b[0]/l[0][0];
00185     for(i = 1; i < n; i++) {
00186         sum = b[i];
00187         for(k = 0; k <= i-1; k++) 
00188             sum -= l[k][i]*y[k];
00189         y[i] = sum/l[i][i];
00190     }
00191 
00192     /* solve L(T)x = y */
00193 
00194     b[n-1] = y[n-1]/l[n-1][n-1];
00195     for(i = n-2; i >= 0; i--) {
00196         sum = y[i];
00197         for(k = i+1; k < n; k++) 
00198             sum -= l[i][k]*b[k];
00199         b[i] = sum/l[i][i];
00200     }
00201 }
00202 
00203 /* returns fraction of pixel bounded by 0 -  r_out
00204  * x,y coordinates relative to centre
00205  * Uses linear approximation ok if pixel located >>1 away from centre */
00206 
00207 static float fraction (float x, float y, float r_out) {
00208     float r,t,x_a,x_b,frac,tanao2,cosa,tanp2a,sqrt2o2;
00209 
00210     r = sqrtf(x*x + y*y);
00211     sqrt2o2 = 0.5*M_SQRT2;
00212 
00213     /* is it worth bothering? */
00214 
00215     if(r > r_out+sqrt2o2) 
00216         return(0.0);
00217 
00218     /* is it trivially all in? */
00219 
00220     if(r < r_out-sqrt2o2) 
00221         return(1.0);
00222 
00223     /* bugger - have to do some work then ... ok first ...
00224      * use 8-fold symmetry to convert to 0-45 degree range */
00225 
00226     x = fabsf(x);
00227     y = fabsf(y);
00228     if(y > x) {
00229         t = x;
00230         x = y;
00231         y = t;
00232     }
00233 
00234     /* If the angles are too close to cardinal points, then fudge something */
00235 
00236     if (x > 0.0 && y > 0.0) {
00237         tanao2 = 0.5*y/x;
00238         tanp2a = x/y;
00239         cosa = x/sqrt(x*x + y*y);
00240     } else {
00241         tanao2 = 0.00005;
00242         tanp2a = 10000.0;
00243         cosa = 1.0;
00244     }
00245 
00246     /* only outer radius - compute linear intersections top and bot of pixel */
00247 
00248     x_a = x - tanao2 + (r_out - r)/cosa;
00249     if(x_a < x+0.5) {
00250 
00251         /* intersects */
00252 
00253         x_b = x + tanao2 + (r_out - r)/cosa;
00254 
00255         /* three cases to consider */
00256 
00257         if(x_a < x-0.5)
00258             frac = 0.5*MAX(0.0,x_b-(x-0.5))*MAX(0.0,x_b-(x-0.5))*tanp2a;
00259         else {
00260             if(x_b > x+0.5)
00261                 frac = 1.0 - 0.5*(x+0.5-x_a)*(x+0.5-x_a)*tanp2a;
00262             else
00263                 frac = 0.5-(x-x_a)+0.5*(x_b-x_a);
00264         }
00265     } else  /* missed entirely */
00266         frac = 1.0;
00267 
00268     return(frac);
00269 }
00270 
00271 /*
00272 
00273 $Log: imcore_phopt.c,v $
00274 Revision 1.1  2005/09/13 13:25:29  jim
00275 Initial entry after modifications to make cpl compliant
00276 
00277 
00278 */

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