Changeset 4215

Timestamp:

Jun 11, 2005, 5:43:47 AM (21 years ago)

Author:

eugene

Message:

adding mask & noise images

Location:

trunk/psphot/src

Files:

• LocalSky.c (added)
• choose_psf_model.c (modified) (1 diff)
• find_peaks.c (modified) (2 diffs)
• image_stats.c (modified) (1 diff)
• psphot-utils.c (modified) (1 diff)
• psphot.c (modified) (2 diffs)
• psphot.h (modified) (1 diff)
• setup.c (modified) (1 diff)
• source_moments.c (modified) (3 diffs)

trunk/psphot/src/choose_psf_model.c

@@ -26 +26 @@
     float RADIUS   = psMetadataLookupF32 (&status, config, "FIT_RADIUS");
 
-    float XBORDER  = psMetadataLookupF32 (&status, config, "XBORDER");
-    float YBORDER  = psMetadataLookupF32 (&status, config, "YBORDER");
-    psRegion *keep = psRegionAlloc (XBORDER, -XBORDER, YBORDER, -YBORDER);
-    keep           = psRegionForImage (keep, image, keep);
-
     for (int i = 0; i < stars->n; i++) {
         psSource *source = stars->data[i];
         pmSourceSetPixelCircle (source, image, RADIUS);
-        pmSourceMaskRegion (source, keep);
-        // XXX EAM - this is silly: reconsider how we use the masks
-        // XXX - we are needlessly re-defining the pixels here
     }
 

trunk/psphot/src/find_peaks.c

@@ -1 +1 @@
 # include "psphot.h"
 
-psArray *find_peaks (psImage *image, psMetadata *config, psStats *sky) 
+psArray *find_peaks (psImageData *imdata, psMetadata *config, psStats *sky) 
 {
     bool  status = false;
@@ -16 +16 @@
     NSIGMA = psMetadataLookupF32 (&status, config, "SMOOTH_NSIGMA");
 
+    psImage *image  = imdata->image;
     psImage *smooth = psImageCopy (NULL, image, PS_TYPE_F32);
     psImageSmooth (smooth, SIGMA, NSIGMA);

trunk/psphot/src/image_stats.c

@@ -1 +1 @@
 # include "psphot.h"
 
-psStats *image_stats (psImage *image, psMetadata *config) 
+psStats *image_stats (psImageData *imdata, psMetadata *config) 
 { 
     psStats *stats = NULL;
     psStats *sky   = NULL;
+
+    psImage *image = imdata->image;
 
     // get image stats on a subset of the image (random 1e5 pts)

trunk/psphot/src/psphot-utils.c

@@ -160 +160 @@
 }
 
+// mask the area contained by the region
+void psImageMaskRegion (psImage *image, psRegion *region, int maskValue) {
+
+    for (int iy = 0; iy < image->numRows; iy++) {
+        for (int ix = 0; ix < image->numCols; ix++) {
+            if (ix + image->col0 <  region->x0) continue;
+            if (ix + image->col0 >= region->x1) continue;
+            if (iy + image->row0 <  region->y0) continue;
+            if (iy + image->row0 >= region->y1) continue;
+            image->data.U8[iy][ix] |= maskValue;
+        }
+    }
+}
+
+// mask the area not contained by the region
+void psImageKeepRegion (psImage *image, psRegion *region, int maskValue) {
+
+    for (int iy = 0; iy < image->numRows; iy++) {
+        for (int ix = 0; ix < image->numCols; ix++) {
+            if (ix + image->col0 <  region->x0) goto maskit;
+            if (ix + image->col0 >= region->x1) goto maskit;
+            if (iy + image->row0 <  region->y0) goto maskit;
+            if (iy + image->row0 >= region->y1) goto maskit;
+            continue;
+        maskit:
+            image->data.U8[iy][ix] |= maskValue;
+        }
+    }
+}
+
 // create a subset of sources for the PS_SOURCE_PSFSTAR entries
 // XXX deprecated

trunk/psphot/src/psphot.c

@@ -4 +4 @@
 int main (int argc, char **argv) {
 
-    psMetadata *config  = NULL;         // user-provided configuration information
-    psMetadata *header  = NULL;         // input image header
-    psImage    *image   = NULL;         // input image data
-    psImage    *noise   = NULL;         // input image data
-    psImage    *mask    = NULL;         // input image data
-    psArray    *sources = NULL;
-    psArray    *peaks   = NULL;         // a list of pmPeaks 
-    pmPSF      *psf     = NULL;
-    psStats    *sky     = NULL;
+    psMetadata  *config  = NULL;                // user-provided configuration information
+    psImageData *imdata  = NULL;
+    psArray     *sources = NULL;
+    psArray     *peaks   = NULL;                // a list of pmPeaks 
+    pmPSF       *psf     = NULL;
+    psStats     *sky     = NULL;
 
     psLogArguments (&argc, argv);
@@ -20 +17 @@
     // load input data (image and config)
     // create or load mask and noise images
-    setup (&image, &mask, &noise, &header, config);
+    imdata = setup (config);
 
     // measure image stats for initial guess 
-    sky = image_stats (image, config);
+    sky = image_stats (imdata, config);
 
     // find the peaks in the image
-    peaks = find_peaks (image, config, sky);
+    peaks = find_peaks (imdata, config, sky);
 
     // construct sources and measure basic stats
-    sources = source_moments (image, config, peaks);
+    sources = source_moments (imdata, config, peaks);
 
     // source analysis is done in S/N order (brightest first)

trunk/psphot/src/psphot.h

@@ -5 +5 @@
 # include <unistd.h>
 # include <stdlib.h>
+
+typedef struct {
+  psImage *image;
+  psImage *mask;
+  psImage *noise;
+  psMetadata *header;
+} psImageData;
 
 // data to test a given PSF model type

trunk/psphot/src/setup.c

@@ -1 +1 @@
 # include "psphot.h"
 
-bool setup (psImage **image, psImage **mask, psImage **noise, psMetadata **header, psMetadata **config) 
+psImageData *setup (psMetadata *config) 
 {
-    int       Nfail;
-    psFits   *file = NULL;              // FITS pointer to input image
     psRegion  region = {0,0,0,0};       // a region representing the entire array
 
     psTimerStart ("psphot");
 
-    char *input = pmConfigLookupSTR (&status, *config, *header, "INPUT");
-    file = psFitsAlloc (input);
+    char *input = psMetadataLookupSTR (&status, config, "INPUT");
+    psFits *file = psFitsAlloc (input);
 
-    *header = psMetadataAlloc ();
-    psMetadataAdd (*header, PS_LIST_HEAD, "COMMENT", PS_META_MULTI, "folder for comment", NULL);
-    psMetadataAdd (*header, PS_LIST_HEAD, "HISTORY", PS_META_MULTI, "folder for history", NULL);
-    *header = psFitsReadHeader (*header, file);
-    *image  = psFitsReadImage  (*image,  file, region, 0);
+    psImageData *imdata = psAlloc(sizeof(psImageData));
+    psMetadata *header = psMetadataAlloc ();
+    psMetadataAdd (header, PS_LIST_HEAD, "COMMENT", PS_META_MULTI, "folder for comment", NULL);
+    psMetadataAdd (header, PS_LIST_HEAD, "HISTORY", PS_META_MULTI, "folder for history", NULL);
+    header = psFitsReadHeader (header, file);
+    psImage *image = psFitsReadImage (NULL, file, region, 0);
     psFree (file);
 
     // grab these values from the approrpiate location (image header if necessary)
     bool  status   = false;
-    float RDNOISE  = pmConfigLookupF32 (&status, *config, *header, "RDNOISE");
-    float GAIN     = pmConfigLookupF32 (&status, *config, *header, "GAIN");
-    psMetadataAdd (*config, PS_LIST_TAIL, "RDNOISE", PS_META_F32 | PS_META_REPLACE, "read noise value used", RDNOISE);
-    psMetadataAdd (*config, PS_LIST_TAIL, "GAIN",    PS_META_F32 | PS_META_REPLACE, "gain value used", GAIN);
+    float RDNOISE  = pmConfigLookupF32 (&status, config, header, "RDNOISE");
+    float GAIN     = pmConfigLookupF32 (&status, config, header, "GAIN");
+    psMetadataAdd (config, PS_LIST_TAIL, "RDNOISE", PS_META_F32 | PS_META_REPLACE, "read noise value used", RDNOISE);
+    psMetadataAdd (config, PS_LIST_TAIL, "GAIN",    PS_META_F32 | PS_META_REPLACE, "gain value used", GAIN);
 
-    // load the noise image if it exists, otherwise build from input
-    char *noiseName = pmConfigLookupSTR (&status, *config, *header, "NOISE");
+    imdata->image = image;
+    imdata->header = header;
+
+    // load the noise image if it is supplied, otherwise build from input
+    char *noiseName = psMetadataLookupSTR (&status, config, "NOISE");
     if (status == true) {
       file = psFitsAlloc (noiseName);
-      *noise  = psFitsReadImage  (*noise, file, region, 0);
+      imdata->noise = psFitsReadImage  (NULL, file, region, 0);
       psFree (file);
     } else {
-      *value = psScalar (1.0 / GAIN);
-      *noise = BinaryOp (NULL, *image, '/', value);
+      psScalar *value = psScalarAlloc (1.0 / GAIN, PS_TYPE_F64);
+      imdata->noise = BinaryOp (NULL, imdata->image, '/', value);
       psFree (value);
 
-      *value = psScalar (PS_SQR(RDNOISE/GAIN));
-      *noise = BinaryOp (*noise, *noise, '+', value);
+      psScalar *value = psScalarAlloc (PS_SQR(RDNOISE/GAIN), PS_TYPE_F64);
+      imdata->noise = BinaryOp (imdata->noise, imdata->noise, '+', value);
       psFree (value);
     }
 
-    char *maskName = pmConfigLookupSTR (&status, *config, *header, "MASK");
+    char *maskName = psMetadataLookupSTR (&status, config, "MASK");
     if (status == true) {
       file = psFitsAlloc (maskName);
-      *mask  = psFitsReadImage  (*mask, file, region, 0);
+      imdata->mask  = psFitsReadImage  (NULL, file, region, 0);
+      // check if U8 & force?
       psFree (file);
     } else {
-      *mask = psImageAlloc (image->numCols, image->numRows, PS_TYPE_U8);
+      imdata->mask = psImageAlloc (image->numCols, image->numRows, PS_TYPE_U8);
     }
+
+    float XBORDER  = psMetadataLookupF32 (&status, config, "XBORDER");
+    float YBORDER  = psMetadataLookupF32 (&status, config, "YBORDER");
+    psRegion *keep = psRegionAlloc (XBORDER, -XBORDER, YBORDER, -YBORDER);
+    keep           = psRegionForImage (keep, image, keep);
+    psImageKeepRegion (imdata->mask, keep, 0x01);
 
     psLogMsg ("psphot", 3, "load data: %f sec\n", psTimerMark ("psphot"));

trunk/psphot/src/source_moments.c

@@ -1 +1 @@
 # include "psphot.h"
 
-psArray *source_moments (psImage *image, psMetadata *config, psArray *allpeaks) 
+psArray *source_moments (psImageData *imdata, psMetadata *config, psArray *allpeaks) 
 {
     bool     status  = false;
@@ -8 +8 @@
 
     psTimerStart ("psphot");
-
-    // reduce the initial peaks to a subset of the better peaks
-    // is this the appropriate place for this culling?
-
-    float XBORDER  = psMetadataLookupF32 (&status, config, "XBORDER");
-    float YBORDER  = psMetadataLookupF32 (&status, config, "YBORDER");
-    float SATURATE = psMetadataLookupF32 (&status, config, "SATURATE");
-
-    psRegion *keep = psRegionAlloc (XBORDER, -XBORDER, YBORDER, -YBORDER);
-    keep           = psRegionForImage (keep, image, keep);
-    // I need to exclude the saturated objects here, but include them again later...?
-    peaks          = pmPeaksSubset (allpeaks, 100000, keep);
 
     // determine properties (sky, moments) of initial sources
@@ -32 +20 @@
     for (int i = 0; i < peaks->n; i++) {
         // provide the user arguments as metadata?
-        psSource *source  = pmSourceLocalSky (image, (psPeak *)peaks->data[i], PS_STAT_SAMPLE_MEDIAN, INNER, OUTER);
+        psSource *source  = pmSourceDefinePixels (imdata, (psPeak *)peaks->data[i], OUTER);
         if (source == NULL) continue;
-        pmSourceMaskRegion (source, keep);
+
+        pmSourceLocalSky_EAM (source, PS_STAT_SAMPLE_MEDIAN, INNER);
         pmSourceMoments (source, RADIUS);
         psArrayAdd (sources, 100, source);