Changeset 34800 for trunk/ppSkycell/src/ppSkycellLoop.c

Timestamp:

Dec 11, 2012, 2:04:31 PM (13 years ago)

Author:

watersc1

Message:

Merge from my branch including background restoration and stack stage projection cell binned images.

Location:

trunk

Files:

• . (modified) (1 prop)
• ppSkycell/src/ppSkycellLoop.c (modified) (12 diffs)

trunk/ppSkycell/src/ppSkycellLoop.c

@@ -15 +15 @@
                          )
 {
-    base->x0 = PS_MIN(base->x0, compare->x0);
-    base->x1 = PS_MAX(base->x1, compare->x1);
+    base->x0 = PS_MAX(base->x0, compare->x0);
+    base->x1 = PS_MIN(base->x1, compare->x1);
     base->y0 = PS_MIN(base->y0, compare->y0);
     base->y1 = PS_MAX(base->y1, compare->y1);
@@ -27 +27 @@
 {
     psPlane *fromCoords = psPlaneAlloc(), *toCoords = psPlaneAlloc(); // Coordinates for transforms
-    psRegion *region = psRegionAlloc(INFINITY, -INFINITY, INFINITY, -INFINITY); // Region for skycell
+    psRegion *region = psRegionAlloc(-INFINITY, INFINITY, INFINITY, -INFINITY); // Region for skycell
 
 // Limit the region using the nominated point (X,Y)
@@ -38 +38 @@
     toCoords->x += wcs->crpix1; \
     toCoords->y += wcs->crpix2; \
-    region->x0 = PS_MIN(region->x0, toCoords->x); \
-    region->x1 = PS_MAX(region->x1, toCoords->x); \
+    region->x0 = PS_MAX(region->x0, toCoords->x); \
+    region->x1 = PS_MIN(region->x1, toCoords->x); \
     region->y0 = PS_MIN(region->y0, toCoords->y); \
     region->y1 = PS_MAX(region->y1, toCoords->y);
@@ -132 +132 @@
     psVector *target = psVectorAlloc(data->numInputs, PS_TYPE_S32); // Target for each input
     psArray *imageRegions = psArrayAlloc(data->numInputs); // Region for image
-
+    psArray *regionHDUs = psArrayAlloc(data->numInputs);
+    // Determine which projection cells we have to deal with.
     for (int i = 0; i < data->numInputs; i++) {
         pmFPAfile *file = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.IMAGE", i); // File to examine
@@ -178 +179 @@
             target->data.S32[i] = numProj;
             psTrace("ppSkycell", 3, "Image %d uses new projection\n", i);
-
+            psArrayAdd(regionHDUs,1,hdu);
             numProj++;
         }
@@ -184 +185 @@
 
     pmFPAfileActivate(data->config->files, false, NULL);
-
+    // Loop over projections
     for (int i = 0; i < numProj; i++) {
         psRegion *projRegion = projRegions->data[i]; // Region for skycell projection
         psTrace("ppSkycell", 2, "Projection %d: [%.0f:%.0f,%.0f:%.0f]\n",
                 i, projRegion->x0, projRegion->x1, projRegion->y0, projRegion->y1);
-        int xSize = projRegion->x1 - projRegion->x0 + 1; // Size of unbinned image
+        int xSize = -projRegion->x1 + projRegion->x0 + 1; // Size of unbinned image
         int ySize = projRegion->y1 - projRegion->y0 + 1; // Size of unbinned image
         // Size of binned image 1
@@ -198 +199 @@
         psImage *image1 = psImageAlloc(numCols1, numRows1, PS_TYPE_F32); // Binned image
         psImage *image2 = psImageAlloc(numCols2, numRows2, PS_TYPE_F32); // Binned image
-        psImageInit(image1, 0);
-        psImageInit(image2, 0);
-
+        psImageInit(image1,NAN);
+        psImageInit(image2,NAN);
+        
         psImage *mask1 = NULL, *mask2 = NULL; // Binned masks
         if (data->masksName) {
@@ -208 +209 @@
             psImageInit(mask2, 0xFF);
         }
-
+        pmHDU *projhdu = NULL;
+        int modify_wcs1 = 1;
+        int modify_wcs2 = 1;
+        // Loop over inputs to this projection.
         for (int j = 0; j < data->numInputs; j++) {
             if (target->data.S32[j] != i) {
                 continue;
             }
-
             pmFPAfileActivateSingle(data->config->files, true, "PPSKYCELL.IMAGE", j);
             if (data->masksName) {
@@ -227 +230 @@
 
             pmFPAfile *file = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.IMAGE", j);
+            if (!projhdu) {
+              projhdu = file->fpa->hdu;
+              //              psMetadataCopy(projhdu->header,file->fpa->hdu->header);
+            }
+#if 0
+            else {
+              if ((psMetadataLookupF32(NULL,file->fpa->hdu->header,"CRPIX1") >=
+                   psMetadataLookupF32(NULL,projhdu->header,"CRPIX1"))&&
+                  (psMetadataLookupF32(NULL,file->fpa->hdu->header,"CRPIX2") >=
+                   psMetadataLookupF32(NULL,projhdu->header,"CRPIX1"))) {
+                projhdu = file->fpa->hdu;
+                psMetadataCopy(projhdu->header,file->fpa->hdu->header);
+              }
+            }
+#endif    
+            
             pmReadout *inRO = pmFPAviewThisReadout(view, file->fpa); // Readout with input
             psFree(view);
 
-            // Flip images; no idea why this has to be done, but apparently it does
-            {
-                psImage *rot = psImageRotate(NULL, inRO->image, M_PI, NAN, PS_INTERPOLATE_BILINEAR);
-                psFree(inRO->image);
-                inRO->image = rot;
-            }
-            if (inRO->mask) {
-                psImage *rot = psImageRotate(NULL, inRO->mask, M_PI, 0, PS_INTERPOLATE_FLAT);
-                psFree(inRO->mask);
-                inRO->mask = rot;
-            }
-
+            //      data->maskVal = 0xffff;
             pmReadout *bin1RO = pmReadoutAlloc(NULL), *bin2RO = pmReadoutAlloc(NULL); // Binned readouts
             if (!pmReadoutRebin(bin1RO, inRO, data->maskVal, data->bin1, data->bin1)) {
@@ -256 +264 @@
             psRegion *imageRegion = imageRegions->data[j]; // Region for image
             // Offsets for image on skycell
-            int xOffset1 = (imageRegion->x0 - projRegion->x0) / (float)data->bin1;
-            int yOffset1 = (imageRegion->y0 - projRegion->y0) / (float)data->bin1;
+            int xOffset1 = (-imageRegion->x0 + projRegion->x0) / (float)data->bin1;
+            int yOffset1 = (-imageRegion->y1 + projRegion->y1) / (float)data->bin1;
             int xOffset2 = xOffset1 / (float)data->bin2, yOffset2 = yOffset1 / (float)data->bin2;
-
             // XXX Completely neglecting rotations
             // The skycells are divided up neatly with them all having the same orientation
-            psImageOverlaySection(image1, bin1RO->image, xOffset1, yOffset1, "=");
-            psImageOverlaySection(image2, bin2RO->image, xOffset2, yOffset2, "=");
+            psImageOverlaySection(image1, bin1RO->image, xOffset1, yOffset1, "E");
+            psImageOverlaySection(image2, bin2RO->image, xOffset2, yOffset2, "E");
             if (data->masksName) {
-                psImageOverlaySection(mask1, bin1RO->mask, xOffset1, yOffset1, "=");
-                psImageOverlaySection(mask2, bin2RO->mask, xOffset2, yOffset2, "=");
+                psImageOverlaySection(mask1, bin1RO->mask, xOffset1, yOffset1, "M");
+                psImageOverlaySection(mask2, bin2RO->mask, xOffset2, yOffset2, "M");
             }
 
@@ -287 +294 @@
           pmFPAfileActivate(data->config->files, true, "PPSKYCELL.BIN1");
           pmFPAfileActivate(data->config->files, true, "PPSKYCELL.BIN2");
+          if (data->masksName) {
+            pmFPAfileActivate(data->config->files, true, "PPSKYCELL.BIN1.MASK");
+            pmFPAfileActivate(data->config->files, true, "PPSKYCELL.BIN2.MASK");
+          }
         }
         
@@ -313 +324 @@
 
         if (data->doFits) {
+
+          pmFPAfile *fits1 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.BIN1", 0);
+          // Copy header from projection root hdu
+          fits1->fpa->hdu = pmHDUAlloc(NULL);
+          fits1->fpa->hdu->header = psMetadataAlloc();
+          psMetadataCopy(fits1->fpa->hdu->header,projhdu->header);
+
+          // Change wcs here
+#define WCS_DEBUG 0
+          if (modify_wcs1) {
+            pmAstromWCS *WCS = pmAstromWCSfromHeader(fits1->fpa->hdu->header);
+            double cd1f = 1.0 * data->bin1;
+            double cd2f = 1.0 * data->bin1;
+#if WCS_DEBUG
+            fprintf(stderr,">>> %d %d (%g %g) (%g %g) (%g %g)\n",data->bin1,data->bin2,
+                    cd1f,cd2f,WCS->cdelt1,WCS->cdelt2,
+                    WCS->crpix1,WCS->crpix2
+                    );
+#endif
+            WCS->cdelt1 *= cd1f;
+            WCS->cdelt2 *= cd2f;
+            WCS->crpix1 = WCS->crpix1 / cd1f;
+            WCS->crpix2 = WCS->crpix2 / cd2f;
+#if WCS_DEBUG
+            fprintf(stderr,">>> %d %d (%g %g) (%g %g) (%g %g) %d %d %d %d\n",data->bin1,data->bin2,
+                    cd1f,cd2f,WCS->cdelt1,WCS->cdelt2,
+                    WCS->crpix1,WCS->crpix2,
+                    WCS->trans->x->nX,              WCS->trans->x->nY,
+                    WCS->trans->y->nX,              WCS->trans->y->nY
+                    );
+#endif
+            for (int q = 0; q <= WCS->trans->x->nX; q++) {
+              for (int r = 0; r <= WCS->trans->x->nY; r++) {
+                WCS->trans->x->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r);
+#if WCS_DEBUG
+                fprintf(stderr,"PC: %d %d %g %g\n",
+                        q,r,pow(cd1f,q) * pow(cd2f,r),
+                        WCS->trans->x->coeff[q][r]);
+#endif
+              }
+            }
+            for (int q = 0; q <= WCS->trans->y->nX; q++) {
+              for (int r = 0; r <= WCS->trans->y->nY; r++) {
+                WCS->trans->y->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r);
+#if WCS_DEBUG
+                fprintf(stderr,"PC: %d %d %g %g\n",
+                         q,r,pow(cd1f,q) * pow(cd2f,r),
+                         WCS->trans->y->coeff[q][r]);
+#endif
+              }
+            }
+            pmAstromWCStoHeader (fits1->fpa->hdu->header,WCS);
+            WCS = pmAstromWCSfromHeader(fits1->fpa->hdu->header);
+#if WCS_DEBUG
+            fprintf(stderr,">>> %d %d (%g %g) (%g %g) (%g %g)\n",data->bin1,data->bin2,
+                    cd1f,cd2f,WCS->cdelt1,WCS->cdelt2,
+                    WCS->crpix1,WCS->crpix2
+                    );
+#endif
+            modify_wcs1 = 0;
+          }
+
+          
+          pmChip *Fchip1 = pmFPAfileThisChip(data->config->files, view, "PPSKYCELL.BIN1");
+          psMetadataAddS32(Fchip1->concepts,PS_LIST_TAIL,"CHIP.XPARITY", PS_META_REPLACE,"",1);
+          psMetadataAddS32(Fchip1->concepts,PS_LIST_TAIL,"CHIP.YPARITY", PS_META_REPLACE,"",1);
+
           pmCell *Fcell1 = pmFPAfileThisCell(data->config->files, view, "PPSKYCELL.BIN1"); // Rebinned cell 1
-          pmCell *Fcell2 = pmFPAfileThisCell(data->config->files, view, "PPSKYCELL.BIN2"); // Rebinned cell 2
-
-          // This is a hack to get a functioning header created so the fits images can be written out.
+/*        // This is a hack to get a functioning header created so the fits images can be written out. */
           psMetadataAddS32(Fcell1->concepts,PS_LIST_TAIL,"CELL.XPARITY", PS_META_REPLACE,"",1);
           psMetadataAddS32(Fcell1->concepts,PS_LIST_TAIL,"CELL.YPARITY", PS_META_REPLACE,"",1);
           psMetadataAddS32(Fcell1->concepts,PS_LIST_TAIL,"CELL.READDIR", PS_META_REPLACE,"",1);
+
+          pmReadout *Fro1 = pmReadoutAlloc(Fcell1);
+          Fro1->image = image1;
+          Fro1->mask = mask1;
+          Fro1->data_exists = Fcell1->data_exists = Fcell1->parent->data_exists = true;
+          
+          fits1->save = true;
+          fits1->fileIndex = i;
+
+          // Do the mask if we need to
+          if (data->masksName) {
+/*          pmFPAfile *Mask1 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.BIN1.MASK", 0); */
+/*          //      Mask1->fpa->hdu = pmHDUAlloc(NULL); */
+/*          //      Mask1->fpa->hdu->header = psMetadataAlloc(); */
+/*          //      psMetadataCopy(Mask1->fpa->hdu->header,fits1->fpa->hdu->header); */
+/*          pmChip *Mchip1 = pmFPAfileThisChip(data->config->files, view, "PPSKYCELL.BIN1.MASK"); */
+/*          psMetadataAddS32(Mchip1->concepts,PS_LIST_TAIL,"CHIP.XPARITY", PS_META_REPLACE,"",1); */
+/*          psMetadataAddS32(Mchip1->concepts,PS_LIST_TAIL,"CHIP.YPARITY", PS_META_REPLACE,"",1); */
+            
+/*          pmCell *Mcell1 = pmFPAfileThisCell(data->config->files, view, "PPSKYCELL.BIN1.MASK"); // Rebinned cell 1 */
+/*          /\*           // This is a hack to get a functioning header created so the fits images can be written out. *\/ */
+/*          psMetadataAddS32(Mcell1->concepts,PS_LIST_TAIL,"CELL.XPARITY", PS_META_REPLACE,"",1); */
+/*          psMetadataAddS32(Mcell1->concepts,PS_LIST_TAIL,"CELL.YPARITY", PS_META_REPLACE,"",1); */
+/*          psMetadataAddS32(Mcell1->concepts,PS_LIST_TAIL,"CELL.READDIR", PS_META_REPLACE,"",1); */
+            
+/*          pmReadout *Mro1 = pmReadoutAlloc(Mcell1); */
+/*          Mro1->image = image1; */
+/*          Mro1->mask = mask1; */
+/*          Mro1->data_exists = Mcell1->data_exists = Mcell1->parent->data_exists = true; */
+            
+/*          Mask1->save = true; */
+/*          Mask1->fileIndex = i; */
+          }
+            
+            
+          
+          // Repeat with second binned image
+          pmFPAfile *fits2 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.BIN2", 0);
+          fits2->fpa->hdu = pmHDUAlloc(NULL);
+          fits2->fpa->hdu->header = psMetadataAlloc();
+          psMetadataCopy(fits2->fpa->hdu->header,projhdu->header);
+
+          if (modify_wcs2) {
+            pmAstromWCS *WCS = pmAstromWCSfromHeader(fits2->fpa->hdu->header);
+            double cd1f = 1.0 * data->bin2 * data->bin1;
+            double cd2f = 1.0 * data->bin2 * data->bin1;
+            
+            WCS->cdelt1 *= cd1f;
+            WCS->cdelt2 *= cd2f;
+            WCS->crpix1 = WCS->crpix1 / cd1f;
+            WCS->crpix2 = WCS->crpix2 / cd2f;
+            
+            for (int q = 0; q < WCS->trans->x->nX; q++) {
+              for (int r = 0; r < WCS->trans->x->nY; r++) {
+                WCS->trans->x->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r);
+              }
+            }
+            for (int q = 0; q < WCS->trans->y->nX; q++) {
+              for (int r = 0; r < WCS->trans->y->nY; r++) {
+                WCS->trans->y->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r);
+              }
+            }
+            pmAstromWCStoHeader (fits2->fpa->hdu->header,WCS);
+            modify_wcs2 = 0;
+          }
+          
+          pmChip *Fchip2 = pmFPAfileThisChip(data->config->files, view, "PPSKYCELL.BIN2");
+          psMetadataAddS32(Fchip2->concepts,PS_LIST_TAIL,"CHIP.XPARITY", PS_META_REPLACE,"",1);
+          psMetadataAddS32(Fchip2->concepts,PS_LIST_TAIL,"CHIP.YPARITY", PS_META_REPLACE,"",1);
+
+          
+          pmCell *Fcell2 = pmFPAfileThisCell(data->config->files, view, "PPSKYCELL.BIN2"); // Rebinned cell 2
           psMetadataAddS32(Fcell2->concepts,PS_LIST_TAIL,"CELL.XPARITY", PS_META_REPLACE,"",1);
           psMetadataAddS32(Fcell2->concepts,PS_LIST_TAIL,"CELL.YPARITY", PS_META_REPLACE,"",1);
           psMetadataAddS32(Fcell2->concepts,PS_LIST_TAIL,"CELL.READDIR", PS_META_REPLACE,"",1);
-
-          psMetadataAddS32(Fcell1->parent->concepts,PS_LIST_TAIL,"CHIP.XPARITY", PS_META_REPLACE,"",1);
-          psMetadataAddS32(Fcell1->parent->concepts,PS_LIST_TAIL,"CHIP.YPARITY", PS_META_REPLACE,"",1);
-          psMetadataAddS32(Fcell2->parent->concepts,PS_LIST_TAIL,"CHIP.XPARITY", PS_META_REPLACE,"",1);
-          psMetadataAddS32(Fcell2->parent->concepts,PS_LIST_TAIL,"CHIP.YPARITY", PS_META_REPLACE,"",1);
-
-          pmReadout *Fro1 = pmReadoutAlloc(Fcell1), *Fro2 = pmReadoutAlloc(Fcell2); // Binned readouts
-          
-          Fro1->image = image1;
+          
+          pmReadout *Fro2 = pmReadoutAlloc(Fcell2); 
           Fro2->image = image2;
-          Fro1->mask = mask1;
           Fro2->mask = mask2;
-          
-          Fro1->data_exists = Fcell1->data_exists = Fcell1->parent->data_exists = true;
           Fro2->data_exists = Fcell2->data_exists = Fcell2->parent->data_exists = true;
-          
-          pmFPAfile *fits1 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.BIN1", 0);
-          fits1->save = true;
-          fits1->fileIndex = i;
-          pmFPAfile *fits2 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.BIN2", 0);
+
           fits2->save = true;
           fits2->fileIndex = i;
+
+          
         }