Changeset 7385 for trunk/psphot/src/psphotImageMedian.c

Timestamp:

Jun 6, 2006, 5:37:11 PM (20 years ago)

Author:

Paul Price

Message:

Moving pslib additional functions from psModules into psLib proper; some API changes fixed.

File:

• trunk/psphot/src/psphotImageMedian.c (modified) (5 diffs)

trunk/psphot/src/psphotImageMedian.c

@@ -3 +3 @@
 // generate the median in NxN boxes, clipping heavily
 // linear interpolation to generate full-scale model
-bool psphotImageMedian (pmConfig *config, pmFPAview *view) 
-{ 
+bool psphotImageMedian (pmConfig *config, pmFPAview *view)
+{
     bool status;
     psRegion region;
@@ -15 +15 @@
 
     MAX_SAMPLE_PIXELS = psMetadataLookupF32 (&status, recipe, "IMSTATS_NPIX");
-    if (!status) MAX_SAMPLE_PIXELS = 1000;
+    if (!status) {
+        MAX_SAMPLE_PIXELS = 1000;
+    }
     unsigned long seed = psMetadataLookupS32 (&status, recipe, "IMSTATS_SEED");
-    if (!status) seed = 0;
-    psImageClippedStatsInit(MAX_SAMPLE_PIXELS, seed);
+    if (!status) {
+        seed = 0;
+    }
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, seed);
 
     // subtract this amount extra from the sky
     float SKY_BIAS = psMetadataLookupF32 (&status, recipe, "SKY_BIAS");
-    if (!status) SKY_BIAS = 0;
+    if (!status) {
+        SKY_BIAS = 0;
+    }
 
     // find the currently selected readout
@@ -49 +55 @@
     pmReadout *model = pmFPAfileThisReadout (config->files, view, "PSPHOT.BACKMDL");
     if (model == NULL) {
-        // select model pixels, from output background model file, or create
-        model = pmFPAfileCreateInternal (config->files, "PSPHOT.BACKMDL", nx, ny, PS_TYPE_F32);
+        // select model pixels, from output background model file, or create
+        model = pmFPAfileCreateInternal (config->files, "PSPHOT.BACKMDL", nx, ny, PS_TYPE_F32);
     } else {
-        // replace the supplied image data with an image of the desired size
-        psImageRecycle (model->image, nx, ny, PS_TYPE_F32);
+        // replace the supplied image data with an image of the desired size
+        psImageRecycle (model->image, nx, ny, PS_TYPE_F32);
     }
     psF32 **modelData = model->image->data.F32;
 
-    // measure clipped median for subimages 
+    // measure clipped median for subimages
     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);
+        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);
 
-            // XXX the value of the upper and lower cuts probably should be studied...
-            psStats *stats = psImageClippedStats (subset, submask, 0xff, 0.25, 0.75);
-            modelData[iy][ix] = stats->robustMedian + SKY_BIAS;
+            // XXX the value of the upper and lower cuts probably should be studied...
+            psStats *stats = psImageBackground(subset, submask, 0xff, 0.25, 0.75, MAX_SAMPLE_PIXELS, rng);
+            modelData[iy][ix] = stats->robustMedian + SKY_BIAS;
 
-            psFree (stats);
-            psFree (subset);
-            psFree (submask);
+            psFree (stats);
+            psFree (subset);
+            psFree (submask);
 
-            // XXX psImageStats is still very inefficient and poorly coded...
-            // psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
-            // stats = psImageStats (stats, subset, maskset, 0xff);
-            // backMdl->data.F32[iy][ix] = stats->clippedMean;
-        }
+            // XXX psImageStats is still very inefficient and poorly coded...
+            // psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+            // stats = psImageStats (stats, subset, maskset, 0xff);
+            // backMdl->data.F32[iy][ix] = stats->clippedMean;
+        }
     }
     psLogMsg ("psphot", 3, "build median image: %f sec\n", psTimerMark ("psphot"));
@@ -98 +104 @@
     pmReadout *background = pmFPAfileThisReadout (config->files, view, "PSPHOT.BACKGND");
     if (background == NULL) {
-        background = pmFPAfileCreateInternal (config->files, "PSPHOT.BACKGND", Nx, Ny, PS_TYPE_F32);
-    } 
+        background = pmFPAfileCreateInternal (config->files, "PSPHOT.BACKGND", Nx, Ny, PS_TYPE_F32);
+    }
     psF32 **backData = background->image->data.F32;
 
@@ -115 +121 @@
     // subtract the background model (save in backSub, if requested)
     for (int j = 0; j < image->numRows; j++) {
-        for (int i = 0; i < image->numCols; i++) {
-            if (!mask->data.U8[j][i]) {
-                image->data.F32[j][i] -= backData[j][i];
-                if (backSub) {
-                    backSub->image->data.F32[j][i] = image->data.F32[j][i];
-                }
-            }
-        }
+        for (int i = 0; i < image->numCols; i++) {
+            if (!mask->data.U8[j][i]) {
+                image->data.F32[j][i] -= backData[j][i];
+                if (backSub) {
+                    backSub->image->data.F32[j][i] = image->data.F32[j][i];
+                }
+            }
+        }
     }
 
     psLogMsg ("psphot", 3, "subtracted background model: %f sec\n", psTimerMark ("psphot"));
-    psImageClippedStatsCleanup();
+    psFree(rng);
 
     // the pmReadout selected in this function are all view on entries in config->files