Changeset 24691

Timestamp:

Jul 6, 2009, 4:24:20 PM (17 years ago)

Author:

bills

Message:

when destreaking warp and diff skyfiles, censor any sources whose coordinates are masked

Location:

trunk

Files:

• ippScripts/scripts/magic_destreak.pl (modified) (4 diffs)
• magic/remove/src/streaksio.c (modified) (7 diffs)
• magic/remove/src/streaksremove.c (modified) (6 diffs)
• magic/remove/src/streaksremove.h (modified) (1 diff)

trunk/ippScripts/scripts/magic_destreak.pl

@@ -194 +194 @@
     }
 
-    my ($image, $mask, $ch_mask, $weight, $astrom);
+    my ($image, $mask, $ch_mask, $weight, $astrom, $sources);
 
     if ($stage eq "raw") {
@@ -211 +211 @@
         $mask   = $ipprc->filename("PSWARP.OUTPUT.MASK", $path_base);
         $weight = $ipprc->filename("PSWARP.OUTPUT.VARIANCE", $path_base);
+        $sources    = $ipprc->filename("PSWARP.OUTPUT.SOURCES", $path_base);
     } elsif ($stage eq "diff") {
         my $name = $inverse ? "PPSUB.INVERSE" : "PPSUB.OUTPUT"; # Base name for images
@@ -216 +217 @@
         $mask   = $ipprc->filename("$name.MASK", $path_base);
         $weight = $ipprc->filename("$name.VARIANCE", $path_base);
+        $sources    = $ipprc->filename("PPSUB.OUTPUT.SOURCES", $path_base);
     }
 
@@ -227 +229 @@
         $command .= " -chip_mask $ch_mask" if defined $ch_mask;
         $command .= " -weight $weight" if defined $weight;
+        $command .= " -sources $sources" if defined $sources;
         $command .= " -skycelllist $skycell_list" if defined $skycell_list;
         $command .= " -replace" if $replace;

trunk/magic/remove/src/streaksio.c

@@ -20 +20 @@
     memset(sf, 0, sizeof(*sf));
 
+    sf->config = config;
+    sf->program_name = basename(program_name);
+
     if (remove) {
         // remember pointer so that streaksExit can delete temps
@@ -25 +28 @@
     }
 
-    sf->config = config;
-    sf->program_name = basename(program_name);
 
     // error checking is done by sFileOpen. If a file can't be opened we just exit
@@ -80 +81 @@
         } else {
             sf->recWeight = sFileOpen(config, stage, "RECOVERY.WEIGHT", NULL, true);
+        }
+    }
+    if (remove) {
+        sf->inSources = sFileOpen(config, stage, "SOURCES", NULL, false);
+        if (sf->inSources) {
+            inputBasename = basename(sf->inSources->name);
+            sf->outSources = sFileOpen(config, stage, "OUTPUT", inputBasename, true);
         }
     }
@@ -607 +615 @@
         streaksExit("", PS_EXIT_DATA_ERROR);
     }
-
+ 
     bool status;
     in->extname = psMetadataLookupStr(&status, in->header, "EXTNAME");
@@ -929 +937 @@
         if (!replicate(sfiles->outWeight, sfiles->inWeight)) {
             psError(PM_ERR_SYS, false, "failed to replicate outImage.");
+            return false;
+        }
+    }
+
+    if (sfiles->outSources) {
+        if (!replicate(sfiles->outSources, sfiles->inSources)) {
+            psError(PM_ERR_SYS, false, "failed to replicate outSources.");
             return false;
         }
@@ -997 +1012 @@
     }
 
+    if (sfiles->outSources) {
+        if (!swapOutputToInput(sfiles->inSources, sfiles->outSources)) {
+            psError(PM_ERR_SYS, false, "failed to swap instances for sources.");
+            return false;
+        }
+    }
+
     if (!swapOutputToInput(sfiles->inImage, sfiles->outImage)) {
         psError(PM_ERR_SYS, false, "failed to swap instances for Image.");
         return false;
     }
+
 
     return true;
@@ -1042 +1065 @@
     if (sfiles->outImage) {
         deleteFile(sfiles->outImage);
+    }
+
+    if (sfiles->outSources) {
+        deleteFile(sfiles->outSources);
     }
 

trunk/magic/remove/src/streaksremove.c

@@ -17 +17 @@
 static void writeImages(streakFiles *sf, bool exciseImageCube);
 static void updateAstrometry(streakFiles *sfiles);
+static void censorSources(streakFiles *sfiles, psU32 maskStreak);
 
 int
@@ -188 +189 @@
         }
 
+        censorSources(sfiles, maskStreak);
+
         // write the destreaked "temporary" images and the recovery images
         writeImages(sfiles, exciseImageCube);
+
 
         psLogMsg("streaksremove", PS_LOG_INFO, "time to process component %d: %f\n", sfiles->extnum, psTimerClear("PROCESS_COMPONENT"));
@@ -259 +263 @@
     fprintf(stderr, "\t-weight WEIGHT.fits: weight file to de-streak\n");
     fprintf(stderr, "\t-replace: replace the input images with the output\n");
-    fprintf(stderr, "\t-remove: remove the original image after processing (requires -replace)\n");
     fprintf(stderr, "\t-keepnonwarped: do not exise pixels that were not part of difference processing\n");
     fprintf(stderr, "\t-transparent val: instead of setting excicsed pixel to NAN add val\n");
@@ -335 +338 @@
     bool gotReplace = false;
     if ((argnum = psArgumentGet(argc, argv, "-replace"))) {
+        if (!nebulousImage) {
+            psError(PS_ERR_UNKNOWN, true, "replace is only supported for nebulous files");
+            usage();
+        }
         gotReplace = true;
         psArgumentRemove(argnum, &argc, argv);
@@ -416 +423 @@
         psMetadataAddStr(config->arguments, PS_LIST_TAIL, "INPUT.WEIGHT", 0,
                 "name of input weight image", argv[argnum]);
+        psArgumentRemove(argnum, &argc, argv);
+    }
+    if ((argnum = psArgumentGet(argc, argv, "-sources"))) {
+        psArgumentRemove(argnum, &argc, argv);
+        bool nebulousSources = IN_NEBULOUS(argv[argnum]);
+        if (nebulousSources != nebulousImage) {
+            psError(PS_ERR_UNKNOWN, true, "sources file must have %snebulous path with %s image path\n",
+                nebulousImage ? "" : "non ", nebulousImage ? "nebulous" : "regular");
+            usage();
+        }
+        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "SOURCES", 0,
+                "name of input sources file", argv[argnum]);
         psArgumentRemove(argnum, &argc, argv);
     }
@@ -757 +776 @@
 }
 
+// read a sources file (.cmf) and remove any rows whose coordinate is convered by a
+// streak mask
+static void
+censorSources(streakFiles *sfiles, psU32 maskStreak)
+{
+    if ((!sfiles->inSources) || (!sfiles->outMask)) {
+        return;
+    }
+    psImage *maskImage = sfiles->outMask->image;
+    if (!maskImage) {
+        psError(PS_ERR_IO, false, "maskImage is null!");
+        streaksExit("", PS_EXIT_PROG_ERROR);
+    }
+
+    sFile *in = sfiles->inSources;
+    sFile *out = sfiles->outSources;
+
+    in->header = psFitsReadHeader(NULL, in->fits);
+    if (!in->header) {
+        psError(PS_ERR_IO, false, "failed to read header from %s", in->resolved_name);
+        streaksExit("", PS_EXIT_DATA_ERROR);
+    }
+
+    bool status;
+    psString extname = psMetadataLookupStr(&status, in->header, "EXTNAME");
+    if (!extname) {
+        psError(PS_ERR_IO, false, "failed to find extname in header of %s", in->resolved_name);
+        streaksExit("", PS_EXIT_DATA_ERROR);
+    }
+
+    psArray *inTable = psFitsReadTable(in->fits);
+    if (!inTable->n) {
+        psError(PS_ERR_IO, false, "table in %s is empty", in->resolved_name);
+        streaksExit("", PS_EXIT_DATA_ERROR);
+    }
+
+    psArray *outTable = psArrayAllocEmpty(inTable->n);
+    int j = 0;
+    int numCensored = 0;
+    for (int i = 0 ; i < inTable->n; i++) {
+        psMetadata *row = inTable->data[i];
+
+        psF32 x = psMetadataLookupF32 (&status, row, "X_PSF");
+        psF32 y = psMetadataLookupF32 (&status, row, "Y_PSF");
+        
+        psU32 mask = psImageGet(maskImage, x, y);
+
+        // Key the source if the center pixel is not masked with maskStreak
+        if (! (mask & maskStreak) ) {
+            psArraySet(outTable, j++, row);
+        } else {
+            numCensored++;
+        }
+    }
+
+    printf("Censored %d sources\n", numCensored);
+
+    // get rid of unused elements (don't know if this is necessary)
+    psArrayRealloc(outTable, j);
+
+    addDestreakKeyword(in->header);
+    if (! psFitsWriteTable(out->fits, in->header, outTable, extname)) {
+        psError(PS_ERR_IO, false, "failed to write table to %s", out->resolved_name);
+        streaksExit("", PS_EXIT_DATA_ERROR);
+    }
+
+    if (!psFitsClose(out->fits)) {
+        psError(PS_ERR_IO, false, "failed to close table %s", out->resolved_name);
+        streaksExit("", PS_EXIT_DATA_ERROR);
+    }
+}

trunk/magic/remove/src/streaksremove.h

@@ -62 +62 @@
     sFile *outChMask;
     sFile *recChMask;
+    sFile *inSources;
+    sFile *outSources;
     psString class_id;
     pmFPAfile  *inAstrom;