IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26870


Ignore:
Timestamp:
Feb 10, 2010, 4:13:14 PM (16 years ago)
Author:
eugene
Message:

updates from trunk

Location:
branches/eam_branches/20091201/magic
Files:
1 deleted
14 edited

Legend:

Unmodified
Added
Removed
  • branches/eam_branches/20091201/magic/Makefile

    r26810 r26870  
    22default: all
    33
    4 all: ssa-core-cpp magic remove
     4all: ssa-core-cpp magic remove verify
    55
    6 install: magic.install remove.install
     6install: magic.install remove.install verify.install
    77
    88clean:
    99        if [ -d magic ]; then (cd magic && make -f Makefile.magic clean); fi
    1010        if [ -d ssa-core-cpp ]; then (cd ssa-core-cpp && make clean); fi
    11         cd $(REMOVE) && make -f Makefile.simple clean
    1211
    1312update:
  • branches/eam_branches/20091201/magic/censorObjects/src/censor.c

    r24695 r26870  
    1414
    1515static void usage(void) {
    16     fprintf(stderr, "USAGE: censor [-file image(s)] [-list imagelist] (output)\n");
     16    fprintf(stderr, "USAGE: censorObjects -file (sources) -masklist (masklist) (output)\n");
     17    fprintf(stderr, "       censorObjects -file (sources) -mask (mask.fits) (output)\n");
     18    fprintf(stderr, "       censorObjects -list (sourceslist) -masklist (masklist) (output)\n");
     19    fprintf(stderr, "       censorObjects -list (sourceslist) -mask (mask.fits) (output)\n");
    1720    exit(PS_EXIT_CONFIG_ERROR);
    1821}
  • branches/eam_branches/20091201/magic/remove

    • Property svn:ignore set to
      configure
      Makefile.in
      Doxyfile
      config.log
      depcomp
      config.status
      config.guess
      ltmain.sh
      config.sub
      autom4te.cache
      libtool
      missing
      Makefile
      aclocal.m4
      install-sh

  • branches/eam_branches/20091201/magic/remove/src

    • Property svn:ignore
      •  

        old new  
        44streakscompare
        55streaksrelease
        6 makefile
         6Makefile
         7Makefile.in
         8config.h
         9.deps
         10streaksVersionDefinitions.h
         11config.h.in
         12stamp-h1
         13
  • branches/eam_branches/20091201/magic/remove/src/Makefile.simple

    r26256 r26870  
    1212    streaksremove.o \
    1313    streaksextern.o \
    14     warpedpixels.o \
     14    diffedpixels.o \
    1515    Line.o
    1616
  • branches/eam_branches/20091201/magic/remove/src/isdestreaked.c

    r25082 r26870  
    11#include "streaksremove.h"
     2
     3char *streaksProgram = "isdestreaked";
    24
    35int
  • branches/eam_branches/20091201/magic/remove/src/streaksastrom.c

    r25082 r26870  
    379379{
    380380    if (!pmAstromLinearizeTransforms((pmFPA *) astrom->fpa, (pmChip *) astrom->chip)) {
    381         // streaksExit("linear fit to astrometry failed\n", PS_EXIT_UNKNOWN_ERROR);
    382         psError(PS_ERR_UNKNOWN, false, "linear fit to astrometry failed ignoring");
     381        psErrorStackPrint(stderr, "linear fit to astrometry failed. ignoring\n");
    383382        return false;
    384383    }
  • branches/eam_branches/20091201/magic/remove/src/streakscompare.c

    r25082 r26870  
    11#include "streaksremove.h"
     2
     3char *streaksProgram = "streakscompare";
    24
    35static pmConfig * parseArguments(int, char **);
  • branches/eam_branches/20091201/magic/remove/src/streaksio.c

    r26041 r26870  
    121121    sf->transparentStreaks = psMetadataLookupF64(&status, config->arguments, "TRANSPARENT_STREAKS");
    122122
     123    sf->stats = psMetadataAlloc();
     124    psString statsFileName= psMetadataLookupStr(&status, config->arguments, "STATS");
     125    if (statsFileName) {
     126        sf->statsFile = fopen(statsFileName, "w");
     127        if (!sf->statsFile) {
     128            psError(PS_ERR_IO, true, "failed to open stats file %s", statsFileName);
     129            streaksExit("", PS_EXIT_CONFIG_ERROR);
     130        }
     131    }
     132
    123133    return sf;
    124134}
     
    128138{
    129139    freeImages(sf);
    130     psFree(sf->warpedPixels);
     140    psFree(sf->diffedPixels);
    131141    psFree(sf->tiles);
    132142    psFree(sf->view);
     
    800810        return;
    801811    }
     812    streaksVersionHeaderFull(sfile->header);
    802813    if (!psFitsWriteImage(sfile->fits, sfile->header, sfile->image, 0, extname)) {
    803814        psError(PS_ERR_IO, false, "failed to write image to %s extnum: %d",
     
    817828        return;
    818829    }
     830    streaksVersionHeaderFull(sfile->header);
    819831    if (!psFitsWriteImageCube(sfile->fits, sfile->header, imagecube, extname)) {
    820832        psError(PS_ERR_IO, false, "failed to write image to %s extnum: %d",
     
    851863    sfile->header = NULL;
    852864    psFree(sfile->image);
    853     sfile->header = NULL;
     865    sfile->image = NULL;
    854866    psFree(sfile->imagecube);
    855     sfile->header = NULL;
     867    sfile->imagecube = NULL;
    856868}
    857869
     
    11811193}
    11821194
    1183 void
    1184 strkGetMaskValues(streakFiles *sfiles, psU32 *maskStreak, psU32 *maskMask)
    1185 {
     1195// Get the mask values that we need.
     1196// If an input mask file is provided get them from there.
     1197// Otherwise get them from the recipe MASKS
     1198void
     1199strkGetMaskValues(streakFiles *sfiles)
     1200{
     1201    if (sfiles->maskStreak) {
     1202        // already done
     1203        return;
     1204    }
    11861205    if (sfiles->inMask && sfiles->inMask->header) {
    11871206        if (!pmConfigMaskReadHeader(sfiles->config, sfiles->inMask->header)) {
     
    11961215        streaksExit("", PS_EXIT_CONFIG_ERROR);
    11971216    }
    1198     *maskStreak = psMetadataLookupU32(&status, masks, "STREAK");
     1217    sfiles->maskStreak = psMetadataLookupU32(&status, masks, "STREAK");
    11991218    if (!status) {
    12001219        psError(PM_ERR_CONFIG, false, "failed to lookup mask value for STREAK in recipes\n");
     
    12121231        }
    12131232    }
    1214     *maskMask = ~convPoor;
     1233    sfiles->maskMask = ~convPoor;
    12151234}
    12161235
  • branches/eam_branches/20091201/magic/remove/src/streaksio.h

    r24853 r26870  
    1616void copyFitsOptions(sFile *out, sFile *rec, sFile *in, psVector *tiles);
    1717void setupImageRefs(sFile *out, sFile *recoveryOut, sFile *in, int extnum, bool exciseAll);
    18 void strkGetMaskValues(streakFiles *sfiles, psU32 *maskStreak, psU32 *maskMask);
     18void strkGetMaskValues(streakFiles *sfiles);
    1919void copyExtraExtensions(streakFiles *sfiles);
    2020
  • branches/eam_branches/20091201/magic/remove/src/streaksrelease.c

    r25082 r26870  
    44static bool readAndCopyToOutput(streakFiles *sf, bool exciseAll);
    55static void writeImages(streakFiles *sf, bool exciseImageCube);
     6
     7char *streaksProgram = "streaksrelease";
    68
    79int
     
    1618        return PS_EXIT_CONFIG_ERROR;
    1719    }
    18 
    19     // Values to set for masked pixels
    20     psU32 maskStreak = 0;           // for the image and weight (usually NAN, MAXINT for integer images)
    21     psU32 maskMask = 0;             // value looked up for MASK.STREAK
    2220
    2321    // Does true work here?
     
    5048
    5149        // now that we've read the input files, lookup the mask values that we read
    52         if (maskStreak == 0) {
    53             strkGetMaskValues(sfiles, &maskStreak, &maskMask);
    54         }
    55 
    56         setMaskedToNAN(sfiles, maskMask, true);
     50        strkGetMaskValues(sfiles);
     51
     52        setMaskedToNAN(sfiles, sfiles->maskMask, true);
    5753
    5854        // write out the destreaked temporary images and the recovery images
  • branches/eam_branches/20091201/magic/remove/src/streaksremove.c

    r25434 r26870  
    1212static pmConfig *parseArguments(int argc, char **argv);
    1313static bool readAndCopyToOutput(streakFiles *sf, bool exciseAll);
    14 static void exciseNonWarpedPixels(streakFiles *sfiles, psImageMaskType newMaskValue);
    15 static bool warpedPixel(streakFiles *sfiles, int x, int y);
     14static void exciseNonDiffedPixels(streakFiles *sfiles, psImageMaskType newMaskValue);
     15static bool diffedPixel(streakFiles *sfiles, int x, int y);
    1616static void excisePixel(streakFiles *sfiles, unsigned int x, unsigned int y, bool streak, psImageMaskType newMaskValue);
    1717static void writeImages(streakFiles *sf, bool exciseImageCube);
    1818static void updateAstrometry(streakFiles *sfiles);
    1919static void censorSources(streakFiles *sfiles, psImageMaskType maskStreak);
    20 static long censorPixels(streakFiles *sfiles, psImage * pixels, bool checkNonWarpedPixels, psU16 maskStreak);
    21 
     20static long censorPixels(streakFiles *sfiles, psImage * pixels, bool checkNonDiffedPixels, psU16 maskStreak);
     21
     22char *streaksProgram = "streaksremove";
     23
     24// Note: For clarity the flow of this program is in main().
     25// There is not a lot of error checking is done in main.
     26// Until the end, where we might be doing Nebulous operations, called functions exit when an error
     27// is encountered.
    2228int
    2329main(int argc, char *argv[])
     
    2632
    2733    psLibInit(NULL);
    28     psTimerStart("STREAKSREMOVE");
     34    psTimerStart("TOTAL_TIME");
    2935
    3036    pmConfig *config = parseArguments(argc, argv);
     
    3339        streaksExit("", PS_EXIT_CONFIG_ERROR);
    3440    }
    35 
    36     // Values to set for masked pixels
    37     psU32 maskStreak = 0;           // for the image and weight (usually NAN, MAXINT for integer images)
    38     psU32 maskMask = 0;             // value looked up for MASK.STREAK
    3941
    4042    psString streaksFileName = psMetadataLookupStr(NULL, config->arguments, "STREAKS");
     
    5052    streakFiles *sfiles = openFiles(config, true, argv[0]);
    5153    setupAstrometry(sfiles);
     54    sfiles->stats = psMetadataAlloc();
    5255
    5356    // Optionally we can set pixels that are masked to NAN since they couldn't have been
     
    6164
    6265    bool exciseAll = false;
    63     // --keepnonwarped is a test and debug mode
    64     bool keepNonWarpedPixels = psMetadataLookupBool(&status, config->arguments, "KEEP_NON_WARPED");
    65 
    66     // we need to check for non warped pixels unless we've been asked not to or the stage is diff
    67     // (By definition pixels in diff images were included in a diff)
    68     bool checkNonWarpedPixels = ! (keepNonWarpedPixels || (sfiles->stage == IPP_STAGE_DIFF) );
    69 
    70     if (checkNonWarpedPixels ) {
     66    // --keepnondiffed is a test and debug mode
     67    bool keepNonDiffedPixels = psMetadataLookupBool(&status, config->arguments, "KEEP_NON_DIFFED");
     68
     69    // we need to check for non diffed pixels unless we've been asked not to or the stage is diff
     70    // (By definition pixels in diff images were included in the difference images)
     71    bool checkNonDiffedPixels = ! (keepNonDiffedPixels || (sfiles->stage == IPP_STAGE_DIFF) );
     72
     73    if (checkNonDiffedPixels ) {
    7174        // From magic ICD:
    7275        // In the raw and detrended images, the pixels which were not
     
    7679        // if no skycells are provided sfiles->exciseAll is set to true
    7780
    78         psTimerStart("COMPUTE_WARPED_PIXELS");
    79         if (! computeWarpedPixels(sfiles) ) {
     81        psTimerStart("COMPUTE_DIFFED_PIXELS");
     82        if (! computeDiffedPixels(sfiles) ) {
    8083            // we have no choice to excise all pixels
    8184            exciseAll = true;
    8285        }
    83         psF64 cwp_t = psTimerClear("COMPUTE_WARPED_PIXELS");
    84         psLogMsg("streaksremove", PS_LOG_INFO, "time to compute warped pixels: %f\n", cwp_t);
     86        psF64 cdp_t = psTimerClear("COMPUTE_DIFFED_PIXELS");
     87        psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "COMPUTE_UNDIFFED_PIXELS", PS_META_REPLACE, "time to compute non-diffedpixels", cdp_t);
     88        psLogMsg("streaksremove", PS_LOG_INFO, "time to compute diffed pixels: %f\n", cdp_t);
    8589    }
    8690
     
    101105    long totalStreakPixels = 0;
    102106
     107    // accumulators for the various timers
     108    psF64 gsp_t = 0;
     109    psF64 enw_t = 0;
     110    psF64 rms_t = 0;
     111    psF64 cs_t = 0;
     112    psF64 wi_t = 0;
     113    psF64 ua_t = 0;
    103114    // Iterate through each component of the input (except for raw images there is only one)
    104115    do {
     
    118129
    119130        // now that we've read the input files, lookup the mask values
    120         if (maskStreak == 0) {
    121             strkGetMaskValues(sfiles, &maskStreak, &maskMask);
    122         }
     131        strkGetMaskValues(sfiles);
    123132
    124133        totalPixels += sfiles->inImage->numRows * sfiles->inImage->numCols;
     
    132141            StreakPixels *pixels = streak_on_component(streaks, sfiles->astrom, sfiles->inImage->numCols,
    133142                                                        sfiles->inImage->numRows);
     143            gsp_t +=  psTimerClear("GET_STREAK_PIXELS");
     144            psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "GET_STREAK_PIXELS", PS_META_REPLACE, "", gsp_t);
    134145            psLogMsg("streaksremove", PS_LOG_INFO, "time to get streak pixels: %f\n", psTimerClear("GET_STREAK_PIXELS"));
    135146
     
    137148            // otherwise it contained an image cube (video cell) which is handled in the if block
    138149            if (sfiles->inImage->image) {
    139                 if (checkNonWarpedPixels) {
    140                     psTimerStart("EXCISE_NON_WARPED");
    141 
    142                     // set non-warped pixels and variance to NAN, mask to maskStreak (since the pixel
     150                if (checkNonDiffedPixels) {
     151                    psTimerStart("EXCISE_NON_DIFFED");
     152
     153                    // set non-diffed pixels and variance to NAN, mask to maskStreak (since the pixel
    143154                    // is excised as part of the destreaking process)
    144                     exciseNonWarpedPixels(sfiles, maskStreak);
    145 
    146                     psLogMsg("streaksremove", PS_LOG_INFO, "time to excise non warped pixels: %f\n", psTimerClear("EXCISE_NON_WARPED"));
     155                    exciseNonDiffedPixels(sfiles, sfiles->maskStreak);
     156
     157                    enw_t +=  psTimerClear("EXCISE_NON_DIFFED");
     158                    psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "EXCISE_NON_DIFFED", PS_META_REPLACE, "", enw_t);
     159                    psLogMsg("streaksremove", PS_LOG_INFO, "time to excise non diffed pixels: %f\n", enw_t);
    147160                }
    148161
    149162                psTimerStart("REMOVE_STREAKS");
    150163
    151                 totalStreakPixels += censorPixels(sfiles, pixels, checkNonWarpedPixels, maskStreak);
    152 
    153                 psLogMsg("streaksremove", PS_LOG_INFO, "time to remove streak pixels: %f\n", psTimerClear("REMOVE_STREAKS"));
     164                totalStreakPixels += censorPixels(sfiles, pixels, checkNonDiffedPixels, sfiles->maskStreak);
     165
     166                rms_t += psTimerClear("REMOVE_STREAKS");
     167                psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "REMOVE_STREAKS", PS_META_REPLACE, "", enw_t);
     168                psLogMsg("streaksremove", PS_LOG_INFO, "time to remove streak pixels: %f\n", rms_t);
    154169
    155170                if (nanForRelease) {
    156171                    // set any pixels that were masked, to NAN (unless they are already NAN)
    157                     setMaskedToNAN(sfiles, maskMask, true);
     172                    setMaskedToNAN(sfiles, sfiles->maskMask, true);
    158173                }
    159174
     
    170185            // chip processed files with the data calcuated by psastro at the camera stage
    171186            // (actually we use a linear approximation)
     187            psTimerStart("UPDATE_ASTROMETRY");
    172188            updateAstrometry(sfiles);
    173         }
    174 
    175         censorSources(sfiles, maskStreak);
     189            ua_t += psTimerClear("UPDATE_ASTROMETRY");
     190            psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "UPDATE_ASTROMETRY", PS_META_REPLACE, "", ua_t);
     191            psLogMsg("streaksremove", PS_LOG_INFO, "time to update astrometry: %f\n", ua_t);
     192        }
     193
     194        psTimerStart("CENSOR_SOURCES");
     195        censorSources(sfiles, sfiles->maskStreak);
     196        cs_t += psTimerClear("CENSOR_SOURCES");
     197        psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "CENSOR_SOURCES", PS_META_REPLACE, "", cs_t);
    176198
    177199        // write the destreaked "temporary" images and the recovery images
     200        psTimerStart("WRITE_IMAGES");
    178201        writeImages(sfiles, exciseImageCube);
    179 
     202        wi_t += psTimerClear("WRITE_IMAGES");
     203        psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "WRITE_IMAGES", PS_META_REPLACE, "", wi_t);
    180204
    181205        psLogMsg("streaksremove", PS_LOG_INFO, "time to process component %d: %f\n", sfiles->extnum, psTimerClear("PROCESS_COMPONENT"));
     206
    182207    } while (streakFilesNextExtension(sfiles));
    183208
     
    185210    psFree(streaks);
    186211
    187     psLogMsg("streaksremove", PS_LOG_INFO, "pixels: %ld streak pixels: %ld %4.2f%%\n", totalPixels, totalStreakPixels, 100. * totalStreakPixels / totalPixels);
     212    if (exciseAll) {
     213        totalStreakPixels = totalPixels;
     214    }
     215
     216    psF64 streakFraction = (double) totalStreakPixels / totalPixels;
     217    psLogMsg("streaksremove", PS_LOG_INFO, "pixels: %ld streak pixels: %ld %4.2f%%\n", totalPixels, totalStreakPixels, streakFraction * 100);
     218    psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "STREAK_FRACTION", PS_META_REPLACE, "", streakFraction);
    188219
    189220    // check the weight and mask files for extra extensions that might be in files
     
    194225
    195226    psTimerStart("CLOSE_IMAGES");
    196 
    197227    closeImages(sfiles);
    198 
    199     psLogMsg("streaksremove", PS_LOG_INFO, "time to close images: %f\n", psTimerClear("CLOSE_IMAGES"));
    200 
    201 
     228    psF64 ci_t = psTimerClear("CLOSE_IMAGES");
     229    psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "CLOSE_IMAGES", PS_META_REPLACE, "", ci_t);
     230
     231    psLogMsg("streaksremove", PS_LOG_INFO, "time to close images: %f\n", ci_t);
     232
     233    psTimerStart("REPLICATE_OUTPUTS");
    202234    if (!replicateOutputs(sfiles)) {
    203235        psErrorStackPrint(stderr, "failed to replicate output files");
     
    205237        exit(PS_EXIT_UNKNOWN_ERROR);
    206238    }
     239    psF64 ro_t = psTimerClear("REPLICATE_OUTPUTS");
     240    psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "REPLICATE_OUTPUTS", PS_META_REPLACE, "", ro_t);
    207241
    208242    // NOTE: from here on we can't just quit if something goes wrong.
     
    214248        //     swap the instances for the input and output
    215249        //     Note this is a nebulous database operation. No file I/O is performed
     250        psTimerStart("SWAP_INSTANCES");
    216251        if (!swapOutputsToInputs(sfiles)) {
    217             // XXX: Now what?
    218252            // It is up to the program that reverts failed destreak runs to insure that
    219             // any input files that have been swapped are restored and that the de-streaked
    220             // versions are deleted
     253            // any original non-destreaked input files that have been swapped are restored and that the de-streaked
     254            // versions are deleted.
    221255
    222256            psErrorStackPrint(stderr, "failed to swap files");
     
    224258            // XXX: pick a specific error code for this failure
    225259            exit(PS_EXIT_UNKNOWN_ERROR);
     260        }
     261        psF64 si_t = psTimerClear("SWAP_INSTANCES");
     262        psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "SWAP_INSTANCES", PS_META_REPLACE, "", si_t);
     263    }
     264
     265    psF64 total_time = psTimerClear("TOTAL_TIME");
     266    psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "TOTAL_TIME", PS_META_REPLACE, "", total_time);
     267    psLogMsg("streaksremove", PS_LOG_INFO, "time to run streaksremove: %f\n", total_time);
     268
     269    if (sfiles->statsFile) {
     270        const char *statsMDC = psMetadataConfigFormat(sfiles->stats);
     271        if (!statsMDC || strlen(statsMDC) == 0) {
     272            psError(PS_ERR_IO, false, "Unable to get statistics MDC file.\n");
     273        } else {
     274            fprintf(sfiles->statsFile, "%s", statsMDC);
     275            psFree(statsMDC);
     276            fclose(sfiles->statsFile);
     277            sfiles->statsFile = NULL;
     278            psFree(sfiles->stats);
     279            sfiles->stats = NULL;
    226280        }
    227281    }
     
    234288    streaksNebulousCleanup();
    235289    pmConfigDone();
    236     psLogMsg("streaksremove", PS_LOG_INFO, "time to run streaksremove: %f\n", psTimerClear("STREAKSREMOVE"));
    237290    psLibFinalize();
    238291
     
    243296
    244297static long
    245 censorPixels(streakFiles *sfiles, psImage *pixels, bool checkNonWarpedPixels, psU16 maskStreak)
     298censorPixels(streakFiles *sfiles, psImage *pixels, bool checkNonDiffedPixels, psU16 maskStreak)
    246299{
    247300    long streakPixels = 0;
     
    251304            if (psImageGet(pixels, x, y)) {
    252305                ++streakPixels;
    253                 if (!checkNonWarpedPixels || warpedPixel(sfiles, x, y)) {
     306                if (!checkNonDiffedPixels || diffedPixel(sfiles, x, y)) {
    254307
    255308                    excisePixel(sfiles, x, y, true, maskStreak);
     
    257310                } else {
    258311                    // This pixel was not included in any warp and has thus already excised
    259                     // by exciseNonWarpedPixels
     312                    // by exciseNonDiffedPixels
    260313                }
    261314            }
     
    274327    fprintf(stderr, "\t-weight WEIGHT.fits: weight file to de-streak\n");
    275328    fprintf(stderr, "\t-replace: replace the input images with the output\n");
    276     fprintf(stderr, "\t-keepnonwarped: do not exise pixels that were not part of difference processing\n");
     329    fprintf(stderr, "\t-keepnondiffed: do not exise pixels that were not part of difference processing\n");
    277330    fprintf(stderr, "\t-transparent val: instead of setting excicsed pixel to NAN add val\n");
    278331    fprintf(stderr, "\t-chip_mask MASK.fits: name of mask for chip stage (camera stage mask is passed tih -mask)\n");
     
    365418    }
    366419
    367     if ((argnum = psArgumentGet(argc, argv, "-keepnonwarped"))) {
    368         psArgumentRemove(argnum, &argc, argv);
    369         psMetadataAddBool(config->arguments, PS_LIST_TAIL, "KEEP_NON_WARPED", 0,
    370             "skip excising of non warped pixels", true);
     420    if ((argnum = psArgumentGet(argc, argv, "-keepnondiffed"))) {
     421        psArgumentRemove(argnum, &argc, argv);
     422        psMetadataAddBool(config->arguments, PS_LIST_TAIL, "KEEP_NON_DIFFED", 0,
     423            "skip excising of non diffed pixels", true);
    371424    }
    372425
     
    378431        psArgumentRemove(argnum, &argc, argv);
    379432        double transparentStreaks = atof(argv[argnum]);
    380         psMetadataAddF64(config->arguments, PS_LIST_TAIL, "TRANSPARENT_STREAKS", 0,
     433        psMetadataAddF32(config->arguments, PS_LIST_TAIL, "TRANSPARENT_STREAKS", 0,
    381434            "value to adjust excised pixels", transparentStreaks);
    382435        psArgumentRemove(argnum, &argc, argv);
    383436    }
    384437
    385     // if skycells are not provided then we have to execise all pixels  unless -keepnonwarped
     438    // if skycells are not provided then we have to execise all pixels  unless -keepnondiffed
    386439    pmConfigFileSetsMD(config->arguments, &argc, argv, "SKYCELLS", "-skycell", "-skycelllist");
    387440
     
    449502    }
    450503
     504    if ((argnum = psArgumentGet(argc, argv, "-stats"))) {
     505        psArgumentRemove(argnum, &argc, argv);
     506        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "STATS", 0,
     507                "name of input stats file", argv[argnum]);
     508        psArgumentRemove(argnum, &argc, argv);
     509    }
     510
    451511    if ((argnum = psArgumentGet(argc, argv, "-tmproot"))) {
    452512        psArgumentRemove(argnum, &argc, argv);
     
    491551updateAstrometry(streakFiles *sf)
    492552{
    493     // XXX: why do I check this here? Shouldn't it be just around the call to linearizeTransforms?
    494553    if (sf->bilevelAstrometry) {
    495 
    496554        if (!linearizeTransforms(sf->astrom)) {
    497             // fit failed, leave the astrometry unchanged
     555            // fit failed, leave the transform in the file unchanged
    498556            return;
    499557        }
    500 
    501         if (!pmAstromWriteWCS(sf->outImage->header, sf->inAstrom->fpa, sf->chip, 0.001)) {
    502             psError(PS_ERR_UNKNOWN, false, "failed to update astrometry for extension %d", sf->extnum);
    503             streaksExit("", PS_EXIT_UNKNOWN_ERROR);
    504         }
    505         if (sf->outMask) {
    506             pmAstromWriteWCS(sf->outMask->header, sf->inAstrom->fpa, sf->chip, 0.001);
    507         }
    508         if (sf->outWeight) {
    509             pmAstromWriteWCS(sf->outWeight->header, sf->inAstrom->fpa, sf->chip, 0.001);
     558    }
     559    if (!pmAstromWriteWCS(sf->outImage->header, sf->inAstrom->fpa, sf->chip, 0.001)) {
     560        psError(PS_ERR_UNKNOWN, false, "failed to update astrometry for extension %d", sf->extnum);
     561        streaksExit("", PS_EXIT_UNKNOWN_ERROR);
     562    }
     563    if (sf->outMask) {
     564        pmAstromWriteWCS(sf->outMask->header, sf->inAstrom->fpa, sf->chip, 0.001);
     565    }
     566    if (sf->outWeight) {
     567        pmAstromWriteWCS(sf->outWeight->header, sf->inAstrom->fpa, sf->chip, 0.001);
     568    }
     569}
     570
     571static void
     572setStreakBits(psImage *maskImage, psU32 maskStreak)
     573{
     574    for (int y=0 ; y < maskImage->numRows; y++) {
     575        for (int x=0 ; x < maskImage->numCols; x++) {
     576            maskImage->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskStreak;
    510577        }
    511578    }
     
    516583readAndCopyToOutput(streakFiles *sf, bool exciseAll)
    517584{
    518     bool    updateAstrometry = false;
    519585    if (sf->inImage->pmfile) {
    520586        // image data from pmFPAfile (diff or warp file)
     
    534600        streaksExit("", PS_EXIT_UNKNOWN_ERROR);
    535601    }
    536     // For the chip level files, copy the WCS from the astrometry file to the header
    537     // XXX: do we want to do this for raw images as well?
    538     if (sf->stage == IPP_STAGE_CHIP) {
    539         if (!sf->bilevelAstrometry) {
    540             updateAstrometry = true;
    541             if (!pmAstromWriteWCS(sf->inImage->header, sf->inAstrom->fpa, sf->chip, 0.001)) {
    542                 psError(PS_ERR_UNKNOWN, false, "failed to update astrometry for extension %d", sf->extnum);
    543                 streaksExit("", PS_EXIT_UNKNOWN_ERROR);
    544             }
    545         }
    546     }
    547602    sf->outImage->header =  psMemIncrRefCounter(sf->inImage->header);
    548603    if (sf->recImage) {
     
    577632            }
    578633            addDestreakKeyword(sf->outMask->header);
    579             if (updateAstrometry) {
    580                 pmAstromWriteWCS(sf->outMask->header, sf->inAstrom->fpa, sf->chip, 0.001);
    581             }
    582             setupImageRefs(sf->outMask, sf->recMask, sf->inMask, sf->extnum, exciseAll);
     634            // Note: we don't excise the mask pixels even if exciseAll is true.
     635            setupImageRefs(sf->outMask, sf->recMask, sf->inMask, sf->extnum, false);
     636            if (exciseAll) {
     637                strkGetMaskValues(sf);
     638               
     639                // add the STREAK bit to the mask image pixels
     640                setStreakBits(sf->inMask->image, sf->maskStreak);
     641            }
    583642            if (sf->outChMask) {
    584643                sf->outChMask->header = (psMetadata *) psMemIncrRefCounter(sf->outMask->header);
     
    605664        }
    606665        addDestreakKeyword(sf->outWeight->header);
    607         if (updateAstrometry) {
    608             pmAstromWriteWCS(sf->inWeight->header, sf->inAstrom->fpa, sf->chip, 0.001);
    609         }
    610666        setupImageRefs(sf->outWeight, sf->recWeight, sf->inWeight, sf->extnum, exciseAll);
    611667
     
    630686    } else {
    631687        // we have an image cube
    632         double initValue;
    633688        if (exciseImageCube) {
    634689            // copy the entire input image to the recovery image
    635690            writeImageCube(sf->recImage, sf->inImage->imagecube, extname, sf->extnum);
    636             initValue = NAN;
    637691        } else {
    638692            // otherwise write it to the output
    639693            writeImageCube(sf->outImage, sf->inImage->imagecube, extname, sf->extnum);
    640             initValue = 0;
    641         }
    642 
    643         // borrow one of the images from the imagecube and set it to init value
    644         psImage *image = psArrayGet (sf->inImage->imagecube, 0);
    645         psMemIncrRefCounter(image);
    646         psImageInit(image, initValue);
     694        }
     695
     696        // Now deal with the other output image
    647697        if (exciseImageCube) {
    648             sf->outImage->image = image;
    649             writeImage(sf->outImage, extname, sf->extnum);
     698            // Set the values in the imagecube images to NAN and write them to the output image
     699            for (int i = 0; i < psArrayLength(sf->inImage->imagecube); i++) {
     700                psImage *image = psArrayGet (sf->inImage->imagecube, i);
     701                // XXX: NAN isn't right. It should be the integer equivalent. Should use exciseValue
     702                // but it isn't set with this code path. Fix that.
     703                psImageInit(image, 65535);
     704            }
     705            writeImageCube(sf->outImage, sf->inImage->imagecube, extname, sf->extnum);
    650706        } else {
    651             // write zero valued image to reccovery
    652707            if (sf->recImage) {
    653                 sf->recImage->image = image;
    654                 writeImage(sf->recImage, extname, sf->extnum);
     708                // Set the values in the imagecube images to zero
     709                for (int i = 0; i < psArrayLength(sf->inImage->imagecube); i++) {
     710                    psImage *image = psArrayGet (sf->inImage->imagecube, i);
     711                    psImageInit(image, 0);
     712                }
     713                // copy the entire zeroed image to the recovery image
     714                writeImageCube(sf->recImage, sf->inImage->imagecube, extname, sf->extnum);
    655715            }
    656716        }
     
    683743    }
    684744
    685     float imageValue  = sfiles->inImage->image->data.F32[y][x];
    686     if (sfiles->recImage && !isExciseValue(imageValue, sfiles->inImage->exciseValue) ) {
    687         sfiles->recImage->image->data.F32[y][x] = imageValue;
    688     }
    689 
    690     if (sfiles->transparentStreaks == 0) {
    691         sfiles->outImage->image->data.F32[y][x] = exciseValue;
     745    if (sfiles->inImage->image->type.type == PS_TYPE_U16) {
     746        psU16 imageValue  = sfiles->inImage->image->data.U16[y][x];
     747        if (sfiles->recImage && !isExciseValue(imageValue, sfiles->inImage->exciseValue) ) {
     748            sfiles->recImage->image->data.U16[y][x] = imageValue;
     749        }
     750
     751        if (sfiles->transparentStreaks == 0) {
     752            sfiles->outImage->image->data.U16[y][x] = exciseValue;
     753        } else {
     754            if (streak) {
     755                // as a visualization aid don't mask the pixel, just change the intensity
     756                sfiles->outImage->image->data.U16[y][x] = imageValue + sfiles->transparentStreaks;
     757            } else {
     758                sfiles->outImage->image->data.U16[y][x] = exciseValue;
     759            }
     760        }
    692761    } else {
    693         if (streak) {
    694             // as a visualization aid don't mask the pixel, just change the intensity
    695             sfiles->outImage->image->data.F32[y][x] = imageValue + sfiles->transparentStreaks;
     762        float imageValue  = sfiles->inImage->image->data.F32[y][x];
     763        if (sfiles->recImage && !isExciseValue(imageValue, sfiles->inImage->exciseValue) ) {
     764            sfiles->recImage->image->data.F32[y][x] = imageValue;
     765        }
     766
     767        if (sfiles->transparentStreaks == 0) {
     768            sfiles->outImage->image->data.F32[y][x] = exciseValue;
    696769        } else {
    697             sfiles->outImage->image->data.F32[y][x] = exciseValue;
     770            if (streak) {
     771                // as a visualization aid don't mask the pixel, just change the intensity
     772                sfiles->outImage->image->data.F32[y][x] = imageValue + sfiles->transparentStreaks;
     773            } else {
     774                sfiles->outImage->image->data.F32[y][x] = exciseValue;
     775            }
    698776        }
    699777    }
     
    711789                sfiles->inMask->image->data.PS_TYPE_IMAGE_MASK_DATA[y][x];
    712790        }
    713         sfiles->outMask->image->data.PS_TYPE_IMAGE_MASK_DATA[y][x] =
    714             sfiles->inMask->image->data.PS_TYPE_IMAGE_MASK_DATA[y][x] | newMaskValue;
     791        sfiles->outMask->image->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= newMaskValue;
    715792    }
    716793}
    717794
    718795static void
    719 exciseNonWarpedPixels(streakFiles *sfiles, psImageMaskType newMaskValue)
     796exciseNonDiffedPixels(streakFiles *sfiles, psImageMaskType newMaskValue)
    720797{
    721798    int cell_x0 = sfiles->astrom->cell_x0;
     
    739816        }
    740817
    741         psU8 *pixels = sfiles->warpedPixels->data.U8[yChip];
     818        psU8 *pixels = sfiles->diffedPixels->data.U8[yChip];
    742819
    743820        if (xParity == 1) {
     
    762839
    763840static bool
    764 warpedPixel(streakFiles *sfiles, int x, int y)
     841diffedPixel(streakFiles *sfiles, int x, int y)
    765842{
    766843    PixelPos chipCoord;
    767844
    768845    if (!CHIP_LEVEL_INPUT(sfiles->stage)) {
    769         // if we're here on a skycell image by definition this pixel was warped
     846        // if we're here on a skycell image by definition this pixel was diffed
    770847        return true;
    771848    }
     
    780857    cellToChipInt(&chipCoord.x, &chipCoord.y, sfiles->astrom, x, y);
    781858
    782     if (chipCoord.x < 0 || chipCoord.x >= sfiles->warpedPixels->numCols) {
     859    if (chipCoord.x < 0 || chipCoord.x >= sfiles->diffedPixels->numCols) {
    783860        return false;
    784861    }
    785     if (chipCoord.y < 0 || chipCoord.y >= sfiles->warpedPixels->numRows) {
     862    if (chipCoord.y < 0 || chipCoord.y >= sfiles->diffedPixels->numRows) {
    786863        return false;
    787864    }
    788865
    789     return psImageGet(sfiles->warpedPixels, chipCoord.x, chipCoord.y) ? true : false;
     866    return psImageGet(sfiles->diffedPixels, chipCoord.x, chipCoord.y) ? true : false;
    790867}
    791868
  • branches/eam_branches/20091201/magic/remove/src/streaksremove.h

    r25082 r26870  
    3333    int         yParity;
    3434    double      exciseValue;
     35    psString    program;
    3536} sFile;
    3637
     38// used for error messages
     39extern char * streaksProgram;
    3740
    3841typedef enum {
     
    6467    sFile *inSources;
    6568    sFile *outSources;
     69    psMetadata *stats;
     70    FILE *statsFile;
    6671    psString class_id;
    6772    pmFPAfile  *inAstrom;
     
    7176    pmChip  *chip;  // current chip
    7277    pmCell  *cell;  // current cell
    73     psImage *warpedPixels;
     78    psImage *diffedPixels;
    7479    psVector    *tiles;
    7580    double  transparentStreaks;
    7681    psString    program_name;
     82    psU32   maskStreak;
     83    psU32   maskMask;
    7784} streakFiles;
    7885
     
    8693extern void cellToChipInt(unsigned int *xChip, unsigned int *yChip, strkAstrom *astrom, int xCell, int yCell);
    8794
    88 extern bool computeWarpedPixels(streakFiles *sf);
     95extern bool computeDiffedPixels(streakFiles *sf);
    8996extern void streaksExit(psString, int);
    9097extern ippStage parseStage(psString);
    9198extern psString pathToDirectory(char *path);
    9299extern void setStreakFiles( streakFiles *);
     100
     101extern bool streaksVersionHeaderFull(psMetadata *header);
     102extern psString streaksVersionLong(void);
    93103
    94104#define CHIP_LEVEL_INPUT(_stage) ((_stage == IPP_STAGE_RAW) || (_stage == IPP_STAGE_CHIP))
  • branches/eam_branches/20091201/magic/remove/src/streaksreplace.c

    r25082 r26870  
    55static void writeImages(streakFiles *sf, bool exciseImageCube);
    66static bool readAndCopyToOutput(streakFiles *sf, bool restoreImageCube);
     7
     8char *streaksProgram = "streaksreplace";
    79
    810int main(int argc, char *argv[])
Note: See TracChangeset for help on using the changeset viewer.