IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 41413


Ignore:
Timestamp:
Sep 10, 2020, 9:40:15 AM (6 years ago)
Author:
eugene
Message:

gDiff now runs and does not crash on forced photometry

Location:
branches/eam_branches/gDiff.20200202/src
Files:
16 added
12 edited

Legend:

Unmodified
Added
Removed
  • branches/eam_branches/gDiff.20200202/src/Makefile.am

    r41260 r41413  
    88gDiffVersion.c: gDiffVersionDefinitions.h
    99gDiffVersionDefinitions.h: gDiffVersionDefinitions.h.in FORCE
    10         ../pslib-setsvnversion.pl GDIFF gDiffVersionDefinitions.h.in gDiffVersionDefinitions.h
     10        pslib-setsvnversion.pl GDIFF gDiffVersionDefinitions.h.in gDiffVersionDefinitions.h
    1111FORCE: ;
    1212
     
    1717        gDiff.c                         \
    1818        gDiffArguments.c                \
    19         gDiffVersion.c                  \
     19        gDiffBackground.c               \
    2020        gDiffCamera.c                   \
    2121        gDiffData.c                     \
     22        gDiffDefineOutput.c             \
    2223        gDiffErrorCodes.c               \
     24        gDiffExit.c                     \
     25        gDiffExtras.c                   \
     26        gDiffFiles.c                    \
     27        gDiffFlagNeighbors.c            \
     28        gDiffInputDetections.c          \
    2329        gDiffLoop.c                     \
    24         gDiffExit.c                     \
    25         gDiffFiles.c                    \
    26         gDiffDefineOutput.c             
    27 
    28 #       gDiffBackground.c               \
    29 #       gDiffVarianceRescale.c          \
    30 #       gDiffInputDetections.c          \
    31 #       gDiffExtras.c                   \
    32 #       gDiffFlagNeighbors.c            \
    33 #       gDiffMakePSF.c                  \
    34 #       gDiffMatchPSFs.c                \
    35 #       gDiffMaskStats.c                \
    36 #       gDiffReadoutInverse.c           \
    37 #       gDiffReadoutJpeg.c              \
    38 #       gDiffReadoutPhotometry.c        \
    39 #       gDiffReadoutForcedPhot.c        \
    40 #       gDiffReadoutStats.c             \
    41 #       gDiffReadoutSubtract.c          \
    42 #       gDiffSetMasks.c                 \
    43 #       gDiffThreshold.c
     30        gDiffMakePSF.c                  \
     31        gDiffMaskStats.c                \
     32        gDiffMatchPSFs.c                \
     33        gDiffReadoutForcedPhot.c        \
     34        gDiffReadoutInverse.c           \
     35        gDiffReadoutJpeg.c              \
     36        gDiffReadoutPhotometry.c        \
     37        gDiffReadoutStats.c             \
     38        gDiffReadoutSubtract.c          \
     39        gDiffSetMasks.c                 \
     40        gDiffThreshold.c                \
     41        gDiffVarianceRescale.c          \
     42        gDiffVersion.c                 
    4443
    4544noinst_HEADERS = \
  • branches/eam_branches/gDiff.20200202/src/gDiff.c

    r41260 r41413  
    1818    psphotErrorRegister();
    1919
     20    psExit exitValue = PS_EXIT_SUCCESS; // default value
     21
    2022    gDiffData *data = NULL; // Processing data
    2123    pmConfig *config = pmConfigRead(&argc, argv, GDIFF_RECIPE); // Configuration
     
    2628    psphotInit(); // cannot return false
    2729
     30    // the 'data' container carries pointers to all data structures used by gDiff
    2831    data = gDiffDataAlloc(config);
    2932
     33    // parse the command-line options
    3034    if (!gDiffArguments(argc, argv, data)) {
    3135        psError(psErrorCodeLast(), false, "Error reading arguments.");
     
    3337    }
    3438
    35     // determine the camera and setup the I/O file handles
     39    // determine the camera and setup the (many) I/O file handles
    3640    if (!gDiffCamera(data)) goto die;
    3741
     
    4044
    4145 die:
    42     {
    43         psExit exitValue = gDiffExitCode(PS_EXIT_SUCCESS); // Exit code
     46    exitValue = gDiffExitCode(exitValue); // Exit code
    4447
    45         if (data && data->stats && data->statsFile) {
    46             psMetadataAddF32(data->stats, PS_LIST_TAIL, "TIME_DIFF", 0, "Total time (sec)", psTimerMark("gDiff"));
    47             psString stats = psMetadataConfigFormat(data->stats); // Statistics to output
    48             if (!stats || strlen(stats) == 0) {
    49                 psError(GDIFF_ERR_IO, false, "Unable to format statistics file");
    50             } else if (fprintf(data->statsFile, "%s", stats) != strlen(stats)) {
    51                 psError(GDIFF_ERR_IO, true, "Unable to write statistics file");
    52             }
    53             psFree(stats);
    54             if (fclose(data->statsFile) == EOF) {
    55                 psError(GDIFF_ERR_IO, true, "Unable to close statistics file");
    56             }
    57             data->statsFile = NULL;
    58             pmConfigRunFilenameAddWrite(data->config, "STATS", data->statsName);
    59             exitValue = gDiffExitCode(exitValue);
    60         }
     48    if (data && data->stats && data->statsFile) {
     49      psMetadataAddF32(data->stats, PS_LIST_TAIL, "TIME_DIFF", 0, "Total time (sec)", psTimerMark("gDiff"));
     50      psString stats = psMetadataConfigFormat(data->stats); // Statistics to output
     51      if (!stats || strlen(stats) == 0) {
     52        psError(GDIFF_ERR_IO, false, "Unable to format statistics file");
     53      } else if (fprintf(data->statsFile, "%s", stats) != strlen(stats)) {
     54        psError(GDIFF_ERR_IO, true, "Unable to write statistics file");
     55      }
     56      psFree(stats);
     57      if (fclose(data->statsFile) == EOF) {
     58        psError(GDIFF_ERR_IO, true, "Unable to close statistics file");
     59      }
     60      data->statsFile = NULL;
     61      pmConfigRunFilenameAddWrite(data->config, "STATS", data->statsName);
     62      exitValue = gDiffExitCode(exitValue);
     63    }
    6164
    62         if (config && !gDiffFilesIterateUp(config, GDIFF_FILES_ALL)) {
    63             psError(psErrorCodeLast(), false, "Unable to close files.");
    64             exitValue = gDiffExitCode(exitValue);
    65             pmFPAfileFreeSetStrict(false);
    66         }
     65    if (config && !gDiffFilesIterateUp(config, GDIFF_FILES_ALL)) {
     66      psError(psErrorCodeLast(), false, "Unable to close files.");
     67      exitValue = gDiffExitCode(exitValue);
     68      pmFPAfileFreeSetStrict(false);
     69    }
    6770
    68         if (data) {
    69             psString dump_file = psMetadataLookupStr(NULL, data->config->arguments, "DUMP_CONFIG");
    70             if (dump_file) {
    71                 if (!pmConfigDump(data->config, dump_file)) {
    72                     psError(GDIFF_ERR_IO, false, "Unable to dump configuration.");
    73                     exitValue = gDiffExitCode(exitValue);
    74                 }
    75             }
    76             psFree(data);
    77         }
     71    if (data) {
     72      psString dump_file = psMetadataLookupStr(NULL, data->config->arguments, "DUMP_CONFIG");
     73      if (dump_file) {
     74        if (!pmConfigDump(data->config, dump_file)) {
     75          psError(GDIFF_ERR_IO, false, "Unable to dump configuration.");
     76          exitValue = gDiffExitCode(exitValue);
     77        }
     78      }
     79      psFree(data);
     80    }
    7881
    79         psTrace("gDiff", 1, "Finished at %f sec\n", psTimerMark("gDiff"));
    80         psLogMsg("gDiff", PS_LOG_INFO, "Complete gDiff run: %f sec\n", psTimerMark("gDiff"));
    81         psTimerStop();
     82    psTrace("gDiff", 1, "Finished at %f sec\n", psTimerMark("gDiff"));
     83    psLogMsg("gDiff", PS_LOG_INFO, "Complete gDiff run: %f sec\n", psTimerMark("gDiff"));
     84    psTimerStop();
    8285
    83         pmVisualClose(); // close plot windows, if -visual is set
    84         pmModelClassCleanup();
    85         pmConfigDone();
    86         pmVisualCleanup ();
    87         psLibFinalize();
     86    pmVisualClose(); // close plot windows, if -visual is set
     87    pmModelClassCleanup();
     88    pmConfigDone();
     89    pmVisualCleanup ();
     90    psLibFinalize();
    8891
    89         fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "gDiff");
    90         exitValue = gDiffExitCode(exitValue);
    91         exit(exitValue);
    92     }
     92    fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "gDiff");
     93    exitValue = gDiffExitCode(exitValue);
     94    exit(exitValue);
    9395}
  • branches/eam_branches/gDiff.20200202/src/gDiff.h

    r41260 r41413  
    2121#include <pslib.h>
    2222#include <psmodules.h>
     23#include <ppStats.h>
    2324#include <psphot.h>
    2425
  • branches/eam_branches/gDiff.20200202/src/gDiffArguments.c

    r41260 r41413  
     1/** @file gDiffArguments.c
     2 *  @brief: image subtraction based on matched Gaussians profiles
     3 *  @author Eugene Magnier @ IfA
     4 *  @version $Revision: 1.12 $
     5 *  @date $Date: 2020-01-25 $
     6 *  Copyright 2020 Institute for Astronomy, University of Hawaii
     7 */
     8
    19#include "gDiff.h"
    210
     
    6068
    6169    psMetadata *arguments = config->arguments; // Command-line arguments
    62     psMetadataAddStr(arguments, PS_LIST_TAIL, "-inimage", 0, "Input image", NULL);
    63     psMetadataAddStr(arguments, PS_LIST_TAIL, "-inmask", 0, "Input mask image", NULL);
    64     psMetadataAddStr(arguments, PS_LIST_TAIL, "-invariance", 0, "Input variance image", NULL);
    65     psMetadataAddStr(arguments, PS_LIST_TAIL, "-insources", 0, "Input source list", NULL);
    66     psMetadataAddStr(arguments, PS_LIST_TAIL, "-refimage", 0, "Reference image", NULL);
    67     psMetadataAddStr(arguments, PS_LIST_TAIL, "-refmask", 0, "Reference mask image", NULL);
    68     psMetadataAddStr(arguments, PS_LIST_TAIL, "-refvariance", 0, "Reference variance image", NULL);
    69     psMetadataAddStr(arguments, PS_LIST_TAIL, "-refsources", 0, "Reference source list", NULL);
    70     psMetadataAddStr(arguments, PS_LIST_TAIL, "-kernel", 0, "Precalculated kernel to apply", NULL);
    71     psMetadataAddStr(arguments, PS_LIST_TAIL, "-stats", 0, "Statistics file", NULL);
    72     psMetadataAddStr(arguments, PS_LIST_TAIL, "-stamps", 0, "Stamps filename; x,y on each line", NULL);
    73     psMetadataAddS32(arguments, PS_LIST_TAIL, "-convolve", 0, "Image to convolve [1 or 2]", 0);
     70    psMetadataAddStr (arguments, PS_LIST_TAIL, "-inimage", 0, "Input image", NULL);
     71    psMetadataAddStr (arguments, PS_LIST_TAIL, "-inmask", 0, "Input mask image", NULL);
     72    psMetadataAddStr (arguments, PS_LIST_TAIL, "-invariance", 0, "Input variance image", NULL);
     73    psMetadataAddStr (arguments, PS_LIST_TAIL, "-insources", 0, "Input source list", NULL);
     74    psMetadataAddStr (arguments, PS_LIST_TAIL, "-refimage", 0, "Reference image", NULL);
     75    psMetadataAddStr (arguments, PS_LIST_TAIL, "-refmask", 0, "Reference mask image", NULL);
     76    psMetadataAddStr (arguments, PS_LIST_TAIL, "-refvariance", 0, "Reference variance image", NULL);
     77    psMetadataAddStr (arguments, PS_LIST_TAIL, "-refsources", 0, "Reference source list", NULL);
     78    psMetadataAddStr (arguments, PS_LIST_TAIL, "-kernel", 0, "Precalculated kernel to apply", NULL);
     79    psMetadataAddStr (arguments, PS_LIST_TAIL, "-stats", 0, "Statistics file", NULL);
     80    psMetadataAddStr (arguments, PS_LIST_TAIL, "-stamps", 0, "Stamps filename; x,y on each line", NULL);
     81    psMetadataAddS32 (arguments, PS_LIST_TAIL, "-convolve", 0, "Image to convolve [1 or 2]", 0);
    7482    psMetadataAddBool(arguments, PS_LIST_TAIL, "-photometry", 0, "Perform photometry?", NULL);
    7583    psMetadataAddBool(arguments, PS_LIST_TAIL, "-forced-phot", 0, "Perform forced photometry?", NULL);
    7684    psMetadataAddBool(arguments, PS_LIST_TAIL, "-forced-input1", 0, "Perform forced photometry?", NULL);
    7785    psMetadataAddBool(arguments, PS_LIST_TAIL, "-forced-input2", 0, "Perform forced photometry?", NULL);
    78     psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp", 0, "Zero point for photometry", NAN);
     86    psMetadataAddF32 (arguments, PS_LIST_TAIL, "-zp", 0, "Zero point for photometry", NAN);
    7987    psMetadataAddBool(arguments, PS_LIST_TAIL, "-inverse", 0, "Generate inverse subtractions?", false);
    8088    psMetadataAddBool(arguments, PS_LIST_TAIL, "-save-inconv", 0, "Save input convolved images?", false);
  • branches/eam_branches/gDiff.20200202/src/gDiffCamera.c

    r41260 r41413  
     1/** @file gDiffCamera.c
     2 *  @brief: image subtraction based on matched Gaussians profiles
     3 *  @author Eugene Magnier @ IfA
     4 *  @version $Revision: 1.12 $
     5 *  @date $Date: 2020-01-25 $
     6 *  Copyright 2020 Institute for Astronomy, University of Hawaii
     7 */
     8
    19#include "gDiff.h"
    210
     
    207215    }
    208216
     217    /* in ppSub, there is a recipe value 'PHOTOMETRY' which is supposed to enable / disable photometry analysis.
     218       This is set to FALSE in the recipes by default.
     219
     220       In the script diff_skycell.pl, there is a global variable 'do_photom' which is TRUE by default,
     221       but is set to FALSE if a) the reduction class is NOCONVDIFF or b) the run_state is not new (i.e, 'update' mode)
     222       (the script changes the call in an inconsistent way, if do_photom is true, it appends -photometry, but
     223       if it is false, it appends -Db PHOTOMETRY FALSE.  this could be made consistent.
     224       (also note that it is bad practice to have hard-wired reduction classes)
     225
     226       But both of these situations are ignored because the test in ppSubCamera for data->photometry
     227       (set if PHOTOMETRY is true or -photometry is supplied) is ignored and the photometric analysis
     228       is always performed.
     229
     230     */
     231
    209232    // Now that the camera has been determined, we can read the recipe
    210233    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, GDIFF_RECIPE); // Recipe for ppSim
     
    380403    }
    381404
    382     // psPhot input
    383     if (data->photometry || 1) {
    384         psphotModelClassInit();        // load implementation-specific models
    385 
    386         pmFPAfile *psphot = pmFPAfileDefineFromFPA(config, output->fpa, 1, 1, "PSPHOT.INPUT");
    387         if (!psphot) {
    388             psError(psErrorCodeLast(), false, "Failed to build FPA from PSPHOT.INPUT");
    389             return false;
    390         }
    391         if (psphot->type != PM_FPA_FILE_IMAGE) {
    392             psError(psErrorCodeLast(), true, "PSPHOT.INPUT is not of type IMAGE");
    393             return false;
    394         }
    395         pmFPAfileActivate(config->files, false, "PSPHOT.INPUT");
    396 
    397         // Internal file for getting the PSF from the minuend
    398         pmFPAfile *psf = pmFPAfileDefineFromFPA(config, output->fpa, 1, 1, "PSPHOT.PSF.LOAD");
    399         if (!psf) {
    400             psError(psErrorCodeLast(), false, "Failed to build FPA from PSPHOT.PSF.LOAD");
    401             return false;
    402         }
    403         if (psf->type != PM_FPA_FILE_PSF) {
    404             psError(psErrorCodeLast(), true, "PSPHOT.PSF.LOAD is not of type PSF");
    405             return false;
    406         }
    407         pmFPAfileActivate(config->files, false, "PSPHOT.PSF.LOAD");
    408 
    409         if (!psphotDefineFiles(config, psphot)) {
    410             psError(psErrorCodeLast(), false, "Unable to set up psphot files.");
    411             return false;
    412         }
    413 
    414         // Deactivate psphot output sources --- we want to define output source files of our own
    415         pmFPAfile *psphotOutput = pmFPAfileSelectSingle(config->files, "PSPHOT.OUTPUT", 0);
    416         psphotOutput->save = false;  // this one should NOT be set
    417 
    418         pmFPAfile *outSources = defineOutputFile(config, output, false, "GDIFF.OUTPUT.SOURCES",
    419                                                  PM_FPA_FILE_CMF);
    420         if (!outSources) {
    421             psError(psErrorCodeLast(), false, "Unable to set up output source file.");
    422             return false;
    423         }
    424         checkFileruleFileSave(outSources, config);
    425 
    426         if (data->inverse) {
    427             pmFPAfile *invSources = defineOutputFile(config, inverse, false, "GDIFF.INVERSE.SOURCES",
    428                                                      PM_FPA_FILE_CMF);
    429             if (!invSources) {
    430                 psError(psErrorCodeLast(), false, "Unable to set up inverse source file.");
    431                 return false;
    432             }
    433             checkFileruleFileSave(invSources, config);
    434         }
    435 
    436         // files need to do the forced photometry on the positions of sources in the positive images
    437         if (data->forcedPhot1) {
    438             // this pmFPAfile is used to carry sources detected in the positive image #1
    439             pmFPAfile *posSources1 = defineOutputFile(config, input, true, "GDIFF.POS1.SOURCES", PM_FPA_FILE_CMF);
    440             if (!posSources1) {
    441                 psError(psErrorCodeLast(), false, "Unable to set up forced source file.");
    442                 return false;
    443             }
    444             checkFileruleFileSave(posSources1, config);
    445 
    446             // this pmFPAfile is used to carry sources detected in the diff image @ the positions from positive image #1
    447             pmFPAfile *frcSources1 = defineOutputFile(config, input, true, "GDIFF.FORCED1.SOURCES", PM_FPA_FILE_CMF);
    448             if (!frcSources1) {
    449                 psError(psErrorCodeLast(), false, "Unable to set up forced source file.");
    450                 return false;
    451             }
    452             checkFileruleFileSave(frcSources1, config);
     405    /* in ppSub, there is a recipe value 'PHOTOMETRY' which is supposed to enable / disable photometry analysis.
     406       This is set to FALSE in the recipes by default.
     407
     408       In the script diff_skycell.pl, there is a global variable 'do_photom' which is TRUE by default,
     409       but is set to FALSE if a) the reduction class is NOCONVDIFF or b) the run_state is not new (i.e, 'update' mode)
     410       (the script changes the call in an inconsistent way, if do_photom is true, it appends -photometry, but
     411       if it is false, it appends -Db PHOTOMETRY FALSE.  this could be made consistent.
     412       (also note that it is bad practice to have hard-wired reduction classes)
     413
     414       But both of these situations are ignored because the test in ppSubCamera for data->photometry
     415       (set if PHOTOMETRY is true or -photometry is supplied) is ignored and the photometric analysis
     416       is always performed.
     417
     418       Should we rationalize this?
     419     */
     420
     421    // define the psphot I/O files
     422    psphotModelClassInit();        // load implementation-specific models
     423
     424    pmFPAfile *psphot = pmFPAfileDefineFromFPA(config, output->fpa, 1, 1, "PSPHOT.INPUT");
     425    if (!psphot) {
     426        psError(psErrorCodeLast(), false, "Failed to build FPA from PSPHOT.INPUT");
     427        return false;
     428    }
     429    if (psphot->type != PM_FPA_FILE_IMAGE) {
     430        psError(psErrorCodeLast(), true, "PSPHOT.INPUT is not of type IMAGE");
     431        return false;
     432    }
     433    pmFPAfileActivate(config->files, false, "PSPHOT.INPUT");
     434
     435    // Internal file for getting the PSF from the minuend
     436    pmFPAfile *psf = pmFPAfileDefineFromFPA(config, output->fpa, 1, 1, "PSPHOT.PSF.LOAD");
     437    if (!psf) {
     438        psError(psErrorCodeLast(), false, "Failed to build FPA from PSPHOT.PSF.LOAD");
     439        return false;
     440    }
     441    if (psf->type != PM_FPA_FILE_PSF) {
     442        psError(psErrorCodeLast(), true, "PSPHOT.PSF.LOAD is not of type PSF");
     443        return false;
     444    }
     445    pmFPAfileActivate(config->files, false, "PSPHOT.PSF.LOAD");
     446
     447    if (!psphotDefineFiles(config, psphot)) {
     448        psError(psErrorCodeLast(), false, "Unable to set up psphot files.");
     449        return false;
     450    }
     451
     452    // Deactivate psphot output sources --- we want to define output source files of our own
     453    pmFPAfile *psphotOutput = pmFPAfileSelectSingle(config->files, "PSPHOT.OUTPUT", 0);
     454    psphotOutput->save = false;  // this one should NOT be set
     455
     456    pmFPAfile *outSources = defineOutputFile(config, output, false, "GDIFF.OUTPUT.SOURCES",
     457                                             PM_FPA_FILE_CMF);
     458    if (!outSources) {
     459        psError(psErrorCodeLast(), false, "Unable to set up output source file.");
     460        return false;
     461    }
     462    checkFileruleFileSave(outSources, config);
     463
     464    if (data->inverse) {
     465        pmFPAfile *invSources = defineOutputFile(config, inverse, false, "GDIFF.INVERSE.SOURCES",
     466                                                 PM_FPA_FILE_CMF);
     467        if (!invSources) {
     468            psError(psErrorCodeLast(), false, "Unable to set up inverse source file.");
     469            return false;
    453470        }
    454 
    455         if (data->forcedPhot2) {
    456             // this pmFPAfile is used to carry sources detected in the positive image #2
    457             pmFPAfile *posSources2 = defineOutputFile(config, ref, true, "GDIFF.POS2.SOURCES", PM_FPA_FILE_CMF);
    458             if (!posSources2) {
    459                 psError(psErrorCodeLast(), false, "Unable to set up forced source file.");
    460                 return false;
    461             }
    462             checkFileruleFileSave(posSources2, config);
    463 
    464             // this pmFPAfile is used to carry sources detected in the diff image @ the positions from positive image #2
    465             pmFPAfile *frcSources2 = defineOutputFile(config, ref, true, "GDIFF.FORCED2.SOURCES", PM_FPA_FILE_CMF);
    466             if (!frcSources2) {
    467                 psError(psErrorCodeLast(), false, "Unable to set up forced source file.");
    468                 return false;
    469             }
    470             checkFileruleFileSave(frcSources2, config);
    471         }
     471        checkFileruleFileSave(invSources, config);
     472    }
     473
     474    // files need to do the forced photometry on the positions of sources in the positive images
     475    if (data->forcedPhot1) {
     476        // this pmFPAfile is used to carry sources detected in the positive image #1
     477        pmFPAfile *posSources1 = defineOutputFile(config, input, true, "GDIFF.POS1.SOURCES", PM_FPA_FILE_CMF);
     478        if (!posSources1) {
     479            psError(psErrorCodeLast(), false, "Unable to set up forced source file.");
     480            return false;
     481        }
     482        checkFileruleFileSave(posSources1, config);
     483
     484        // this pmFPAfile is used to carry sources detected in the diff image @ the positions from positive image #1
     485        pmFPAfile *frcSources1 = defineOutputFile(config, input, true, "GDIFF.FORCED1.SOURCES", PM_FPA_FILE_CMF);
     486        if (!frcSources1) {
     487            psError(psErrorCodeLast(), false, "Unable to set up forced source file.");
     488            return false;
     489        }
     490        checkFileruleFileSave(frcSources1, config);
     491    }
     492
     493    if (data->forcedPhot2) {
     494        // this pmFPAfile is used to carry sources detected in the positive image #2
     495        pmFPAfile *posSources2 = defineOutputFile(config, ref, true, "GDIFF.POS2.SOURCES", PM_FPA_FILE_CMF);
     496        if (!posSources2) {
     497            psError(psErrorCodeLast(), false, "Unable to set up forced source file.");
     498            return false;
     499        }
     500        checkFileruleFileSave(posSources2, config);
     501
     502        // this pmFPAfile is used to carry sources detected in the diff image @ the positions from positive image #2
     503        pmFPAfile *frcSources2 = defineOutputFile(config, ref, true, "GDIFF.FORCED2.SOURCES", PM_FPA_FILE_CMF);
     504        if (!frcSources2) {
     505            psError(psErrorCodeLast(), false, "Unable to set up forced source file.");
     506            return false;
     507        }
     508        checkFileruleFileSave(frcSources2, config);
    472509    }
    473510
  • branches/eam_branches/gDiff.20200202/src/gDiffData.c

    r41260 r41413  
     1/** @file gDiffData.c
     2 *  @brief: image subtraction based on matched Gaussians profiles
     3 *  @author Eugene Magnier @ IfA
     4 *  @version $Revision: 1.12 $
     5 *  @date $Date: 2020-01-25 $
     6 *  Copyright 2020 Institute for Astronomy, University of Hawaii
     7 */
     8
    19#include "gDiff.h"
    210
  • branches/eam_branches/gDiff.20200202/src/gDiffDefineOutput.c

    r41260 r41413  
    11/** @file gDiffDefineOutput.c
    2  *
    3  *  @brief
    4  *
    5  *  @ingroup gDiff
    6  *
    7  *  @author IfA
    8  *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-18 00:31:20 $
    10  *  Copyright 2009 Institute for Astronomy, University of Hawaii
     2 *  @brief: image subtraction based on matched Gaussians profiles
     3 *  @author Eugene Magnier @ IfA
     4 *  @version $Revision: 1.12 $
     5 *  @date $Date: 2020-01-25 $
     6 *  Copyright 2020 Institute for Astronomy, University of Hawaii
    117 */
    12 
    13 #ifdef HAVE_CONFIG_H
    14 #include <config.h>
    15 #endif
    16 
    17 #include <stdio.h>
    18 #include <libgen.h>
    19 #include <pslib.h>
    20 #include <psmodules.h>
    218
    229#include "gDiff.h"
  • branches/eam_branches/gDiff.20200202/src/gDiffErrorCodes.dat

    r41261 r41413  
    33#
    44BASE = 14000            First value we use; lower values belong to psLib
    5 UNKNOWN                 Unknown ppSub error code
     5UNKNOWN                 Unknown gDiff error code
    66NOT_IMPLEMENTED         Desired feature is not yet implemented
    77ARGUMENTS               Incorrect arguments
  • branches/eam_branches/gDiff.20200202/src/gDiffExit.c

    r41260 r41413  
    1 #include <stdio.h>
    2 #include <pslib.h>
    3 #include <psphotErrorCodes.h>
     1/** @file gDiffExit.c
     2 *  @brief: image subtraction based on matched Gaussians profiles
     3 *  @author Eugene Magnier @ IfA
     4 *  @version $Revision: 1.12 $
     5 *  @date $Date: 2020-01-25 $
     6 *  Copyright 2020 Institute for Astronomy, University of Hawaii
     7 */
    48
    59#include "gDiff.h"
     
    1620    // not include the possible psphot values
    1721
    18     // for now, to get around this, we just use an int for the switch
     22    // thus, we cannot do this:
     23    // psErrorCode errorCode = psErrorCodeLast(); // Error code
    1924
    20     // psErrorCode errorCode = psErrorCodeLast(); // Error code
     25    // for now, to get around this, we just use an int for the switch:
    2126    int errorCode = psErrorCodeLast(); // Error code
     27
    2228    if (errorCode != PS_ERR_NONE) {
    2329        pmFPAfileFreeSetStrict(false);
  • branches/eam_branches/gDiff.20200202/src/gDiffFiles.c

    r41260 r41413  
    1 #ifdef HAVE_CONFIG_H
    2 #include <config.h>
    3 #endif
    4 
    5 #include <stdio.h>
    6 #include <pslib.h>
    7 #include <psmodules.h>
    8 #include <psphot.h>
     1/** @file gDiffFiles.c
     2 *  @brief: image subtraction based on matched Gaussians profiles
     3 *  @author Eugene Magnier @ IfA
     4 *  @version $Revision: 1.12 $
     5 *  @date $Date: 2020-01-25 $
     6 *  Copyright 2020 Institute for Astronomy, University of Hawaii
     7 */
    98
    109#include "gDiff.h"
     
    1211
    1312// Input files
    14 static const char *inputFiles[] = { "GDIFF.INPUT", "GDIFF.INPUT.MASK", "GDIFF.INPUT.VARIANCE",
    15                                     "GDIFF.INPUT.SOURCES", "GDIFF.REF", "GDIFF.REF.MASK",
    16                                     "GDIFF.REF.VARIANCE", "GDIFF.REF.SOURCES", NULL };
     13static const char *inputFiles[] = { "GDIFF.INPUT",
     14                                    "GDIFF.INPUT.MASK",
     15                                    "GDIFF.INPUT.VARIANCE",
     16                                    "GDIFF.INPUT.SOURCES",
     17                                    "GDIFF.REF",
     18                                    "GDIFF.REF.MASK",
     19                                    "GDIFF.REF.VARIANCE",
     20                                    "GDIFF.REF.SOURCES", NULL };
    1721
    1822// Convolved files
    1923static const char *convFiles[] = { "GDIFF.INPUT.CONV", "GDIFF.INPUT.CONV.MASK", "GDIFF.INPUT.CONV.VARIANCE",
    20                                    "GDIFF.REF.CONV", "GDIFF.REF.CONV.MASK", "GDIFF.REF.CONV.VARIANCE",
     24                                   "GDIFF.REF.CONV",   "GDIFF.REF.CONV.MASK",  "GDIFF.REF.CONV.VARIANCE",
    2125                                   NULL };
    2226
     
    234238    return true;
    235239}
     240
  • branches/eam_branches/gDiff.20200202/src/gDiffLoop.c

    r41260 r41413  
     1/** @file gDiffLoop.c
     2 *  @brief: image subtraction based on matched Gaussians profiles
     3 *  @author Eugene Magnier @ IfA
     4 *  @version $Revision: 1.12 $
     5 *  @date $Date: 2020-01-25 $
     6 *  Copyright 2020 Institute for Astronomy, University of Hawaii
     7 */
     8
    19#include "gDiff.h"
     10
     11# define HDUCHECK(A)                                                    \
     12    fprintf (stderr, "%s\n", A);                                        \
     13    listhdu (config, "GDIFF.OUTPUT"); listhdu (config, "PSPHOT.INPUT");
     14
     15// internal functions:
     16bool dumpout(pmConfig *config, char *name);
     17bool listhdu(pmConfig *config, char *name);
     18
     19bool gDiffLoop(gDiffData *data)
     20{
     21    bool mdok = false;
     22    bool success = true;
     23
     24    psAssert(data, "Require processing data");
     25    pmConfig *config = data->config;    // Configuration
     26    psAssert(config, "Require configuration.");
     27
     28    HDUCHECK ("0a: start gDiffLoop");
     29
     30    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, GDIFF_RECIPE); // Recipe for ppSim
     31    psAssert(recipe, "We checked this earlier, so it should be here.");
     32    bool noConvolve = psMetadataLookupBool(&mdok, recipe, "NOCONVOLVE"); // Do not use convolved images.
     33
     34    pmConfigCamerasCull(config, NULL);
     35    pmConfigRecipesCull(config, "GDIFF,PPSTATS,PSPHOT,PSASTRO,MASKS,JPEG");
     36
     37    pmFPAfile *input     = psMetadataLookupPtr(NULL, config->files, "GDIFF.INPUT");
     38    pmFPAfile *reference = psMetadataLookupPtr(NULL, config->files, "GDIFF.REF");
     39    pmFPAfile *output    = psMetadataLookupPtr(NULL, config->files, "GDIFF.OUTPUT");
     40    psAssert(input && reference && output, "Require files");
     41
     42    // load input files : we are required to load a complete set of files that describe the input.
     43    // we are not able to process subdivisions of an image (e.g., cells vs chips)
     44    if (!gDiffFilesIterateDown(config, GDIFF_FILES_INPUT | GDIFF_FILES_CONV)) {
     45        psError(GDIFF_ERR_IO, false, "Unable to load files.");
     46        return false;
     47    }
     48
     49    psTimerStart("GDIFF_MATCH");
     50
     51    HDUCHECK ("0b: gDiffSetMasks");
     52
     53    if (!gDiffSetMasks(config)) {
     54        psError(psErrorCodeLast(), false, "Unable to set masks.");
     55        return false;
     56    }
     57
     58    HDUCHECK ("0b.1: gDiffInputDetections");
     59    if (data->forcedPhot1) {
     60        bool foundDetections = false;
     61        if (!gDiffInputDetections(&foundDetections, "GDIFF.POS1.SOURCES", "GDIFF.INPUT", data)) {
     62            psError(psErrorCodeLast(), false, "Unable to measure positive detections (1)");
     63            success = false;
     64            goto ESCAPE;
     65        }
     66        // if nothing was found, don't bother doing the forced photometry below
     67        if (!foundDetections) {
     68            psWarning ("no sources found in positive image 1, skipping forced photometry");
     69            data->forcedPhot1 = false;
     70        }
     71    }
     72    HDUCHECK ("0b.2: gDiffInputDetections");
     73    if (data->forcedPhot2) {
     74        // XXX why is this only done for Phot2 and not Phot1 ??
     75
     76        // Change the recipe to use a higher nsigma limit and quit after pass1
     77        psMetadata *psphotRecipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
     78
     79        psF32 nsigma_peak_save = psMetadataLookupF32 (NULL, psphotRecipe, "PEAKS_NSIGMA_LIMIT");
     80        char *breakPt_save =  psMetadataLookupStr (NULL, psphotRecipe, "BREAK_POINT");
     81
     82        psF32 pos2_nsigma_peak = psMetadataLookupF32 (&mdok, psphotRecipe, "PEAKS_POS2_NSIGMA_LIMIT");
     83        if (!mdok) {
     84            psWarning("PEAKS_POS2_NSIGMA_LIMIT not found in psphotRecipe. Will use 25.\n");
     85            pos2_nsigma_peak = 25.;
     86        }
     87        psMetadataAddF32(psphotRecipe, PS_LIST_TAIL, "PEAKS_NSIGMA_LIMIT", PS_META_REPLACE, "", pos2_nsigma_peak);
     88        psMetadataAddStr(psphotRecipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", "PASS1");
     89
     90        bool foundDetections = false;
     91        if (!gDiffInputDetections(&foundDetections, "GDIFF.POS2.SOURCES", "GDIFF.REF", data)) {
     92            psError(psErrorCodeLast(), false, "Unable to measure positive detections (2)");
     93            psMetadataAddF32(psphotRecipe, PS_LIST_TAIL, "PEAKS_NSIGMA_LIMIT", PS_META_REPLACE, "", nsigma_peak_save);
     94            psMetadataAddStr(psphotRecipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", breakPt_save);
     95            success = false;
     96            goto ESCAPE;
     97        }
     98        psMetadataAddF32(psphotRecipe, PS_LIST_TAIL, "PEAKS_NSIGMA_LIMIT", PS_META_REPLACE, "", nsigma_peak_save);
     99        psMetadataAddStr(psphotRecipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", breakPt_save);
     100        // if nothing was found, don't bother doing the forced photometry below
     101        if (!foundDetections) {
     102            psWarning ("no sources found in positive image 2, skipping forced photometry");
     103            data->forcedPhot2 = false;
     104        }
     105    }
     106
     107    // XXX if it exists, use the POS1, POS2 successs for the FWHMs
     108
     109    /**** The PSF matching process is performed in this function below. ****
     110     **** This is the core of the PSF-matched image subtraction code.   ****/
     111
     112    HDUCHECK ("1: before gDiffMatchPSFs");
     113
     114    if (!gDiffMatchPSFs(data)) {
     115        psError(psErrorCodeLast(), false, "Unable to match PSFs.");
     116        success = false;
     117        goto ESCAPE;
     118    }
     119    if (data->quality) {
     120        // Can't do anything at all
     121        success = false;
     122        goto ESCAPE;
     123    }
     124
     125    // the functions in this block are not critical for initial testing, so block them out for now
     126
     127    // generate the residual stamp grid for visualization
     128    if (!gDiffResidualSampleJpeg(config)) {
     129        psError(psErrorCodeLast(), false, "Unable to update.");
     130        success = false;
     131        goto ESCAPE;
     132    }
     133
     134    // XXX add in a positive image detection step here (if needed)
     135    // XXX EAM what does this mean???
     136
     137    psMetadataAddF32(data->stats, PS_LIST_TAIL, "TIME_MATCH", 0, "Time to match PSFs", psTimerClear("GDIFF_MATCH"));
     138
     139    // Close input files (freeing up space) : for the noConvolve case, we cannot close the
     140    // inputs since we will use them for subtraction below.  wait until later to do the work
     141    if (!noConvolve) {
     142        if (!gDiffFilesIterateUp(config, GDIFF_FILES_INPUT)) {
     143            psError(GDIFF_ERR_IO, false, "Unable to close input files.");
     144            success = false;
     145            goto ESCAPE;
     146        }
     147    }
     148    if (!gDiffLowThreshold(data)) {
     149        psError(psErrorCodeLast(), false, "Unable to threshold images.");
     150        success = false;
     151        goto ESCAPE;
     152    }
     153
     154    HDUCHECK ("2: after gDiffLowThreshold");
     155
     156    // Set up subtraction files
     157    if (!gDiffFilesIterateDown(config, GDIFF_FILES_SUB)) {
     158        psError(GDIFF_ERR_IO, false, "Unable to set up subtraction files.");
     159        success = false;
     160        goto ESCAPE;
     161    }
     162
     163    if (!gDiffDefineOutput("GDIFF.OUTPUT", config)) {
     164        psError(psErrorCodeLast(), false, "Unable to define output.");
     165        success = false;
     166        goto ESCAPE;
     167    }
     168
     169    HDUCHECK ("2a: gDiffMakePSF");
     170
     171    // the functions in this block are not critical for initial testing, so block them out for now
     172    if (!data->quality && !gDiffMakePSF(data)) {
     173        psError(psErrorCodeLast(), false, "Unable to generate PSF.");
     174        success = false;
     175        goto ESCAPE;
     176    }
     177
     178    // Now we've got a PSF, blow away detections in case they're confused with real output detections
     179    {
     180        pmFPAview *view = gDiffViewReadout(); // View to readout
     181        pmReadout *out = pmFPAfileThisReadout(config->files, view, "GDIFF.OUTPUT");
     182        if (psMetadataLookup(out->analysis, "PSPHOT.DETECTIONS")) {
     183            psMetadataRemoveKey(out->analysis, "PSPHOT.DETECTIONS");
     184        }
     185        psFree(view);
     186    }
     187
     188    HDUCHECK ("3: before gDiffReadoutSubtract");
     189
     190    if (!gDiffReadoutSubtract(config)) {
     191        psError(psErrorCodeLast(), false, "Unable to subtract images.");
     192        success = false;
     193        goto ESCAPE;
     194    }
     195
     196    // Close convolved files
     197    if (!gDiffFilesIterateUp(config, GDIFF_FILES_PSF | GDIFF_FILES_CONV)) {
     198        psError(GDIFF_ERR_IO, false, "Unable to close input files.");
     199        success = false;
     200        goto ESCAPE;
     201    }
     202    // Close unconvolved input files if we did the noConvolve subtraction
     203    // (otherwise they are closed above)
     204    if (noConvolve && !gDiffFilesIterateUp(config, GDIFF_FILES_INPUT)) {
     205        psError(GDIFF_ERR_IO, false, "Unable to close input files.");
     206        success = false;
     207        goto ESCAPE;
     208    }
     209   
     210    HDUCHECK ("4: before gDiffBackground");
     211
     212    // Higher order background subtraction using psphot
     213    if (!gDiffBackground(config)) {
     214        psError(psErrorCodeLast(), false, "Unable to subtract background.");
     215        success = false;
     216        goto ESCAPE;
     217    }
     218    // dumpout(config, "diff.2b.fits");
     219
     220    HDUCHECK ("5: before gDiffVarianceRescale");
     221
     222    // Perform Variance correction (rescale within a modest range)
     223    if (!gDiffVarianceRescale(config, data)) {
     224        psError(psErrorCodeLast(), false, "Unable to rescale variance.");
     225        success = false;
     226        goto ESCAPE;
     227    }
     228    // dumpout(config, "diff.2c.fits");
     229
     230    if (data->quality) {
     231        // Done all we can do up to this point
     232        success = false;
     233        goto ESCAPE;
     234    }
     235
     236    if (!gDiffFilesIterateDown(config, GDIFF_FILES_PHOT_SUB)) {
     237        psError(GDIFF_ERR_IO, false, "Unable to set up photometry files.");
     238        success = false;
     239        goto ESCAPE;
     240    }
     241
     242    HDUCHECK ("6: before gDiffReadoutPhotometry");
     243
     244    if (!data->quality && !gDiffReadoutPhotometry("GDIFF.OUTPUT", data)) {
     245        psError(psErrorCodeLast(), false, "Unable to perform photometry.");
     246        success = false;
     247        goto ESCAPE;
     248    }
     249    // dumpout(config, "diff.3.fits");
     250
     251    HDUCHECK ("7: before gDiffReadoutForcedPhot (1)");
     252
     253    // forced photometry for positive image 1
     254    if (data->forcedPhot1 && !data->quality) {
     255        if (!gDiffReadoutForcedPhot("GDIFF.FORCED1.SOURCES", "GDIFF.OUTPUT", "GDIFF.POS1.SOURCES", data)) {
     256            psError(psErrorCodeLast(), false, "Unable to perform photometry.");
     257            success = false;
     258            goto ESCAPE;
     259        }
     260    }
     261
     262    // forced photometry for positive image 2
     263    if (data->forcedPhot2 && !data->quality) {
     264        if (!gDiffReadoutForcedPhot("GDIFF.FORCED2.SOURCES", "GDIFF.OUTPUT", "GDIFF.POS2.SOURCES", data)) {
     265            psError(psErrorCodeLast(), false, "Unable to perform photometry.");
     266            success = false;
     267            goto ESCAPE;
     268        }
     269    }
     270
     271    if (!gDiffFilesIterateUp(config, GDIFF_FILES_PHOT_SUB)) {
     272        psError(GDIFF_ERR_IO, false, "Unable to set up photometry files.");
     273        success = false;
     274        goto ESCAPE;
     275    }
     276
     277    // Perform statistics on the cell
     278    if (!gDiffReadoutStats(data)) {
     279        psError(psErrorCodeLast(), false, "Unable to collect statistics");
     280        success = false;
     281        goto ESCAPE;
     282    }
     283    // Do Mask Stats
     284    {
     285        pmFPAview *view = gDiffViewReadout(); // View to readout
     286        if (!gDiffMaskStats(config, view,data->stats)) {
     287            psError(psErrorCodeLast(), false, "Unable to generate mask statistics");
     288            success = false;
     289            psFree(view);
     290            goto ESCAPE;
     291        }
     292        psFree(view);
     293    }
     294    // dumpout(config, "diff.4.fits");
     295   
     296    // generate the binned image used to write the jpeg
     297    if (!gDiffReadoutJpeg(config)) {
     298        psError(psErrorCodeLast(), false, "Unable to update.");
     299        success = false;
     300        goto ESCAPE;
     301    }
     302
     303    if (data->inverse) {
     304        // Set up inverse subtraction files
     305        if (!gDiffFilesIterateDown(config, GDIFF_FILES_INV)) {
     306            psError(GDIFF_ERR_IO, false, "Unable to set up inverse files.");
     307            success = false;
     308            goto ESCAPE;
     309        }
     310
     311        if (data->inverse && !gDiffDefineOutput("GDIFF.INVERSE", config)) {
     312            psError(psErrorCodeLast(), false, "Unable to define inverse.");
     313            success = false;
     314            goto ESCAPE;
     315        }
     316
     317        if (!gDiffReadoutInverse(config)) {
     318            psError(psErrorCodeLast(), false, "Unable to invert images.");
     319            success = false;
     320            goto ESCAPE;
     321        }
     322
     323        //MEH -- need to also clear out detections or inv.cmf corrupted on update
     324        pmFPAview *viewinv = gDiffViewReadout(); // View to readout
     325        pmReadout *inv = pmFPAfileThisReadout(config->files, viewinv, "GDIFF.INVERSE");
     326        if (psMetadataLookup(inv->analysis, "PSPHOT.DETECTIONS")) {
     327            psMetadataRemoveKey(inv->analysis, "PSPHOT.DETECTIONS");
     328        }
     329        psFree(viewinv);
     330
     331        // Close subtraction files and open inverse photometry files
     332        if (!gDiffFilesIterateUp(config, GDIFF_FILES_SUB)) {
     333            psError(GDIFF_ERR_IO, false, "Unable to close subtraction files.");
     334            success = false;
     335            goto ESCAPE;
     336        }
     337
     338        if (!gDiffFilesIterateDown(config, GDIFF_FILES_PHOT_INV)) {
     339            psError(GDIFF_ERR_IO, false, "Unable to set up inverse files.");
     340            success = false;
     341            goto ESCAPE;
     342        }
     343
     344        if (!data->quality && !gDiffReadoutPhotometry("GDIFF.INVERSE", data)) {
     345            psError(psErrorCodeLast(), false, "Unable to perform photometry.");
     346            success = false;
     347            goto ESCAPE;
     348        }
     349
     350        // Close inverse subtraction files
     351        if (!gDiffFilesIterateUp(config, GDIFF_FILES_INV | GDIFF_FILES_PHOT_INV)) {
     352            psError(GDIFF_ERR_IO, false, "Unable to close subtraction files.");
     353            success = false;
     354            goto ESCAPE;
     355        }
     356    } else {
     357        // Close subtraction files
     358        // dumpout(config, "diff.5.fits");
     359        if (!gDiffFilesIterateUp(config, GDIFF_FILES_SUB)) {
     360            psError(GDIFF_ERR_IO, false, "Unable to close subtraction files.");
     361            success = false;
     362            goto ESCAPE;
     363        }
     364    }
     365
     366 ESCAPE:
     367    pmFPAfileDropInternal(config->files, "PSPHOT.BACKGND");
     368    pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL");
     369    pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL.STDEV");
     370
     371    return success;
     372}
     373
     374/**** internal functions ****/
    2375
    3376bool dumpout(pmConfig *config, char *name)
     
    10383}
    11384
    12 
    13 bool gDiffLoop(gDiffData *data)
     385bool listhdu(pmConfig *config, char *name)
    14386{
    15     bool mdok = false;
    16     bool success = true;
    17 
    18     psAssert(data, "Require processing data");
    19     pmConfig *config = data->config;    // Configuration
    20     psAssert(config, "Require configuration.");
    21 
    22     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, GDIFF_RECIPE); // Recipe for ppSim
    23     psAssert(recipe, "We checked this earlier, so it should be here.");
    24     bool noConvolve = psMetadataLookupBool(&mdok, recipe, "NOCONVOLVE"); // Do not use convolved images.
    25 
    26     pmConfigCamerasCull(config, NULL);
    27     pmConfigRecipesCull(config, "GDIFF,PPSTATS,PSPHOT,PSASTRO,MASKS,JPEG");
    28 
    29     pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "GDIFF.INPUT");
    30     pmFPAfile *reference = psMetadataLookupPtr(NULL, config->files, "GDIFF.REF");
    31     pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "GDIFF.OUTPUT");
    32     psAssert(input && reference && output, "Require files");
    33 
    34     if (!gDiffFilesIterateDown(config, GDIFF_FILES_INPUT | GDIFF_FILES_CONV)) {
    35         psError(GDIFF_ERR_IO, false, "Unable to load files.");
    36         return false;
    37     }
    38 
    39     psTimerStart("GDIFF_MATCH");
    40 
    41     if (!gDiffSetMasks(config)) {
    42         psError(psErrorCodeLast(), false, "Unable to set masks.");
    43         return false;
    44     }
    45 
    46     if (data->forcedPhot1) {
    47         bool foundDetections = false;
    48         if (!gDiffInputDetections(&foundDetections, "GDIFF.POS1.SOURCES", "GDIFF.INPUT", data)) {
    49             psError(psErrorCodeLast(), false, "Unable to measure positive detections (1)");
    50             success = false;
    51             goto ESCAPE;
    52         }
    53         // if nothing was found, don't bother doing the forced photometry below
    54         if (!foundDetections) {
    55             psWarning ("no sources found in positive image 1, skipping forced photometry");
    56             data->forcedPhot1 = false;
    57         }
    58     }
    59     if (data->forcedPhot2) {
    60         // Change the recipe to use a higher nsigma limit and quit after pass1
    61         psMetadata *psphotRecipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
    62 
    63         psF32 nsigma_peak_save = psMetadataLookupF32 (NULL, psphotRecipe, "PEAKS_NSIGMA_LIMIT");
    64         char *breakPt_save =  psMetadataLookupStr (NULL, psphotRecipe, "BREAK_POINT");
    65 
    66         psF32 pos2_nsigma_peak = psMetadataLookupF32 (&mdok, psphotRecipe, "PEAKS_POS2_NSIGMA_LIMIT");
    67         if (!mdok) {
    68             psWarning("PEAKS_POS2_NSIGMA_LIMIT not found in psphotRecipe. Will use 25.\n");
    69             pos2_nsigma_peak = 25.;
    70         }
    71         psMetadataAddF32(psphotRecipe, PS_LIST_TAIL, "PEAKS_NSIGMA_LIMIT", PS_META_REPLACE, "", pos2_nsigma_peak);
    72         psMetadataAddStr(psphotRecipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", "PASS1");
    73 
    74         bool foundDetections = false;
    75         if (!gDiffInputDetections(&foundDetections, "GDIFF.POS2.SOURCES", "GDIFF.REF", data)) {
    76             psError(psErrorCodeLast(), false, "Unable to measure positive detections (2)");
    77             psMetadataAddF32(psphotRecipe, PS_LIST_TAIL, "PEAKS_NSIGMA_LIMIT", PS_META_REPLACE, "", nsigma_peak_save);
    78             psMetadataAddStr(psphotRecipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", breakPt_save);
    79             success = false;
    80             goto ESCAPE;
    81         }
    82         psMetadataAddF32(psphotRecipe, PS_LIST_TAIL, "PEAKS_NSIGMA_LIMIT", PS_META_REPLACE, "", nsigma_peak_save);
    83         psMetadataAddStr(psphotRecipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", breakPt_save);
    84         // if nothing was found, don't bother doing the forced photometry below
    85         if (!foundDetections) {
    86             psWarning ("no sources found in positive image 2, skipping forced photometry");
    87             data->forcedPhot2 = false;
    88         }
    89     }
    90 
    91     // XXX if it exists, use the POS1, POS2 successs for the FWHMs
    92     if (!gDiffMatchPSFs(data)) {
    93         psError(psErrorCodeLast(), false, "Unable to match PSFs.");
    94         success = false;
    95         goto ESCAPE;
    96     }
    97 
    98     if (data->quality) {
    99         // Can't do anything at all
    100         success = false;
    101         goto ESCAPE;
    102     }
    103     // generate the residual stamp grid for visualization
    104     if (!gDiffResidualSampleJpeg(config)) {
    105         psError(psErrorCodeLast(), false, "Unable to update.");
    106         success = false;
    107         goto ESCAPE;
    108     }
    109 
    110     // XXX add in a positive image detection step here (if needed)
    111    
    112 
    113     psMetadataAddF32(data->stats, PS_LIST_TAIL, "TIME_MATCH", 0, "Time to match PSFs",
    114                      psTimerClear("GDIFF_MATCH"));
    115 
    116     // Close input files (freeing up space) : for the noConvolve case, we cannot close the
    117     // inputs since we will use them for subtraction below.  wait until later to do the work
    118     if (!noConvolve) {
    119         if (!gDiffFilesIterateUp(config, GDIFF_FILES_INPUT)) {
    120             psError(GDIFF_ERR_IO, false, "Unable to close input files.");
    121             success = false;
    122             goto ESCAPE;
    123         }
    124     }
    125     if (!gDiffLowThreshold(data)) {
    126         psError(psErrorCodeLast(), false, "Unable to threshold images.");
    127         success = false;
    128         goto ESCAPE;
    129     }
    130 
    131     // Set up subtraction files
    132     if (!gDiffFilesIterateDown(config, GDIFF_FILES_SUB)) {
    133         psError(GDIFF_ERR_IO, false, "Unable to set up subtraction files.");
    134         success = false;
    135         goto ESCAPE;
    136     }
    137 
    138     if (!gDiffDefineOutput("GDIFF.OUTPUT", config)) {
    139         psError(psErrorCodeLast(), false, "Unable to define output.");
    140         success = false;
    141         goto ESCAPE;
    142     }
    143 
    144     if (!data->quality && !gDiffMakePSF(data)) {
    145         psError(psErrorCodeLast(), false, "Unable to generate PSF.");
    146         success = false;
    147         goto ESCAPE;
    148     }
    149 
    150     // Now we've got a PSF, blow away detections in case they're confused with real output detections
    151     {
    152         pmFPAview *view = gDiffViewReadout(); // View to readout
    153         pmReadout *out = pmFPAfileThisReadout(config->files, view, "GDIFF.OUTPUT");
    154         if (psMetadataLookup(out->analysis, "PSPHOT.DETECTIONS")) {
    155             psMetadataRemoveKey(out->analysis, "PSPHOT.DETECTIONS");
    156         }
    157         psFree(view);
    158     }
    159 
    160     if (!gDiffReadoutSubtract(config)) {
    161         psError(psErrorCodeLast(), false, "Unable to subtract images.");
    162         success = false;
    163         goto ESCAPE;
    164     }
    165     // dumpout(config, "diff.1.fits");
    166 
    167     // Close convolved files
    168     if (!gDiffFilesIterateUp(config, GDIFF_FILES_PSF | GDIFF_FILES_CONV)) {
    169         psError(GDIFF_ERR_IO, false, "Unable to close input files.");
    170         success = false;
    171         goto ESCAPE;
    172     }
    173     // dumpout(config, "diff.2a.fits");
    174     if (noConvolve) {
    175         if (!gDiffFilesIterateUp(config, GDIFF_FILES_INPUT)) {
    176             psError(GDIFF_ERR_IO, false, "Unable to close input files.");
    177             success = false;
    178             goto ESCAPE;
    179         }
    180     }
    181    
    182     // Higher order background subtraction using psphot
    183     if (!gDiffBackground(config)) {
    184         psError(psErrorCodeLast(), false, "Unable to subtract background.");
    185         success = false;
    186         goto ESCAPE;
    187     }
    188     // dumpout(config, "diff.2b.fits");
    189 
    190     // Perform Variance correction (rescale within a modest range)
    191     if (!gDiffVarianceRescale(config, data)) {
    192         psError(psErrorCodeLast(), false, "Unable to rescale variance.");
    193         success = false;
    194         goto ESCAPE;
    195     }
    196     // dumpout(config, "diff.2c.fits");
    197 
    198     if (data->quality) {
    199         // Done all we can do up to this point
    200         success = false;
    201         goto ESCAPE;
    202     }
    203 
    204     if (!gDiffFilesIterateDown(config, GDIFF_FILES_PHOT_SUB)) {
    205         psError(GDIFF_ERR_IO, false, "Unable to set up photometry files.");
    206         success = false;
    207         goto ESCAPE;
    208     }
    209 
    210     if (!data->quality && !gDiffReadoutPhotometry("GDIFF.OUTPUT", data)) {
    211         psError(psErrorCodeLast(), false, "Unable to perform photometry.");
    212         success = false;
    213         goto ESCAPE;
    214     }
    215     // dumpout(config, "diff.3.fits");
    216 
    217     // forced photometry for positive image 1
    218     if (data->forcedPhot1 && !data->quality) {
    219         if (!gDiffReadoutForcedPhot("GDIFF.FORCED1.SOURCES", "GDIFF.OUTPUT", "GDIFF.POS1.SOURCES", data)) {
    220             psError(psErrorCodeLast(), false, "Unable to perform photometry.");
    221             success = false;
    222             goto ESCAPE;
    223         }
    224     }
    225 
    226     // forced photometry for positive image 2
    227     if (data->forcedPhot2 && !data->quality) {
    228         if (!gDiffReadoutForcedPhot("GDIFF.FORCED2.SOURCES", "GDIFF.OUTPUT", "GDIFF.POS2.SOURCES", data)) {
    229             psError(psErrorCodeLast(), false, "Unable to perform photometry.");
    230             success = false;
    231             goto ESCAPE;
    232         }
    233     }
    234 
    235     if (!gDiffFilesIterateUp(config, GDIFF_FILES_PHOT_SUB)) {
    236         psError(GDIFF_ERR_IO, false, "Unable to set up photometry files.");
    237         success = false;
    238         goto ESCAPE;
    239     }
    240 
    241     // Perform statistics on the cell
    242     if (!gDiffReadoutStats(data)) {
    243         psError(psErrorCodeLast(), false, "Unable to collect statistics");
    244         success = false;
    245         goto ESCAPE;
    246     }
    247     // Do Mask Stats
    248     {
    249         pmFPAview *view = gDiffViewReadout(); // View to readout
    250         if (!gDiffMaskStats(config, view,data->stats)) {
    251             psError(psErrorCodeLast(), false, "Unable to generate mask statistics");
    252             success = false;
    253             psFree(view);
    254             goto ESCAPE;
    255         }
    256         psFree(view);
    257     }
    258     // dumpout(config, "diff.4.fits");
    259    
    260     // generate the binned image used to write the jpeg
    261     if (!gDiffReadoutJpeg(config)) {
    262         psError(psErrorCodeLast(), false, "Unable to update.");
    263         success = false;
    264         goto ESCAPE;
    265     }
    266 
    267     if (data->inverse) {
    268         // Set up inverse subtraction files
    269         if (!gDiffFilesIterateDown(config, GDIFF_FILES_INV)) {
    270             psError(GDIFF_ERR_IO, false, "Unable to set up inverse files.");
    271             success = false;
    272             goto ESCAPE;
    273         }
    274 
    275         if (data->inverse && !gDiffDefineOutput("GDIFF.INVERSE", config)) {
    276             psError(psErrorCodeLast(), false, "Unable to define inverse.");
    277             success = false;
    278             goto ESCAPE;
    279         }
    280 
    281         if (!gDiffReadoutInverse(config)) {
    282             psError(psErrorCodeLast(), false, "Unable to invert images.");
    283             success = false;
    284             goto ESCAPE;
    285         }
    286 
    287         //MEH -- need to also clear out detections or inv.cmf corrupted on update
    288         pmFPAview *viewinv = gDiffViewReadout(); // View to readout
    289         pmReadout *inv = pmFPAfileThisReadout(config->files, viewinv, "GDIFF.INVERSE");
    290         if (psMetadataLookup(inv->analysis, "PSPHOT.DETECTIONS")) {
    291             psMetadataRemoveKey(inv->analysis, "PSPHOT.DETECTIONS");
    292         }
    293         psFree(viewinv);
    294 
    295         // Close subtraction files and open inverse photometry files
    296         if (!gDiffFilesIterateUp(config, GDIFF_FILES_SUB)) {
    297             psError(GDIFF_ERR_IO, false, "Unable to close subtraction files.");
    298             success = false;
    299             goto ESCAPE;
    300         }
    301 
    302         if (!gDiffFilesIterateDown(config, GDIFF_FILES_PHOT_INV)) {
    303             psError(GDIFF_ERR_IO, false, "Unable to set up inverse files.");
    304             success = false;
    305             goto ESCAPE;
    306         }
    307 
    308         if (!data->quality && !gDiffReadoutPhotometry("GDIFF.INVERSE", data)) {
    309             psError(psErrorCodeLast(), false, "Unable to perform photometry.");
    310             success = false;
    311             goto ESCAPE;
    312         }
    313 
    314         // Close inverse subtraction files
    315         if (!gDiffFilesIterateUp(config, GDIFF_FILES_INV | GDIFF_FILES_PHOT_INV)) {
    316             psError(GDIFF_ERR_IO, false, "Unable to close subtraction files.");
    317             success = false;
    318             goto ESCAPE;
    319         }
    320     } else {
    321         // Close subtraction files
    322         // dumpout(config, "diff.5.fits");
    323         if (!gDiffFilesIterateUp(config, GDIFF_FILES_SUB)) {
    324             psError(GDIFF_ERR_IO, false, "Unable to close subtraction files.");
    325             success = false;
    326             goto ESCAPE;
    327         }
    328     }
    329 
    330  ESCAPE:
    331     pmFPAfileDropInternal(config->files, "PSPHOT.BACKGND");
    332     pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL");
    333     pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL.STDEV");
    334 
    335     return success;
     387    pmFPAfile *file = psMetadataLookupPtr(NULL, config->files, name);
     388    fprintf (stderr, "HDU: %s, %llx\n", name, (long long unsigned int) file->fpa->hdu);
     389    return true;
    336390}
  • branches/eam_branches/gDiff.20200202/src/gDiffVersion.c

    r41260 r41413  
     1/** @file gDiffVersion.c
     2 *  @brief: image subtraction based on matched Gaussians profiles
     3 *  @author Eugene Magnier @ IfA
     4 *  @version $Revision: 1.12 $
     5 *  @date $Date: 2020-01-25 $
     6 *  Copyright 2020 Institute for Astronomy, University of Hawaii
     7 */
     8
    19#include "gDiff.h"
    2 #include <ppStats.h>
    310#include "gDiffVersionDefinitions.h"
    411
Note: See TracChangeset for help on using the changeset viewer.