Changeset 35855 for tags/ipp-20130712

Timestamp:

Jul 25, 2013, 3:21:58 PM (13 years ago)

Author:

watersc1

Message:

port from trunk of improved LAP logic and ppSkycell/coverage map code

Location:

tags/ipp-20130712

Files:

• ippScripts/scripts/lap_science.pl (modified) (5 diffs, 1 prop)
• ippScripts/scripts/skycell_jpeg.pl (modified) (8 diffs, 1 prop)
• ppSkycell/src (modified) (1 prop)
• ppSkycell/src/ppSkycell.h (modified) (1 diff)
• ppSkycell/src/ppSkycellArguments.c (modified) (2 diffs)
• ppSkycell/src/ppSkycellCamera.c (modified) (2 diffs)
• ppSkycell/src/ppSkycellData.c (modified) (2 diffs)
• ppSkycell/src/ppSkycellLoop.c (modified) (14 diffs)

tags/ipp-20130712/ippScripts/scripts/lap_science.pl

@@ -309 +309 @@
     my $magicDS_update_cmd  = "magicdstool -setfiletoupdate -magic_ds_id $chip_magicDS_id -set_label $label";
     $magicDS_update_cmd    .= " -dbname $dbname " if defined $dbname;
+    my $warptool_preupdate_cmd = "warptool -revertwarped -warp_id $warp_id -fault 26";
+    $warptool_preupdate_cmd .= " -dbname $dbname " if defined $dbname;
     my $warptool_update_cmd = "warptool -setskyfiletoupdate -warp_id $warp_id -set_label $label";
     $warptool_update_cmd   .= " -dbname $dbname " if defined $dbname;
@@ -314 +316 @@
     # only need to update the data if the warps are not full
     if ($warp->{state} ne 'full') { 
-        if ($chip->{state} eq 'cleaned') {
+        if (($chip->{state} eq 'cleaned')||  # Exposure is cleaned
+            (($chip->{state} eq 'update')&&($chip->{label} ne $label))) {  # Exposure is owned by an interloper.
             ($success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
                 run(command => $chiptool_update_cmd, verbose => $verbose);
@@ -341 +344 @@
         }
 
+
+        ($success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
+            run(command => $warptool_preupdate_cmd, verbose => $verbose);
+        unless ($success) {
+            $error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR);
+            &my_die("Unable to perform warptool --revertwarped: $error_code", $exposure->{exp_id}, $data_group);
+        }
+
         ($success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
             run(command => $warptool_update_cmd, verbose => $verbose);
@@ -364 +375 @@
     if (S64_IS_NOT_NULL($exposure->{chip_id})) {
         if (($exposure->{chip_state} eq 'cleaned')||
+            ($exposure->{chip_state} eq 'update')||
             ($exposure->{chip_state} eq 'goto_cleaned')||
             ($exposure->{chip_state} eq 'error_cleaned')) {
@@ -1407 +1419 @@
         'difftool -dbname @DBNAME@ -updaterun -set_state goto_cleaned -set_label goto_cleaned -label @LABEL@ -diff_id @DIFF_ID@',
         'magictool -dbname @DBNAME@ -updaterun -set_state full -set_label @LABEL@.old -magic_id @MAGIC_ID@',
-        'magicdstool -dbname @DBNAME@ -updaterun -set_state goto_cleaned -state full -set_label goto_cleaned -stage chip -stage_id @CHIP_ID@',
-        'magicdstool -dbname @DBNAME@ -updaterun -set_state goto_cleaned -state full -set_label goto_cleaned -stage warp -stage_id @WARP_ID@',
-        'magicdstool -dbname @DBNAME@ -updaterun -set_state goto_cleaned -state full -set_label goto_cleaned -stage diff -stage_id @DIFF_ID@');
+#       'magicdstool -dbname @DBNAME@ -updaterun -set_state goto_cleaned -state full -set_label goto_cleaned -stage chip -stage_id @CHIP_ID@',
+#       'magicdstool -dbname @DBNAME@ -updaterun -set_state goto_cleaned -state full -set_label goto_cleaned -stage warp -stage_id @WARP_ID@',
+#       'magicdstool -dbname @DBNAME@ -updaterun -set_state goto_cleaned -state full -set_label goto_cleaned -stage diff -stage_id @DIFF_ID@');
     
     foreach my $exposure (@$exposures) {

tags/ipp-20130712/ippScripts/scripts/skycell_jpeg.pl

@@ -17 +17 @@
 my ($help,$masks);
 my ($tempFile,$tempName,$maskFile,$maskName);
+$verbose = 0;
 my $ppSkycell = can_run('ppSkycell') or (warn "Can't find ppSkycell" and $missing_tools = 1);
 my $warptool  = can_run('warptool') or (warn "Can't find warptool" and $missing_tools = 1);
@@ -38 +39 @@
 pod2usage( -msg => "Unknown option: @ARGV", -exitval => 2) if @ARGV;
 pod2usage(
-    -msg => "Required options: --stage --stage_id --path_base\nHelpful options: --camera --dbname",
+    -msg => "Required options: --stage --stage_id --outroot\nHelpful options: --camera --dbname",
     -exitval => 3,
     ) unless
@@ -100 +101 @@
             }           
         }
-        print "$skycell_id $projection_cell\n"; 
+#       print "$skycell_id $projection_cell\n"; 
         my $file = $ipprc->filename("PSWARP.OUTPUT", $path_base, $skycell_id);
-        print "$file $state $quality\n";
+#       print "$file $state $quality\n";
         my $f_fh = $tangents{$projection_cell}{FILE};
         print $f_fh "$file\n";
         if ($masks) {
             my $mask = $ipprc->filename("PSWARP.OUTPUT.MASK", $path_base, $skycell_id);
-            print "$mask\n";
+#           print "$mask\n";
             my $m_fh = $tangents{$projection_cell}{MFILE};
             print $m_fh "$mask\n";
@@ -118 +119 @@
         }
         $command .= " ${outroot}.${projection_cell} ";
-        print "$command\n";
+#       print "$command\n";
         ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) = run(command => $command, verbose => $verbose);
         unless ($success) {
@@ -197 +198 @@
             }           
         }
-        print "$skycell_id $projection_cell\n"; 
+#       print "$skycell_id $projection_cell\n"; 
         my $file = $ipprc->filename("PPSUB.OUTPUT", $path_base, $skycell_id);
-        print "$file $state $quality\n";
+#       print "$file $state $quality\n";
         my $f_fh = $tangents{$projection_cell}{FILE};
         print $f_fh "$file\n";
         if ($masks) {
             my $mask = $ipprc->filename("PPSUB.OUTPUT.MASK", $path_base, $skycell_id);
-            print "$mask\n";
+#           print "$mask\n";
             my $m_fh = $tangents{$projection_cell}{MFILE};
             print $m_fh "$mask\n";
@@ -215 +216 @@
         }
         $command .= " ${outroot}.${projection_cell} ";
-        print "$command\n";
+#       print "$command\n";
         ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) = run(command => $command, verbose => $verbose);
         unless ($success) {
@@ -297 +298 @@
                 $tangents{$projection_cell}{$key}{FILE} = $tempFile;
                 $tangents{$projection_cell}{$key}{NAME} = $tempName;
+                $tangents{$projection_cell}{$key}{N}    = 0;
             }
         }
         foreach my $key (keys %products) {
-            print "$skycell_id $projection_cell\n";     
+#           print "$skycell_id $projection_cell\n";     
             my $file = $ipprc->filename($products{$key}, $path_base, $skycell_id);
-            print "$file $state $quality\n";
+#           print "$file $state $quality\n";
             my $f_fh = $tangents{$projection_cell}{$key}{FILE};
-            print $f_fh "$file\n";
+            if ($ipprc->file_exists($file)) {
+                print $f_fh "$file\n";
+                $tangents{$projection_cell}{$key}{N} ++;
+            }
         }
     }
@@ -322 +327 @@
                 $command .= " -Di BIN1 1 -Di BIN2 1 ";
             }
-            elsif ($key eq 'image') {
-                $command .= " -masks $tangents{$projection_cell}{mask}{NAME} ";
-            }
-            elsif ($key eq 'mask') { 
-                next; # This should be made with the images.
-            }
-            print "$command\n";
-            ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) = run(command => $command, verbose => $verbose);
-            unless ($success) {
-                $error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR);
-                &my_die("unable to perform ppSkycell: $error_code", $stage_id, $error_code);
+#           elsif ($key eq 'image') {
+#               $command .= " -masks $tangents{$projection_cell}{mask}{NAME} ";
+#           }
+#           elsif ($key eq 'mask') { 
+#               next; # This should be made with the images.
+#           }
+            else {
+                # Append the image list to other objects, in case the WCS information is unpopulated
+                $command .= " -wcsref $tangents{$projection_cell}{image}{NAME} ";
+            }
+            if ($tangents{$projection_cell}{$key}{N} > 0) {
+                print "$command\n";
+                ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) = run(command => $command, verbose => $verbose);
+                unless ($success) {
+                    $error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR);
+#                   print "Was unable to perform ppSkycell $stage_id $error_code $key\n";
+                    &my_die("unable to perform ppSkycell: $error_code", $stage_id, $error_code);
+                }
             }
         }

tags/ipp-20130712/ppSkycell/src/ppSkycell.h

@@ -11 +11 @@
     psString imagesName;                // Filename with images
     psString masksName;                 // Filename with masks
+    psString wcsrefName;                // Filename with WCS references
     psString outRoot;                   // Output root name
     int numInputs;                      // Number of inputs

tags/ipp-20130712/ppSkycell/src/ppSkycellArguments.c

@@ -51 +51 @@
     psMetadataAddStr(arguments, PS_LIST_TAIL, "-images", 0, "Filename with input images", NULL);
     psMetadataAddStr(arguments, PS_LIST_TAIL, "-masks", 0, "Filename with input masks", NULL);
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "-wcsref", 0, "Filename with WCS references", NULL);
     if (argc == 1 || !psArgumentParse(arguments, &argc, argv) || argc != 2) {
         usage(argv[0], arguments, data);
@@ -57 +58 @@
     data->imagesName = psMemIncrRefCounter(psMetadataLookupStr(NULL, arguments, "-images"));
     data->masksName = psMemIncrRefCounter(psMetadataLookupStr(NULL, arguments, "-masks"));
+    data->wcsrefName = psMemIncrRefCounter(psMetadataLookupStr(NULL, arguments, "-wcsref"));
+    
     data->outRoot = psStringCopy(argv[1]);
 

tags/ipp-20130712/ppSkycell/src/ppSkycellCamera.c

@@ -73 +73 @@
     }
 
+    psArray *wcsref = NULL;             // Names of WCS reference images
+    if (data->wcsrefName) {
+      wcsref = fileList(data->wcsrefName);
+      if (!wcsref) {
+        psError(psErrorCodeLast(), false, "No WCSrefs provided.");
+        psFree(images);
+        psFree(masks);
+        return(false);
+      }
+      if (wcsref->n != data->numInputs) {
+        psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Number of images (%ld) and wcsrefs (%ld) do not match",
+                images->n, wcsref->n);
+        psFree(images);
+        psFree(masks);
+        psFree(wcsref);
+        return(false);
+      }
+    }
+    
     psMetadataAddStr(data->config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "Output root", data->outRoot);
 
@@ -93 +112 @@
             }
         }
+
+        if (data->wcsrefName) {
+          fileArguments("WCSREF", wcsref->data[i], "Name of the WCS reference", data->config);
+          pmFPAfile *wcsref = pmFPAfileDefineFromArgs(&status, data->config, "PPSKYCELL.WCSREF", "WCSREF");
+          if (!status || !wcsref) {
+            psError(PS_ERR_IO, false, "Failed to build file from PPSKYCELL.WCSREF");
+            return(false);
+          }
+/*        if (!pmFPAfileBindFromArgs(&status, image, data->config, "PPSKYCELL.WCSREF", "WCSREF") || !status) { */
+/*          psError(PS_ERR_IO, false, "Failed to build file from PPSKYCELL.WCSREF"); */
+/*          return false; */
+/*        } */
+        }
+        
     }
 

tags/ipp-20130712/ppSkycell/src/ppSkycellData.c

@@ -15 +15 @@
     psFree(data->imagesName);
     psFree(data->masksName);
+    psFree(data->wcsrefName);
     psFree(data->outRoot);
     psFree(data->config);
@@ -28 +29 @@
     data->imagesName = NULL;
     data->masksName = NULL;
+    data->wcsrefName = NULL;
     data->outRoot = NULL;
     data->numInputs = 0;

tags/ipp-20130712/ppSkycell/src/ppSkycellLoop.c

@@ -10 +10 @@
 
 #define BUFFER 16                       // Size of buffer for projections
+
+#if (0)
+// This is all junk that is basically useless, because it calculates things we already know.
 
 static void regionMinMax(psRegion *base,// Base region; modified
@@ -22 +25 @@
 }
 
+
 static psRegion *skycellRegion(pmAstromWCS *wcs, // World Coordinate System
                                int numCols, int numRows // Size of image
@@ -34 +38 @@
     fromCoords->y = (Y); \
     psPlaneTransformApply(toCoords, wcs->trans, fromCoords); \
+    printf("%f %f %f %f\n",toCoords->x,toCoords->y,fromCoords->x,fromCoords->y); \
     toCoords->x /= wcs->cdelt1; \
     toCoords->y /= wcs->cdelt2; \
-    toCoords->x += wcs->crpix1; \
-    toCoords->y += wcs->crpix2; \
+    toCoords->x -= wcs->crpix1;         \
+    toCoords->y += wcs->crpix2;         \
     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);
-
+    //
     REGION_RANGE(0, 0);                 // Lower left
     psTrace("ppSkycell", 9, "Start: %.0f %.0f", toCoords->x, toCoords->y);
@@ -49 +54 @@
     REGION_RANGE(numCols, numRows);     // Upper right
     psTrace("ppSkycell", 9, "Stop: %.0f %.0f\n", toCoords->x, toCoords->y);
-
+    psTrace("ppSkycell", 7, "%g %g %g %g\n",wcs->cdelt1,wcs->cdelt2,wcs->crpix1,wcs->crpix2);
     return region;
 }
+#endif
 
 // Activate/deactivate a single element for a list
@@ -150 +156 @@
         pmAstromWCS *wcs = pmAstromWCSfromHeader(hdu->header); // World Coordinate System
         if (!wcs) {
+          if (data->wcsrefName) {
+            pmFPAfile *wcsref = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.WCSREF", i);
+            wcs = pmAstromWCSfromHeader(wcsref->fpa->hdu->header);
+          }
+          if (!wcs) {     
             psError(psErrorCodeLast(), false, "Unable to read WCS for image %d", i);
             return false;
-        }
-
-        psRegion *region = imageRegions->data[i] = skycellRegion(wcs, numCols, numRows); // Region of image
+          }
+        }
+
+        psRegion *region = psRegionAlloc(-1.0 * wcs->crpix1,numCols - wcs->crpix1,
+                                         -1.0 * wcs->crpix2,numRows - wcs->crpix2);
+        imageRegions->data[i] = region;
         psTrace("ppSkycell", 5, "Image region %d is: [%.0f:%.0f,%.0f:%.0f]\n",
                 i, region->x0, region->x1, region->y0, region->y1);
@@ -162 +176 @@
             if (wcs->crval1 == crval1->data.F64[j] && wcs->crval2 == crval2->data.F64[j] &&
                 wcs->cdelt1 == cdelt1->data.F64[j] && wcs->cdelt2 == cdelt1->data.F64[j]) {
-                regionMinMax(projRegions->data[j], region);
+              //                regionMinMax(projRegions->data[j], region);
+              psRegion *proj = projRegions->data[j];
+              proj->x0 = PS_MIN(region->x0,proj->x0);
+              proj->x1 = PS_MAX(region->x1,proj->x1);
+              proj->y0 = PS_MIN(region->y0,proj->y0);
+              proj->y1 = PS_MAX(region->y1,proj->y1);
                 target->data.S32[i] = j;
                 found = true;
@@ -188 +207 @@
     for (int i = 0; i < numProj; i++) {
         psRegion *projRegion = projRegions->data[i]; // Region for skycell projection
+        //      projRegion->x0 += 5760;
+        //projRegion->x1 += 5760;
+        //      projRegion->x0 += 
+        //      projRegion->x1 += 
         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
@@ -209 +232 @@
             psImageInit(mask2, 0xFF);
         }
+        // HDU containing the WCS we plan on using
         pmHDU *projhdu = NULL;
+
+        // Do we need to modify the WCS?  This flips to zero after we've done it once.
         int modify_wcs1 = 1;
         int modify_wcs2 = 1;
+
+        // Because we may have holes, we need to ensure that we set the CRPIX in the binned image correction.
+        // Find the minimum/maximum, so we know where the zero is.
+        float maxCRPIX1   = -99e99;
+        float maxCRPIX2   = -99e99;
         // Loop over inputs to this projection.
         for (int j = 0; j < data->numInputs; j++) {
@@ -229 +260 @@
             }
 
+            // Read the HDU/WCS information from the first entry, and use that as the reference.
             pmFPAfile *file = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.IMAGE", j);
+            pmAstromWCS *wcs = pmAstromWCSfromHeader(file->fpa->hdu->header); // World Coordinate System
+            if (!wcs) {
+              if (data->wcsrefName) {
+                pmFPAfile *wcsref = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.WCSREF", j);
+                wcs = pmAstromWCSfromHeader(wcsref->fpa->hdu->header);
+                if (!projhdu) {
+                  projhdu = wcsref->fpa->hdu;
+                }
+              }
+              if (!wcs) {         
+                psError(psErrorCodeLast(), false, "Unable to read WCS for image %d", j);
+                return false;
+              }
+            }
             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    
-            
+            }
+
+            // However, we need to check to see if we've found the maximum CRPIX
+            if (wcs->crpix1 > maxCRPIX1) {
+              maxCRPIX1 = wcs->crpix1;
+            }
+            if (wcs->crpix2 < maxCRPIX2) {
+              maxCRPIX2 = wcs->crpix2;
+            }
+
+            // Actuall do the binning.
             pmReadout *inRO = pmFPAviewThisReadout(view, file->fpa); // Readout with input
             psFree(view);
 
-            //      data->maskVal = 0xffff;
             pmReadout *bin1RO = pmReadoutAlloc(NULL), *bin2RO = pmReadoutAlloc(NULL); // Binned readouts
             if (!pmReadoutRebin(bin1RO, inRO, data->maskVal, data->bin1, data->bin1)) {
@@ -262 +304 @@
             }
 
+#if (0)
             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->y1 + projRegion->y1) / (float)data->bin1;
+            fprintf(stderr,"%d IM %f %f %f %f PR %f %f %f %f\n",
+                    j,
+                    imageRegion->x0,imageRegion->x1,
+                    imageRegion->y0,imageRegion->y1,
+                    projRegion->x0,projRegion->x1,
+                    projRegion->y0,projRegion->y1);
+#endif
+            // Offsets for image on this projection cell are just differences in CRPIX positions.
+            int xOffset1 = (-1 * wcs->crpix1 - projRegion->x0) / (float)(data->bin1);
+            int yOffset1 = (-1 * wcs->crpix2 - projRegion->y0) / (float)(data->bin1);
+            
             int xOffset2 = xOffset1 / (float)data->bin2, yOffset2 = yOffset1 / (float)data->bin2;
+            psTrace("ppSkycell",5,"Offsets: %d %d : %d %d",
+                    xOffset1,yOffset1,xOffset2,yOffset2);
+
+            // Overlay the data onto the appropriate pixels in the final outputs
             // XXX Completely neglecting rotations
             // The skycells are divided up neatly with them all having the same orientation
@@ -276 +331 @@
             }
 
+            // Cleanup on input loop.
             psFree(bin1RO);
             psFree(bin2RO);
@@ -345 +401 @@
             WCS->cdelt1 *= cd1f;
             WCS->cdelt2 *= cd2f;
-            WCS->crpix1 = WCS->crpix1 / cd1f;
-            WCS->crpix2 = WCS->crpix2 / cd2f;
+            // Fudge the CRPIX incase we have missing corners
+            if (maxCRPIX1 > WCS->crpix1) {
+              WCS->crpix1 = maxCRPIX1 / cd1f;
+            }
+            else {
+              WCS->crpix1 = WCS->crpix1 / cd1f;
+            }
+            if (maxCRPIX2 > WCS->crpix2) {
+              WCS->crpix2 = maxCRPIX2 / cd2f;
+            }
+            else {
+              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,
@@ -376 +443 @@
             }
             pmAstromWCStoHeader (fits1->fpa->hdu->header,WCS);
+#if WCS_DEBUG
             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,
@@ -448 +515 @@
             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++) {
+            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++) {
+            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);
               }