moments.c

00001 /*
00002 
00003 $Id: moments.c,v 1.2 2006/08/21 09:07:29 jim Exp $
00004 
00005 */
00006 
00007 #include <stdio.h>
00008 #include <stdlib.h>
00009 
00010 #include "imcore.h"
00011 #include "util.h"
00012 
00013 extern void moments(ap_t *ap, float results[]) {
00014     int i,np;
00015     float x,y,xoff,yoff,xsum,ysum,xsumsq,ysumsq,tsum,xysum,t,tmax,twelfth;
00016     float xbar,ybar,sxx,syy,sxy,xintmin,w,wsum,xsum_w,ysum_w;
00017     plstruct *plarray;
00018 
00019     /* Initialise a few things */
00020 
00021     xintmin = ap->xintmin;
00022     plarray = ap->plarray;
00023     np = ap->npl_pix;
00024     xoff = (float)plarray[0].x;
00025     yoff = (float)plarray[0].y;
00026     xsum = 0.0;
00027     ysum = 0.0;
00028     xsum_w = 0.0;
00029     ysum_w = 0.0;
00030     wsum = 0.0;
00031     xsumsq = 0.0;
00032     ysumsq = 0.0;
00033     tsum = 0.0;
00034     xysum = 0.0;
00035     tmax = plarray[0].z;
00036     twelfth = 1.0/12.0;
00037 
00038     /* Do a moments analysis on an object */
00039 
00040     for (i = 0; i < np; i++) {
00041         x = (float)plarray[i].x - xoff;
00042         y = (float)plarray[i].y - yoff;
00043         t = plarray[i].z;
00044         w = plarray[i].zsm; 
00045         if (t < 0.0)
00046             continue;
00047         xsum += t*x;
00048         ysum += t*y;
00049         tsum += t;
00050         xsum_w += w*t*x;
00051         ysum_w += w*t*y;
00052         wsum += w*t;
00053         tmax = MAX(tmax,plarray[i].z);
00054         xsumsq += (x*x + twelfth)*t;
00055         ysumsq += (y*y + twelfth)*t;
00056         xysum += x*y*t;
00057     }
00058 
00059     /* Check that the total intensity is enough and if it is, then do
00060        the final results */
00061 
00062     if (tsum >= xintmin) {
00063         xbar = xsum/tsum;
00064         ybar = ysum/tsum;
00065         sxx = MAX(0.0,(xsumsq/tsum-xbar*xbar));
00066         syy = MAX(0.0,(ysumsq/tsum-ybar*ybar));
00067         sxy = xysum/tsum - xbar*ybar;
00068         xbar = xsum_w/wsum;
00069         ybar = ysum_w/wsum;
00070         xbar += xoff;
00071         ybar += yoff;
00072         xbar = MAX(1.0,MIN(xbar,ap->lsiz));
00073         ybar = MAX(1.0,MIN(ybar,ap->csiz));
00074         results[0] = 1.0;
00075         results[1] = xbar;
00076         results[2] = ybar;
00077         results[3] = tsum;
00078         results[4] = sxx;
00079         results[5] = sxy;
00080         results[6] = syy;
00081         results[7] = tmax;
00082     } else {
00083         results[0] = -1.0;
00084     }
00085 }
00086 
00087 /*
00088 
00089 $Log: moments.c,v $
00090 Revision 1.2  2006/08/21 09:07:29  jim
00091 Modified centring algorithm
00092 
00093 Revision 1.1  2005/09/13 13:25:30  jim
00094 Initial entry after modifications to make cpl compliant
00095 
00096 
00097 */

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