Changeset 28089 for trunk/psastro/src/psastroMaskUpdates.c

Timestamp:

May 25, 2010, 1:17:19 PM (16 years ago)

Author:

watersc1

Message:

Hopefully final edits for maskstats code

File:

• trunk/psastro/src/psastroMaskUpdates.c (modified) (11 diffs)

trunk/psastro/src/psastroMaskUpdates.c

@@ -30 +30 @@
 int colEnd(psU16 corner_list,int i,int x_0,int y_0,int R,int max) {
   if ((corner_list == 0x04)||(corner_list == 0x05)) {
-    return((int) fabs(x_0 - sqrt(pow(R,2) - pow(y_0 - i,2))));
+    int v = (int) fabs(x_0 - sqrt(pow(R,2) - pow(y_0 - i,2)));
+    if (v > max) {
+      return(max);
+    }
+    else {
+      return(v);
+    }
   }
   else {
@@ -54 +60 @@
            (corner_list == 0x03)||(corner_list == 0x0b)||
            (corner_list == 0x07)) {
-    return((int) fabs(y_0 - sqrt(pow(R,2) - pow(x_0 - j,2))));
+    int v = (int) fabs(y_0 - sqrt(pow(R,2) - pow(x_0 - j,2)));
+    if (v > max) {
+      return(max);
+    }
+    else {
+      return(v);
+    }
   }
   else {
@@ -62 +74 @@
 }
 
-/**
+/* #define MASK_DEBUG 1 */
+
+/*
  * create a mask or mask regions based on the collection of reference stars that * are in the vicinity of each chip
  */
@@ -112 +126 @@
         return false;
     }
-
     psLogMsg ("psastro", PS_LOG_INFO, "generating a bright-star mask");
-
+    
     bool REFSTAR_MASK_BLEED                = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK_BLEED");
 
@@ -137 +150 @@
     psU16 advisoryMaskVal = psMetadataLookupU32(&status, recipe, "MASKSTAT.ADVISORY");
 
-    psS32 Npix_valid = 0;
-    psS32 Npix_static = 0;
-    psS32 Npix_magic = 0;
-    psS32 Npix_dynamic = 0;
-    psS32 Npix_advisory = 0;
+    psS32 Npix_ref_valid = 0;
+    psS32 Npix_ref_static = 0;
+    psS32 Npix_ref_magic = 0;
+    psS32 Npix_ref_dynamic = 0;
+    psS32 Npix_ref_advisory = 0;
+
+    psS32 Npix_max_valid = 0;
+    psS32 Npix_max_static = 0;
+    psS32 Npix_max_magic = 0;
+    psS32 Npix_max_dynamic = 0;
+    psS32 Npix_max_advisory = 0;
     
     psU16 corner_list = 0x00;
@@ -177 +196 @@
 
     // Get camera specific mask stat options
-    psF32 FOV = psMetadataLookupF32(&status, fpaMask->camera, "FOV"); // Almost certainly doesn't exist. Placeholder for now.
-    psS32 CONSTANT_BLANK = psMetadataLookupS32(&status, fpaMask->camera, "NPIX_INTERCHIP"); // Need to add this as well.
-    
-    Npix_valid += CONSTANT_BLANK;
-    Npix_static += CONSTANT_BLANK;
-
+    psF32 FOV_REF = psMetadataLookupF32(&status, fpaMask->camera, "FOV_REF");
+    psF32 FOV_MAX = psMetadataLookupF32(&status, fpaMask->camera, "FOV_MAX");
+    psS32 NPIX_REF = psMetadataLookupS32(&status, fpaMask->camera, "NPIX_REF");
+    psS32 NPIX_MAX = psMetadataLookupS32(&status, fpaMask->camera, "NPIX_MAX");
     
     // select the reference mask fpa :: we use this to determine cell boundaries
@@ -209 +226 @@
     // open/load files as needed
     if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE;
-
+#ifdef MASK_DEBUG
+    psImage *masktest = psImageAlloc(10500,10500,PS_TYPE_U16);
+    psImageMaskType **maskIData = masktest->data.PS_TYPE_IMAGE_MASK_DATA;
+#endif
     // this loop selects the matched stars for all chips
     while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) {
@@ -435 +455 @@
           while ((cell = pmFPAviewNextCell(view, fpa, 1)) != NULL) {
             if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE;
-            psWarning("In CELL: FOV: %f CONSTANT: %d Mask bits: %d %d %d %d\n",
-                      FOV,CONSTANT_BLANK,staticMaskVal,magicMaskVal,dynamicMaskVal,advisoryMaskVal);
-            
             if (!cell->process || !cell->file_exists) {continue; }
             while ((readout = pmFPAviewNextReadout(view, fpa, 1)) != NULL) {
               pmReadout *readoutMask = pmFPAviewThisReadout (view, outMask->fpa);
               if (!pmFPAfileIOChecks (config, viewMask, PM_FPA_BEFORE)) ESCAPE;
-              psWarning("In READOUT: FOV: %f CONSTANT: %d Mask bits: %d %d %d %d\n",
-                        FOV,CONSTANT_BLANK,staticMaskVal,magicMaskVal,dynamicMaskVal,advisoryMaskVal);
-              
               if (!readoutMask->data_exists) {continue; }
+
               psPlane coordFPA;
               psPlane coordCell;
@@ -453 +468 @@
               psImageMaskType **maskData = mask->data.PS_TYPE_IMAGE_MASK_DATA;
               // Dance coordinates around
-              // Calculate which corners fall within the field of view.
+              // Calculate which corners fall within the field of view.  If this chip is fully contained, we can
+              // do a simple scan instead of checking it falls within the FOV.
               // 0x04   0x08
               // 0x01   0x02
+              corner_list = 0;
               
-              coordFPA.x = 0.0;
-              coordFPA.y = 0.0;
-              psPlaneTransformApply(&coordCell,chip->fromFPA,&coordFPA);
-              if (pow(coordFPA.x,2) + pow(coordFPA.y,2) <= pow(FOV,2)) {
+              coordCell.x = 0.0;
+              coordCell.y = 0.0;
+              psPlaneTransformApply(&coordFPA,chip->toFPA,&coordCell);
+              if (pow(coordFPA.x,2) + pow(coordFPA.y,2) <= pow(FOV_REF,2)) {
                 corner_list = corner_list | 0x01;
               }
-              
-              coordCell.x = mask->numCols - 1;
-              psPlaneTransformApply(&coordCell,chip->fromFPA,&coordFPA);
-              if (pow(coordFPA.x,2) + pow(coordFPA.y,2) <= pow(FOV,2)) {
+              coordCell.x = (1.0 * mask->numCols - 1.0);
+              psPlaneTransformApply(&coordFPA,chip->toFPA,&coordCell);
+              if (pow(coordFPA.x,2) + pow(coordFPA.y,2) <= pow(FOV_REF,2)) {
                 corner_list = corner_list | 0x02;
               }
-              
-              coordCell.y = mask->numRows - 1;
-              psPlaneTransformApply(&coordCell,chip->fromFPA,&coordFPA);
-              if (pow(coordFPA.x,2) + pow(coordFPA.y,2) <= pow(FOV,2)) {
+              coordCell.y = 1.0 * (mask->numRows - 1);
+              psPlaneTransformApply(&coordFPA,chip->toFPA,&coordCell);
+              if (pow(coordFPA.x,2) + pow(coordFPA.y,2) <= pow(FOV_REF,2)) {
                 corner_list = corner_list | 0x08;
               }
-              
-              coordCell.x = 0;
-              psPlaneTransformApply(&coordCell,chip->fromFPA,&coordFPA);
-              if (pow(coordFPA.x,2) + pow(coordFPA.y,2) <= pow(FOV,2)) {
+              coordCell.x = 0.0;
+              psPlaneTransformApply(&coordFPA,chip->toFPA,&coordCell);
+              if (pow(coordFPA.x,2) + pow(coordFPA.y,2) <= pow(FOV_REF,2)) {
                 corner_list = corner_list | 0x04;
               }
-              
-              int x_0 = coordFPA.x;
-              int y_0 = coordFPA.y;
-              
-              
+
               // Scan over the valid regions of the image and count masked pixels
-              if ((corner_list == 0x04)||(corner_list == 0x08)||
-                  (corner_list == 0x05)||(corner_list == 0x0a)) {
-                for (int i = 0; i < mask->numRows - 1; i++) {
-                  for (int j = colStart(corner_list,i,x_0,y_0,FOV,mask->numCols - 1);
-                       j < colEnd(corner_list,i,x_0,y_0,FOV,mask->numCols - 1); j++) {
-                    Npix_valid++;
-                    if (maskData[i][j] & staticMaskVal) {
-                      Npix_static++;
-                      continue;
-                    }
-                    if (maskData[i][j] & dynamicMaskVal) {
-                      Npix_dynamic++;
-                      continue;
-                    }
-                    if (maskData[i][j] & magicMaskVal) {
-                      Npix_magic++;
-                      continue;
-                    }
-                    if (maskData[i][j] & advisoryMaskVal) {
-                      Npix_advisory++;
-                      continue;
-                    }
+              for (int i = 0; i < mask->numRows - 1; i++) {
+                for (int j = 0; j < mask->numCols - 1; j++) {
+                  coordCell.x = j;
+                  coordCell.y = i;
+                  coordFPA.x = 0.0;
+                  coordFPA.y = 0.0;
+                  int region = 0;
+
+                  if (corner_list == 0x0f) {
+                    Npix_ref_valid++;
+                    Npix_max_valid++;
+                    region = 1;
+/* #ifdef MASK_DEBUG */
+/*                  psPlaneTransformApply(&coordFPA,chip->toFPA,&coordCell); */
+/*                  maskIData[(int) (coordFPA.y/4 + 5250)][(int) (coordFPA.x/4 + 5250)] = */
+/*                    maskIData[(int) (coordFPA.y/4 + 5250)][(int) (coordFPA.x/4 + 5250)] + 0x01; */
+/* #endif */
                   }
-                }
-              }
-              else {
-                for (int j = 0; j < mask->numCols - 1; j++) {
-                  for (int i = rowStart(corner_list,j,x_0,y_0,FOV,mask->numRows - 1);
-                       i < rowEnd(corner_list,j,x_0,y_0,FOV,mask->numRows - 1); i++) {
-                    Npix_valid++;
-                    if (maskData[i][j] & staticMaskVal) {
-                      Npix_static++;
-                      continue;
-                    }
-                    if (maskData[i][j] & dynamicMaskVal) {
-                      Npix_dynamic++;
-                      continue;
-                    }
-                    if (maskData[i][j] & magicMaskVal) {
-                      Npix_magic++;
-                      continue;
-                    }
-                    if (maskData[i][j] & advisoryMaskVal) {
-                      Npix_advisory++;
-                      continue;
-                    }
+                  if (!region) {
+                    psPlaneTransformApply(&coordFPA,chip->toFPA,&coordCell);
+                    if (pow(coordFPA.x,2) + pow(coordFPA.y,2) <= pow(FOV_REF,2)) {
+                      Npix_ref_valid++;
+                      Npix_max_valid++;
+                      region = 1;
+/* #ifdef MASK_DEBUG */
+/*                    maskIData[(int) (coordFPA.y/4 + 5250)][(int) (coordFPA.x/4 + 5250)] = */
+/*                      maskIData[(int) (coordFPA.y/4 + 5250)][(int) (coordFPA.x/4 + 5250)] +  0x01; */
+/* #endif */
+                    }
+                    else if (pow(coordFPA.x,2) + pow(coordFPA.y,2) <= pow(FOV_MAX,2)) {
+                      Npix_max_valid++;
+                      region = 2;
+/* #ifdef MASK_DEBUG  */
+/*                    maskIData[(int) (coordFPA.y/4 + 5250)][(int) (coordFPA.x/4 + 5250)] = */
+/*                      maskIData[(int) (coordFPA.y/4 + 5250)][(int) (coordFPA.x/4 + 5250)] + 0x01; */
+/* #endif */
+                    }
+                  }
+                  if (!region) {
+                    continue;
+                  }
+                  
+                  if (maskData[i][j] & staticMaskVal) {
+#ifdef MASK_DEBUG
+                    psPlaneTransformApply(&coordFPA,chip->toFPA,&coordCell);
+                    maskIData[(int) (coordFPA.y/4 + 5250)][(int) (coordFPA.x/4 + 5250)] =
+                      maskIData[(int) (coordFPA.y/4 + 5250)][(int) (coordFPA.x/4 + 5250)] + 1;
+#endif
+                    if (region == 1) {
+                      Npix_ref_static++;
+                      Npix_max_static++;
+                    }
+                    if (region == 2) {
+                      Npix_max_static++;
+                    }
+                    continue;
+                  }
+                  if (maskData[i][j] & dynamicMaskVal) {
+                    if (region == 1) {
+                      Npix_ref_dynamic++;
+                      Npix_max_dynamic++;
+                    }
+                    if (region == 2) {
+                      Npix_max_dynamic++;
+                    }
+                    continue;
+                  }
+                  if (maskData[i][j] & magicMaskVal) {
+                    if (region == 1) {
+                      Npix_ref_magic++;
+                      Npix_max_magic++;
+                    }
+                    if (region == 2) {
+                      Npix_max_magic++;
+                    }
+                    continue;
+                  }
+                  if (maskData[i][j] & advisoryMaskVal) {
+                    if (region == 1) {
+                      Npix_ref_advisory++;
+                      Npix_max_advisory++;
+                    }
+                    if (region == 2) {
+                      Npix_max_advisory++;
+                    }
+                    continue;
                   }
                 }
@@ -555 +605 @@
     if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE;
 
+#ifdef MASK_DEBUG
+    psFits *maskFits = psFitsOpen("/data/ipp007.0/watersc1/mask.test.fits","w");
+    psFitsWriteImage(maskFits,NULL,masktest,1,"mask");
+
+    psFree(maskFits);
+    psFree(masktest);
+#endif
     if (COUNT_GHOSTS) {
         // save nGhosts to update header.
@@ -566 +623 @@
     }
 
-    psMetadataAddS32(stats,PS_LIST_TAIL, "MASKFRAC_NPIX", 0,
-                     "Number of valid pixels", Npix_valid);
-    psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_STATIC", 0,
-                     "Fraction of pixels statically masked", (float) Npix_static / Npix_valid);
-    psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_DYNAMIC", 0,
-                     "Fraction of pixels dynamically masked", (float) Npix_dynamic / Npix_valid);
-    psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_MAGIC", 0,
-                     "Fraction of pixels magically masked", (float) Npix_magic / Npix_valid);
-    psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_ADVISORY", 0,
-                     "Fraction of pixels masked as an advisory", (float) Npix_advisory / Npix_valid);
+    Npix_ref_static += (NPIX_REF - Npix_ref_valid);
+    Npix_max_static += (NPIX_MAX - Npix_max_valid);
+    psMetadataAddS32(stats,PS_LIST_TAIL, "MASKFRAC_REF_NPIX", 0,
+                     "Number of valid pixels", Npix_ref_valid);
+    psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_REF_STATIC", 0,
+                     "Fraction of pixels statically masked", (float) Npix_ref_static / NPIX_REF);
+    psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_REF_DYNAMIC", 0,
+                     "Fraction of pixels dynamically masked", (float) Npix_ref_dynamic / NPIX_REF);
+    psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_REF_MAGIC", 0,
+                     "Fraction of pixels magically masked", (float) Npix_ref_magic / NPIX_REF);
+    psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_REF_ADVISORY", 0,
+                     "Fraction of pixels masked as an advisory", (float) Npix_ref_advisory / NPIX_REF);
+
+    psMetadataAddS32(stats,PS_LIST_TAIL, "MASKFRAC_MAX_NPIX", 0,
+                     "Number of valid pixels", Npix_max_valid);
+    psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_MAX_STATIC", 0,
+                     "Fraction of pixels statically masked", (float) Npix_max_static / NPIX_MAX);
+    psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_MAX_DYNAMIC", 0,
+                     "Fraction of pixels dynamically masked", (float) Npix_max_dynamic / NPIX_MAX);
+    psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_MAX_MAGIC", 0,
+                     "Fraction of pixels magically masked", (float) Npix_max_magic / NPIX_MAX);
+    psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_MAX_ADVISORY", 0,
+                     "Fraction of pixels masked as an advisory", (float) Npix_max_advisory / NPIX_MAX);
     
     // deactivate all files