Changeset 6399

Timestamp:

Feb 9, 2006, 9:40:40 AM (20 years ago)

Author:

eugene

Message:

added background median map

Location:

trunk/psphot/src

Files:

• Makefile.am (modified) (1 diff)
• psphot.h (modified) (1 diff)
• psphotCleanup.c (added)
• psphotImageMedian.c (modified) (2 diffs)
• psphotReadout.c (modified) (1 diff)

trunk/psphot/src/Makefile.am

@@ -39 +39 @@
         psphotMagnitudes.c      \
         psphotRadiusChecks.c    \
+        psphotImageMedian.c     \
         psLine.c                \
         psSparse.c              \

trunk/psphot/src/psphot.h

@@ -39 +39 @@
 pmPSF          *psphotChoosePSF (psMetadata *config, psArray *sources, psStats *sky);
 void            psphotOutput (pmReadout *readout, psMetadata *arguments);
+psImage        *psphotImageMedian (pmReadout *readout, psMetadata *config, psStats *sky);
 
 // optional object analysis steps

trunk/psphot/src/psphotImageMedian.c

@@ -1 +1 @@
 # include "psphot.h"
+double psImageClippedStats (psImage *image, psImage *mask, psU8 maskValue, double fmin, double fmax);
+void fsort (float *value, int N);
+
+// random number seed to select a fraction of the image pixels
+static psRandom *rnd;
+
+// use no more than MAX_SAMPLE_PIXELS pixels for each median box
+static int MAX_SAMPLE_PIXELS;
 
 // generate the median in NxN boxes, clipping heavily
-// build a 2D clipped spline to the data (toss out the outliers?
-
+// linear interpolation to generate full-scale model
 psImage *psphotImageMedian (pmReadout *readout, psMetadata *config, psStats *sky) 
 { 
-
+    bool status;
     psRegion region;
+
+    rnd = psRandomAlloc (PS_RANDOM_TAUS, 0);
+    MAX_SAMPLE_PIXELS = psMetadataLookupF32 (&status, config, "IMSTATS_NPIX");
+    if (!status) MAX_SAMPLE_PIXELS = 1000;
 
     // fit background to subset of image points within 3 sigma of sky
     psTimerStart ("psphot");
-    psRandom *rnd     = psRandomAlloc (PS_RANDOM_TAUS, 0);
 
     psImage *image = readout->image;
     psImage *mask  = readout->mask;
 
-    bool      status  = false;
-    int       Nx = image->numCols + image->col0;
-    int       Ny = image->numRows + image->row0;
-    int       xbin = psMetadataLookupS32 (&status, config, "BACKGROUND_XBIN");
-    if (!status) {xbin = 64;}
-    int       ybin = psMetadataLookupS32 (&status, config, "BACKGROUND_YBIN");
-    if (!status) {ybin = 64;}
-
-    // we slide the box in half-steps from -0.5 off the left edge to +0.5 off the right edge 
-    int nx = 2 + 2 * Nx / xbin;
-    int ny = 2 + 2 * Ny / ybin;
-
-    psStats *stats = psStatsAlloc (PS_STATS_CLIPPED_MEAN, PS_STATS_CLIPPED_STDEV);
-    psImage *medimage = psImageAlloc (nx, ny, PS_TYPE_F32);
-
-    psVector *xval = psVectorAlloc (nx, PS_TYPE_F32);
-    psVector *yval = psVectorAlloc (ny, PS_TYPE_F32);
+    // dimensions of input & output image
+    int Nx = image->numCols;
+    int Ny = image->numRows;
+
+    // scaling factor
+    int DX = psMetadataLookupS32 (&status, config, "BACKGROUND_XBIN");
+    if (!status) {DX = 64;}
+    int DY = psMetadataLookupS32 (&status, config, "BACKGROUND_YBIN");
+    if (!status) {DY = 64;}
+
+    // overhang : we will balance this evenly
+    int dx = (Nx % DX) / 2;
+    int dy = (Ny % DY) / 2;
+
+    // dimensions of binned image
+    int nx = (Nx % DX) ? (int)(Nx / DX) + 1 : Nx / DX;
+    int ny = (Ny % DY) ? (int)(Ny / DY) + 1 : Ny / DY;
+
+    // psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+    psImage *model = psImageAlloc (nx, ny, PS_TYPE_F32);
 
     // measure clipped median for subimages 
-    for (int iy = -1; iy < ny - 1; iy++) {
-        for (int ix = -1; ix < nx - 1; ix++) {
-            int sx = 0.5*ix*xbin;
-            int sy = 0.5*iy*ybin;
-            region = psRegionSet (sx, sx + xbin, sy, sy + ybin);
-            region = psRegionForImage (readout->image, region);
-            psImage *subset = psImageSubset (readout->image, region);
-            stats = psImageStats (stats, subset);
-            
-            // XXX filter these on the clippedStdev?
-            medimage->data.F32[iy][ix] = stats->clippedMean;
-
-        }
-    }
-    for (int ix = -1; ix < nx - 1; ix++) {
-        xval->data.F32[ix] = 0.5*ix*xbin + 0.5*xbin;
-    }
-    for (int iy = -1; iy < ny - 1; iy++) {
-        yval->data.F32[iy] = 0.5*iy*ybin + 0.5*ybin;
-    }
-
-    // x-dir / y-dir
-    psImage *devimage = psImageFitSpline (medimage, true);
-    psImage *background = psImageEvalSpline (medimage, devimage, xval, yval, Nx, Ny);
-
-    // this is a very inefficient way to evaluate the function..
-    // lame, temporary step to zero masked pixels
+    for (int iy = 0; iy < ny; iy++) {
+        for (int ix = 0; ix < nx; ix++) {
+            // sx, sy are in parent coords
+            int sx = ix*DX - dx + image->col0;
+            int sy = iy*DY - dy + image->row0;
+            region = psRegionSet (sx, sx + 2*DX, sy, sy + 2*DY);
+            region = psRegionForImage (image, region);
+            psImage *subset  = psImageSubset (image, region);
+            psImage *submask = psImageSubset (mask, region);
+
+            model->data.F32[iy][ix] = psImageClippedStats (subset, submask, 0xff, 0.25, 0.50);
+
+            psFree (subset);
+            psFree (submask);
+
+            // XXX psImageStats is still very inefficient and poorly coded...
+            // stats = psImageStats (stats, subset, maskset, 0xff);
+            // model->data.F32[iy][ix] = stats->clippedMean;
+        }
+    }
+    psLogMsg ("psphot", 3, "build median image: %f sec\n", psTimerMark ("psphot"));
+
+    // linear interpolation to full-scale
+    
+    psImage *background = psImageAlloc (Nx, Ny, PS_TYPE_F32);
+
+    // XXX this code skips the initial pixels
+    for (int Iy = 0; Iy < ny-1; Iy ++) {
+        for (int Ix = 0; Ix < nx-1; Ix ++) {
+
+            float V00 = model->data.F32[Iy+0][Ix+0];
+            float V01 = model->data.F32[Iy+0][Ix+1];
+            float V10 = model->data.F32[Iy+1][Ix+0];
+            float V11 = model->data.F32[Iy+1][Ix+1];
+
+            // a single binned pixel quad
+            // (Xs,Ys) : (Xe,Ye) : binned pixel centers in unbinned coords
+            // corresponding to (Ix,Iy), (Ix+1,Iy+1)
+            int Xs = (Ix + 1)*DX - dx;
+            int Ys = (Iy + 1)*DY - dy;
+            int Xe = Xs + DX;
+            int Ye = Ys + DY;
+
+            for (int iy = Ys; (iy < Ye) && (iy < Ny); iy++) {
+                float Vxs = (V10 - V00)*(iy - Ys) / DY + V00;
+                float Vxe = (V11 - V01)*(iy - Ys) / DY + V01;
+                float dV = (Vxe - Vxs) / DX;
+                float V  = Vxs;
+                for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) {
+                    background->data.F32[iy][ix] = V;
+                    V += dV;
+                }
+            }
+        }
+    }
+
+    // side pixels
+    int Xs = DX - dx;
+    int Xe = nx*DX - dx;
+    for (int Iy = 0; Iy < ny - 1; Iy++) {
+
+        int Ys = (Iy + 1)*DY - dy;
+        int Ye = Ys + DY;
+
+        // leading edge
+        float V0 = model->data.F32[Iy+0][0];
+        float V1 = model->data.F32[Iy+1][0];
+        float dV = (V1 - V0) / DY;
+        float V = V0;
+        for (int iy = Ys; (iy < Ye) && (iy < Ny); iy++) {
+            for (int ix = 0; ix < Xs; ix++) {
+                background->data.F32[iy][ix] = V;
+            }
+            V += dV;
+        }
+
+        // trailing edge
+        V0 = model->data.F32[Iy+0][nx-1];
+        V1 = model->data.F32[Iy+1][nx-1];
+        dV = (V1 - V0) / DY;
+        V = V0;
+        for (int iy = Ys; (iy < Ye) && (iy < Ny); iy++) {
+            for (int ix = Xe; ix < Nx; ix++) {
+                background->data.F32[iy][ix] = V;
+            }
+            V += dV;
+        }
+    }
+
+    // top and bottom pixels
+    int Ys = DY - dy;
+    int Ye = ny*DY - dy;
+    for (int Ix = 0; Ix < nx - 1; Ix++) {
+
+        int Xs = (Ix + 1)*DX - dx;
+        int Xe = Xs + DX;
+
+        // top edge
+        float V0 = model->data.F32[0][Ix+0];
+        float V1 = model->data.F32[0][Ix+1];
+        float dV = (V1 - V0) / DX;
+        for (int iy = 0; iy < Ys; iy++) {
+            float V = V0;
+            for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) {
+                background->data.F32[iy][ix] = V;
+                V += dV;
+            }
+        }
+
+        // bottom edge
+        V0 = model->data.F32[ny-1][Ix+0];
+        V1 = model->data.F32[ny-1][Ix+1];
+        dV = (V1 - V0) / DX;
+        for (int iy = Ye; iy < Ny; iy++) {
+            float V = V0;
+            for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) {
+                background->data.F32[iy][ix] = V;
+                V += dV;
+            }
+        }
+    }
+
+    // the four corners
+    {
+        float V;
+        // 0,0
+        V = model->data.F32[0][0];
+        for (int iy = 0; iy < DY - dy; iy++) {
+            for (int ix = 0; ix < DX - dx; ix++) {
+                background->data.F32[iy][ix] = V;
+            }
+        }
+        // Nx,0
+        V = model->data.F32[0][nx-1];
+        for (int iy = 0; iy < DY - dy; iy++) {
+            for (int ix = nx*DX - dx; ix < Nx; ix++) {
+                background->data.F32[iy][ix] = V;
+            }
+        }
+        // 0,Ny
+        V = model->data.F32[ny-1][0];
+        for (int iy = ny*DY - dy; iy < Ny; iy++) {
+            for (int ix = 0; ix < DX - dx; ix++) {
+                background->data.F32[iy][ix] = V;
+            }
+        }
+        // Nx,Ny
+        V = model->data.F32[nx-1][ny-1];
+        for (int iy = ny*DY - dy; iy < Ny; iy++) {
+            for (int ix = nx*DX - dx; ix < Nx; ix++) {
+                background->data.F32[iy][ix] = V;
+            }
+        }
+    }   
+    psLogMsg ("psphot", 3, "build resampled image: %f sec\n", psTimerMark ("psphot"));
+
+    // subtract the background model
     for (int j = 0; j < image->numRows; j++) {
         for (int i = 0; i < image->numCols; i++) {
@@ -69 +212 @@
         }
     }
-    psFree (medimage);
-    psFree (devimage);
-    psFree (rnd);  // XXX should this be made available at a higher level?
-
-    psLogMsg ("psphot", 3, "fit background model: %f sec\n", psTimerMark ("psphot"));
+    psLogMsg ("psphot", 3, "subtracted background model: %f sec\n", psTimerMark ("psphot"));
+    psFree (rnd);
     return (background);
 }
+
+double psImageClippedStats (psImage *image, psImage *mask, psU8 maskValue, double fmin, double fmax) {
+
+    double value;
+    int nx = image->numCols;
+    int ny = image->numRows;
+    
+    if (nx*ny <= 0) return 0.0;
+
+    int Nsubset = PS_MIN (MAX_SAMPLE_PIXELS, nx*ny);
+    int Npixels = nx*ny;
+
+    psVector *values = psVectorAlloc (Nsubset, PS_TYPE_F32);
+
+    float min = values->data.F32[0];
+    float max = values->data.F32[0];
+
+    int n = 0;
+    for (int i = 0; i < Nsubset; i++) {
+        double frnd = psRandomUniform (rnd);
+        int pixel = Npixels * frnd;
+        int ix = pixel % nx;
+        int iy = pixel / nx;
+
+        if (mask->data.U8[iy][ix] & maskValue) continue;
+
+        value = image->data.F32[iy][ix];
+        min = PS_MIN (value, min);
+        max = PS_MIN (value, max);
+        values->data.F32[n] = value;
+        n++;
+    }
+    values->n = n;
+
+    int imin = fmin * n;
+    int imax = fmax * n;
+    int npts = imax - imin + 1;
+
+    fsort (values->data.F32, n);
+
+    value = 0;
+    for (int i = imin; (i <= imax) && (i < n); i++) {
+        value += values->data.F32[i];
+    }
+    value = value / npts;
+    
+    psFree (values);
+    return value;
+}
+
+void fsort (float *value, int N) {
+
+    int l,j,ir,i;
+    float temp;
+  
+    if (N < 2) return;
+    l = N >> 1;
+    ir = N - 1;
+    for (;;) {
+        if (l > 0) {
+            temp = value[--l];
+        }
+        else {
+            temp = value[ir];
+            value[ir] = value[0];
+            if (--ir == 0) {
+                value[0] = temp;
+                return;
+            }
+        }
+        i = l;
+        j = (l << 1) + 1;
+        while (j <= ir) {
+            if (j < ir && value[j] < value[j+1]) ++j;
+            if (temp < value[j]) {
+                value[i]=value[j];
+                j += (i=j) + 1;
+            }
+            else j = ir + 1;
+        }
+        value[i] = temp;
+    }
+}
+

trunk/psphot/src/psphotReadout.c

@@ -11 +11 @@
     psPolynomial2D *model = NULL;
 
-    psphotSaveImage (NULL, readout->mask, "mask.fits");
-    psphotSaveImage (NULL, readout->image, "image.fits");
-    psphotSaveImage (NULL, readout->weight, "weight.fits");
+    psphotImageMedian (readout, config, sky);
 
     // measure image stats for initial guess