IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 21524


Ignore:
Timestamp:
Feb 17, 2009, 2:31:20 PM (17 years ago)
Author:
Paul Price
Message:

Cleaning up. Had trouble getting the output photometry file written out, but works now.

Location:
trunk/ppSub/src
Files:
16 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppSub/src/ppSub.c

    r21424 r21524  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-09 21:26:05 $
     8 *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-18 00:31:20 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
     12
     13#ifdef HAVE_CONFIG_H
     14#include <config.h>
     15#endif
     16
     17#include <stdio.h>
     18#include <pslib.h>
     19#include <psmodules.h>
     20#include <psphot.h>
    1221
    1322#include "ppSub.h"

  • trunk/ppSub/src/ppSub.h

    r21360 r21524  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-06 01:37:17 $
     8 *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-18 00:31:20 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
    … …  
    1414#define PP_SUB_H
    1515
    16 #ifdef HAVE_CONFIG_H
    17 #include <config.h>
    18 #endif
    19 
    20 #include <stdio.h>
    21 #include <string.h>
    2216#include <pslib.h>
    2317#include <psmodules.h>
    24 #include <psphot.h>
    25 #include <ppStats.h>
    2618
    2719/// @addtogroup ppSub
    … …  
    3224/// Setup the arguments parsing
    3325bool ppSubArgumentsSetup(int argc, char *argv[], ///< Command-line arguments
    34                         pmConfig *config    ///< Configuration
     26                        pmConfig *config    ///< Configuration
    3527    );
    3628
    … …  
    5345    );
    5446
    55 /// generate (if needed) and set or update the masks for input and reference images
    56 bool ppSubSetMasks (pmConfig *config,     ///< Configuration
    57                     const pmFPAview *view ///< View of active readout
     47/// Generate (if needed) and set or update the masks for input and reference images
     48bool ppSubSetMasks(pmConfig *config,     ///< Configuration
     49                   const pmFPAview *view ///< View of active readout
    5850    );
    5951
    60 /// Generate the PSF-matching kernel and convolve the images as needed.  Most of this function
    61 /// involves looking up the parameters in the recipe and supplying them to the function
    62 /// pmSubtractionMatch()
    63 bool ppSubMatchPSFs (pmConfig *config,    ///< Configuration
    64                      const pmFPAview *view ///< View of active readout
     52/// Generate the PSF-matching kernel and convolve the images as needed.  Most of this function involves
     53/// looking up the parameters in the recipe and supplying them to the function pmSubtractionMatch()
     54bool ppSubMatchPSFs(pmConfig *config,    ///< Configuration
     55                    const pmFPAview *view ///< View of active readout
    6556    );
    6657
    67 /// generate the output readout and pass the kernel info to the header
    68 bool ppSubDefineOutput(pmConfig *config,          ///< Configuration
    69                        const pmFPAview *view ///< View of active readout
     58/// Generate the output readout and pass the kernel info to the header
     59bool ppSubDefineOutput(pmConfig *config, ///< Configuration
     60                       const pmFPAview *view ///< View of active readout
    7061    );
    7162
    72 /// Calculate the variance factor for the output image based on the input images
    73 bool ppSubVarianceFactors(pmConfig *config,       ///< Configuration
    74                           psMetadata *stats,    ///< Statistics, for output
    75                           const pmFPAview *view ///< View of active readout
    76     );
    77 
    78 /// Photometry stage 1: measure the PSF from the minuend image
    79 bool ppSubMakePSF(pmConfig *config,       ///< Configuration
    80                   const pmFPAview *view ///< View of active readout
     63/// Photometry stage 1: measure the PSF from the minuend image
     64bool ppSubMakePSF(pmConfig *config,       ///< Configuration
     65                  const pmFPAview *view ///< View of active readout
    8166    );
    8267
    8368/// Perform the actual image subtraction, update output concepts
    84 bool ppSubReadoutSubtract(pmConfig *config,       ///< Configuration
    85                           const pmFPAview *view ///< View of active readout
     69bool ppSubReadoutSubtract(pmConfig *config,       ///< Configuration
     70                          const pmFPAview *view ///< View of active readout
    8671    );
    8772
    8873
    8974/// Photometry stage 2: find and measure sources on the subtracted image
    90 bool ppSubReadoutPhotometry(pmConfig *config,     ///< Configuration
    91                             psMetadata *stats,    ///< Statistics, for output
    92                             const pmFPAview *view ///< View of active readout
     75bool ppSubReadoutPhotometry(pmConfig *config,     ///< Configuration
     76                            psMetadata *stats,    ///< Statistics, for output
     77                            const pmFPAview *view ///< View of active readout
    9378    );
    9479
    9580/// Renormalize, update headers and generate JPEGs
    96 bool ppSubReadoutUpdate(pmConfig *config,         ///< Configuration
    97                         const pmFPAview *view ///< View of active readout
     81bool ppSubReadoutUpdate(pmConfig *config, ///< Configuration
     82                        psMetadata *stats, ///< Statistics for output, or NULL
     83                        const pmFPAview *view ///< View of active readout
    9884    );
    9985
    … …  
    10793    );
    10894
    109 /// Generate and Set the masks if needed
    110 bool ppSubSetMasks (
    111     pmConfig *config,                   ///< Configuration
    112     const pmFPAview *view               ///< view to readout
    113     );
    114 
    115 /// Renormalize readout for peak pixels
    116 bool ppSubReadoutRenormPixels (
    117     pmConfig *config,                   ///< Configuration
    118     psMetadata *recipe,                         ///< Recipe
    119     pmReadout *readout                  ///< Readout
    120     );
    121 
    122 /// Renormalize readout for photometry analysis
    123 bool ppSubReadoutRenormPhot (
    124     pmConfig *config,                   ///< Configuration
    125     psMetadata *recipe,                         ///< Recipe
    126     pmReadout *readout                  ///< Readout
    127     );
    12895
    12996// Copy every instance of a single keyword from one metadata to another

  • trunk/ppSub/src/ppSubArguments.c

    r21424 r21524  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.58 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-09 21:26:05 $
     8 *  @version $Revision: 1.59 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-18 00:31:20 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
     12
     13#ifdef HAVE_CONFIG_H
     14#include <config.h>
     15#endif
     16
     17#include <stdio.h>
     18#include <pslib.h>
     19#include <psmodules.h>
     20#include <psphot.h>
    1221
    1322#include "ppSub.h"
    … …  
    240249    psMetadataAddS32(arguments, PS_LIST_TAIL, "-opt-order", 0, "Maximum order for optimum kernel search", -1);
    241250    psMetadataAddBool(arguments, PS_LIST_TAIL, "-dual", 0, "Dual convolution", false);
    242     psMetadataAddBool(arguments, PS_LIST_TAIL, "-renorm", 0, "Renormalise variance maps?", false);
    243     psMetadataAddStr(arguments, PS_LIST_TAIL, "-renorm-mean", 0,
    244                      "Statistic for mean in renormalisation", NULL);
    245     psMetadataAddStr(arguments, PS_LIST_TAIL, "-renorm-stdev", 0,
    246                      "Statistic for stdev in renormalisation", NULL);
    247     psMetadataAddF32(arguments, PS_LIST_TAIL, "-renorm-width", 0, "Gaussian width for renormalisation", NAN);
    248     psMetadataAddS32(arguments, PS_LIST_TAIL, "-renorm-num", 0, "Number of samples for renormalisation", 0);
    249251    psMetadataAddBool(arguments, PS_LIST_TAIL, "-photometry", 0, "Perform photometry?", false);
    250252    psMetadataAddS32(arguments, PS_LIST_TAIL, "-threads", 0, "Number of threads", 0);
    … …  
    352354                      psMetadataLookupBool(NULL, arguments, "-dual"));
    353355
    354     if (psMetadataLookupBool(NULL, arguments, "-renorm") ||
    355         psMetadataLookupBool(NULL, recipe, "RENORM")) {
    356         psMetadataAddBool(arguments, PS_LIST_TAIL, "RENORM", 0, "Renormalise variance maps?", true);
    357     }
    358     VALUE_ARG_RECIPE_INT("-renorm-num", "RENORM.NUM", S32, 0);
    359     VALUE_ARG_RECIPE_FLOAT("-renorm-width", "RENORM.WIDTH", F32);
    360     valueArgRecipeStr(arguments, recipe, "-renorm-mean", "RENORM.MEAN", recipe);
    361     valueArgRecipeStr(arguments, recipe, "-renorm-stdev", "RENORM.STDEV", recipe);
    362 
    363356    // Need to update this because it could have been overwritten by the camera's own recipe
    364357    if (psMetadataLookupBool(NULL, arguments, "-photometry")) {

  • trunk/ppSub/src/ppSubBackground.c

    r21360 r21524  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-06 01:37:17 $
     8 *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-18 00:31:20 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
    1212
     13#ifdef HAVE_CONFIG_H
     14#include <config.h>
     15#endif
     16
     17#include <stdio.h>
     18#include <pslib.h>
     19#include <psmodules.h>
     20#include <psphot.h>
     21
    1322#include "ppSub.h"
    1423
    15 /**
    16  * Based on ppImageSubtractBackground()
    17  */
    1824bool ppSubBackground(pmConfig *config, const pmFPAview *view)
    1925{
    20     psAssert(config, "Need configuration");
    21     psAssert(view, "Need view to chip");
     26    psAssert(config, "Require configuration");
     27    psAssert(view, "Require view");
    2228
    23     bool status; // Status of metadata lookups
     29    bool mdok; // Status of metadata lookups
    2430
    25     psMetadata *ppSubRecipe = psMetadataLookupPtr(NULL, config->recipes, PPSUB_RECIPE);
     31    psMetadata *ppSubRecipe = psMetadataLookupPtr(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSub
    2632    psAssert(ppSubRecipe, "Need PPSUB recipe");
    27 
    28     psMetadata *psphotRecipe = psMetadataLookupPtr(NULL, config->recipes, PSPHOT_RECIPE);
     33    psMetadata *psphotRecipe = psMetadataLookupPtr(NULL, config->recipes, PSPHOT_RECIPE); // Recipe for psphot
    2934    psAssert(psphotRecipe, "Need PSPHOT recipe for binning");
    3035
    31     psImageMaskType maskBad = pmConfigMaskGet("BLANK", config);
     36    psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask
    3237
    33     // select the output readout
    34     pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT");
     38    pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT"); // Output image
     39    pmReadout *modelRO = pmFPAfileThisReadout(config->files, view, "PSPHOT.BACKMDL"); // Background model
    3540
    36     // select the model readout, if it exist already; if not, generate it
    37     pmReadout *modelRO = pmFPAfileThisReadout(config->files, view, "PSPHOT.BACKMDL");
    38 
    39     // if necessary, generate the background model
     41    // Generate the background model, if required
    4042    if (!modelRO) {
    4143        // Create the background model
    … …  
    4446            return false;
    4547        }
    46         // select the model readout (should now exist)
    47         modelRO = pmFPAfileThisReadout(config->files, view, "PSPHOT.BACKMDL");
     48        // select the model readout (should now exist)
     49        modelRO = pmFPAfileThisReadout(config->files, view, "PSPHOT.BACKMDL");
    4850        if (!modelRO) {
    4951            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find background model");
    … …  
    5153        }
    5254    }
    53     psImageBinning *binning = psMetadataLookupPtr(&status, psphotRecipe,
     55    psImageBinning *binning = psMetadataLookupPtr(&mdok, psphotRecipe,
    5456                                                  "PSPHOT.BACKGROUND.BINNING"); // Binning for model
    5557    psImage *modelImage = modelRO->image; // Background model
    56 
    5758    psImage *image = outRO->image; // Image of interest
    5859    psImage *mask = outRO->mask; // Mask of interest

  • trunk/ppSub/src/ppSubCamera.c

    r21374 r21524  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-06 03:20:02 $
     8 *  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-18 00:31:20 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
     12
     13#ifdef HAVE_CONFIG_H
     14#include <config.h>
     15#endif
     16
     17#include <stdio.h>
     18#include <pslib.h>
     19#include <psmodules.h>
     20#include <psphot.h>
    1221
    1322#include "ppSub.h"
    … …  
    257266    // psPhot input
    258267    if (psMetadataLookupBool(NULL, recipe, "PHOTOMETRY")) {
    259         psphotModelClassInit ();        // load implementation-specific models
    260 
    261 
    262         // Internal-ish file for getting the PSF from the matched addition
     268        psphotModelClassInit();        // load implementation-specific models
     269
     270        // Internal-ish file for getting the PSF from the minuend
    263271        pmFPAfile *psf = pmFPAfileDefineOutputFromFile(config, output, "PSPHOT.PSF.LOAD");
    264272        if (!psf) {
    … …  
    272280        pmFPAfileActivate(config->files, false, "PSPHOT.PSF.LOAD");
    273281
    274         pmFPAfile *psphotResid = pmFPAfileDefineOutputForFormat (config, NULL, "PSPHOT.RESID", output->cameraName, output->formatName);
    275         if (!psphotResid) {
    276             psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.RESID");
    277             return false;
    278         }
    279         // make this optional: psphotResid->save = true;
    280     }
    281 
    282     // Always do this, since we're using psphot for the background model.
    283     // XXX Not certain that I need to generate a separate pmFPAfile for PSPHOT.INPUT
    284     if (0) {
    285282        pmFPAfile *psphot = pmFPAfileDefineFromFPA(config, output->fpa, 1, 1, "PSPHOT.INPUT");
    286283        if (!psphot) {
    … …  
    297294            return false;
    298295        }
    299 
    300296    }
    301297

  • trunk/ppSub/src/ppSubDefineOutput.c

    r21360 r21524  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-06 01:37:17 $
     8 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-18 00:31:20 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
    1212
     13#ifdef HAVE_CONFIG_H
     14#include <config.h>
     15#endif
     16
     17#include <stdio.h>
     18#include <pslib.h>
     19#include <psmodules.h>
     20
    1321#include "ppSub.h"
    1422
    15 bool ppSubDefineOutput (pmConfig *config, const pmFPAview *view) {
     23bool ppSubDefineOutput(pmConfig *config, const pmFPAview *view)
     24{
     25    psAssert(config, "Require configuration");
     26    psAssert(view, "Require view");
    1627
    17     bool mdok;
    18 
    19     // generate an output readout
     28    // generate an output readout
    2029    pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT"); // Output cell
    2130    pmReadout *outRO = pmReadoutAlloc(outCell); // Output readout: subtraction
    … …  
    3039    pmReadout *refConv = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); // Reference readout
    3140
    32     // Add kernel descrption to header.  We don't know which readout gets the kernels because
    33     // it depends on which is larger (the choice is set in ppSubMatchPSFs)
    34     bool inputHasKernel = true;
    35     pmSubtractionKernels *kernels = psMetadataLookupPtr(&mdok, inConv->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); // Subtraction kernel
     41    // Add kernel descrption to header.
     42    // We don't know which readout has the kernels because it depends on which PSF is larger
     43    bool mdok;                          // Status of MD lookup
     44    psMetadata *analysis = inConv->analysis; // Analysis metadata with kernel information
     45    pmSubtractionKernels *kernels = psMetadataLookupPtr(&mdok, analysis,
     46                                                        PM_SUBTRACTION_ANALYSIS_KERNEL); // Subtraction kernel
    3647    if (!kernels) {
    37         kernels = psMetadataLookupPtr(&mdok, refConv->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL);
    38         inputHasKernel = false;
     48        analysis = refConv->analysis;
     49        kernels = psMetadataLookupPtr(&mdok, analysis, PM_SUBTRACTION_ANALYSIS_KERNEL);
    3950    }
    4051    if (!kernels) {
    … …  
    4556        return false;
    4657    }
    47     psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.KERNEL", 0, "Subtraction kernel", kernels->description);
     58    psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.KERNEL", 0, "Subtraction kernel",
     59                     kernels->description);
    4860
    49     if (inputHasKernel) {
    50         outRO->analysis = psMetadataCopy(outRO->analysis, inConv->analysis);
    51     } else {
    52         outRO->analysis = psMetadataCopy(outRO->analysis, refConv->analysis);
    53     }
     61    outRO->analysis = psMetadataCopy(outRO->analysis, analysis);
    5462
    5563#ifdef TESTING
    56     psImage *kernelImage = NULL;
    57 
    58     if (inputHasKernel) {
    59         kernelImage = psMetadataLookupPtr(&mdok, inConv->analysis, "SUBTRACTION.KERNEL.IMAGE"); // Image of the kernels
    60     } else {
    61         kernelImage = psMetadataLookupPtr(&mdok, refConv->analysis, "SUBTRACTION.KERNEL.IMAGE");
     64    {
     65        psImage *kernelImage = psMetadataLookupPtr(&mdok, analysis,
     66                                                   "SUBTRACTION.KERNEL.IMAGE"); // Image of kernel
     67        psFits *fits = psFitsOpen("kernel.fits", "w");
     68        psFitsWriteImage(fits, NULL, kernelImage, 0, NULL);
     69        psFitsClose(fits);
    6270    }
    63 
    64     psFits *fits = psFitsOpen("kernel.fits", "w");
    65     psFitsWriteImage(fits, NULL, kernelImage, 0, NULL);
    66     psFitsClose(fits);
    6771#endif
    6872

  • trunk/ppSub/src/ppSubLoop.c

    r21360 r21524  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-06 01:37:17 $
     8 *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-18 00:31:20 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
     12
     13#ifdef HAVE_CONFIG_H
     14#include <config.h>
     15#endif
     16
     17#include <pslib.h>
     18#include <psmodules.h>
     19#include <ppStats.h>
    1220
    1321#include "ppSub.h"
    … …  
    138146                    return false;
    139147                }
    140                 psImageMaskType maskValue = pmConfigMaskGet("MASK.VALUE", config);
     148                psImageMaskType maskValue = pmConfigMaskGet("MASK.VALUE", config);
    141149                ppStatsFPA(stats, output->fpa, view, maskValue, config);
    142150            }

  • trunk/ppSub/src/ppSubMakePSF.c

    r21360 r21524  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-06 01:37:17 $
     8 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-18 00:31:20 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
    1212
     13#ifdef HAVE_CONFIG_H
     14#include <config.h>
     15#endif
     16
     17#include <stdio.h>
     18#include <pslib.h>
     19#include <psmodules.h>
     20#include <psphot.h>
     21
    1322#include "ppSub.h"
    1423
    15 /**
    16  * Photometry stage 1: measure the PSF from the minuend image
    17  */
    18 bool ppSubMakePSF (pmConfig *config, const pmFPAview *view) {
     24bool ppSubMakePSF(pmConfig *config, const pmFPAview *view)
     25{
     26    psAssert(config, "Require configuration");
     27    psAssert(view, "Require view");
    1928
    20     bool mdok = false;
    21 
    22     // Photometry is to be performed in two stages:
    23     // 1. Measure the PSF using the PSF-matched images
    24     // 2. Find and measure sources on the subtracted image
    2529    psTimerStart("PPSUB_PHOT");
    2630
    27     // Look up recipe values
    2831    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSub
    2932    psAssert(recipe, "We checked this earlier, so it should be here.");
    3033
    31     psMetadata *psphotRecipe = psMetadataLookupMetadata(NULL, config->recipes, PSPHOT_RECIPE); // Recipe for psphot
     34    if (!psMetadataLookupBool(NULL, recipe, "PHOTOMETRY")) {
     35        return true;
     36    }
     37
     38    psMetadata *psphotRecipe = psMetadataLookupMetadata(NULL, config->recipes, PSPHOT_RECIPE);// psphot recipe
    3239    psAssert(recipe, "We checked this earlier, so it should be here.");
    33 
    34     if (!psMetadataLookupBool(NULL, recipe, "PHOTOMETRY")) return true;
    3540
    3641    bool reverse = psMetadataLookupBool(NULL, config->arguments, "REVERSE"); // Reverse sense of subtraction?
    3742
    38     pmReadout *minuend = NULL;
    39     pmFPAfile *minuendFile = NULL;
     43    bool mdok = false;                  // Status of MD lookup
     44    pmReadout *minuend = NULL;          // Image that will be positive following subtraction
     45    pmFPAfile *minuendFile = NULL;      // File for minuend image
    4046    if (reverse) {
    41         minuend = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV");
    42         minuendFile = psMetadataLookupPtr(&mdok, config->files, "PPSUB.REF.CONV");
     47        minuend = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV");
     48        minuendFile = psMetadataLookupPtr(&mdok, config->files, "PPSUB.REF.CONV");
    4349    } else {
    44         minuend = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV");
    45         minuendFile = psMetadataLookupPtr(&mdok, config->files, "PPSUB.INPUT.CONV");
     50        minuend = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV");
     51        minuendFile = psMetadataLookupPtr(&mdok, config->files, "PPSUB.INPUT.CONV");
    4652    }
    4753
    48     // supply the minuend pmFPAfile to psphot as PSPHOT.INPUT:
    49     psMetadataAddPtr (config->files, PS_LIST_TAIL, "PSPHOT.INPUT", PS_DATA_UNKNOWN | PS_META_REPLACE, "psphot input : view on another pmFPAfile", minuendFile);
     54#if 1
     55    pmReadout *template = minuend;
     56    pmFPAfile *photFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.INPUT"); // Photometry file
     57    pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout to photometer
     58    if (!photRO) {
     59        pmCell *cell = pmFPAviewThisCell(view, photFile->fpa); // Cell to photometer
     60        photRO = pmReadoutAlloc(cell); // Output readout: subtraction
     61    }
     62    photRO->image = psImageCopy(photRO->image, template->image, PS_TYPE_F32);
     63    if (template->variance) {
     64        photRO->variance = psImageCopy(photRO->variance, template->variance, PS_TYPE_F32);
     65    } else {
     66        psFree(photRO->variance);
     67        photRO->variance = NULL;
     68    }
     69    if (template->mask) {
     70        photRO->mask = psImageCopy(photRO->mask, template->mask, PS_TYPE_IMAGE_MASK);
     71    } else {
     72        psFree(photRO->mask);
     73        photRO->mask = NULL;
     74    }
     75#else
     76    // Supply the minuend pmFPAfile to psphot as PSPHOT.INPUT:
     77    psMetadataAddPtr(config->files, PS_LIST_TAIL, "PSPHOT.INPUT", PS_DATA_UNKNOWN | PS_META_REPLACE,
     78                     "psphot input: view on another pmFPAfile", minuendFile);
     79#endif
    5080
    51     // old-style variance renormalization
    52     if (!ppSubReadoutRenormPhot (config, recipe, minuend)) {
    53         psError(PS_ERR_UNKNOWN, false, "failure in renormalization");
    54         return false;
    55     }
    56 
    57     // extract the loaded sources from the associated readout
     81    // Extract the loaded sources from the associated readout, and generate PSF
     82    // Here, we assume the image is background-subtracted
    5883    pmReadout *sourcesRO = pmFPAfileThisReadout(config->files, view, "PPSUB.SOURCES");
    5984    psArray *sources = psMetadataLookupPtr(&mdok, sourcesRO->analysis, "PSPHOT.SOURCES");
    60 
    61     // generate PSF from the supplied sources (assumes image is background-subtracted)
    6285    if (!psphotReadoutFindPSF(config, view, sources)) {
    63         psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry on subtracted image.");
    64         return false;
     86        psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry on subtracted image.");
     87        return false;
    6588    }
    6689
    … …  
    7093    psMetadataItemSupplement(hdu->header, psphotRecipe, "FWHM_MAJ");
    7194    psMetadataItemSupplement(hdu->header, psphotRecipe, "FWHM_MIN");
     95
     96    // Get rid of the generated header; it will be regenerated by the real photometry run
     97    psMetadataRemoveKey(photRO->analysis, "PSPHOT.HEADER");
    7298
    7399    return true;
    … …  
    87113
    88114// Blow away the sources psphot found --- they're irrelevant for the subtraction
    89 // XXX is this still needed?  These are now probably not being set 
     115// XXX is this still needed?  These are now probably not being set
    90116// pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout with sources
    91117// psMetadataRemoveKey(photRO->analysis, "PSPHOT.SOURCES");

  • trunk/ppSub/src/ppSubMatchPSFs.c

    r21360 r21524  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-06 01:37:17 $
     8 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-18 00:31:20 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
    1212
     13#ifdef HAVE_CONFIG_H
     14#include <config.h>
     15#endif
     16
     17#include <stdio.h>
     18#include <pslib.h>
     19#include <psmodules.h>
     20#include <psphot.h>
     21
    1322#include "ppSub.h"
    1423
    15 /** Generate the PSF-matching kernel and convolve the images as needed.  Most of this function
    16  * involves looking up the parameters in the recipe and supplying them to the function
    17  * pmSubtractionMatch();
    18  */
    19 
    20 bool ppSubMatchPSFs (pmConfig *config, const pmFPAview *view) {
    21 
    22     bool mdok;                          // Status of MD lookup
     24bool ppSubMatchPSFs(pmConfig *config, const pmFPAview *view)
     25{
     26    psAssert(config, "Require configuration");
     27    psAssert(view, "Require view");
    2328
    2429    // Look up recipe values
    … …  
    2631    psAssert(recipe, "We checked this earlier, so it should be here.");
    2732
    28     // input images
     33    // Input images
    2934    pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); // Input readout
    3035    pmReadout *refRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF"); // Reference readout
    3136
    32     // get or generate output image holders
    33     pmReadout *inConv = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); // Input convolved readout
     37    // Output image holders
     38    pmReadout *inConv = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); // Input convolved
    3439    if (!inConv) {
    35         pmCell *cell = pmFPAfileThisCell(config->files, view, "PPSUB.INPUT.CONV"); // Input convolved readout
    36         inConv = pmReadoutAlloc(cell); // Convolved version of input
     40        pmCell *cell = pmFPAfileThisCell(config->files, view, "PPSUB.INPUT.CONV"); // Cell for convolved input
     41        inConv = pmReadoutAlloc(cell);
    3742    }
    38     pmReadout *refConv = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); // Reference convolved readout
     43    pmReadout *refConv = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); // Reference convolved
    3944    if (!refConv) {
    40         pmCell *cell = pmFPAfileThisCell(config->files, view, "PPSUB.REF.CONV"); // Input convolved readout
    41         refConv = pmReadoutAlloc(cell); // Convolved version of input
     45        pmCell *cell = pmFPAfileThisCell(config->files, view, "PPSUB.REF.CONV"); // Cell for convolved ref.
     46        refConv = pmReadoutAlloc(cell);
    4247    }
    4348
    44     // options which control the generation of optimal parameters
    45     bool optimum = psMetadataLookupBool(&mdok, recipe, "OPTIMUM"); // Derive optimum parameters?
    46     float optMin = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.MIN"); // Minimum width for search
    47     float optMax = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.MAX"); // Maximum width for search
    48     float optStep = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.STEP"); // Step for search
     49    bool mdok;                          // Status of MD lookup
    4950
    50     // Vector with FWHMs for optimum search
    51     psVector *optWidths = NULL;
    52     if (optimum) {
    53         optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32);
    54     }
    55 
    56     // Sources in image : these must be loaded from previous analysis stages
    57     pmReadout *sourcesRO = pmFPAfileThisReadout(config->files, view, "PPSUB.SOURCES"); // Readout with sources
    58     psArray *sources = NULL;           
     51    // Sources in image, used for stamps: these must be loaded from previous analysis stages
     52    pmReadout *sourcesRO = pmFPAfileThisReadout(config->files, view, "PPSUB.SOURCES"); // Readout with sources
     53    psArray *sources = NULL;            // Sources in image; used for stamps
    5954    if (sourcesRO) {
    6055        sources = psMetadataLookupPtr(&mdok, sourcesRO->analysis, "PSPHOT.SOURCES");
    … …  
    6661    float spacing = psMetadataLookupF32(NULL, recipe, "STAMP.SPACING"); // Typical stamp spacing
    6762    float threshold = psMetadataLookupF32(NULL, recipe, "STAMP.THRESHOLD"); // Threshold for stmps
    68 
    6963    const char *stampsName = psMetadataLookupStr(&mdok, config->arguments, "STAMPS"); // Filename for stamps
    7064
    7165    const char *typeStr = psMetadataLookupStr(NULL, recipe, "KERNEL.TYPE"); // Kernel type
    7266    psAssert(typeStr, "We put it here in ppSubArguments.c");
    73 
    7467    pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Type of kernel
    7568    if (type == PM_SUBTRACTION_KERNEL_NONE) {
    … …  
    8780    int binning = psMetadataLookupS32(NULL, recipe, "SPAM.BINNING"); // Binning for SPAM kernel
    8881    float penalty = psMetadataLookupF32(NULL, recipe, "PENALTY"); // Penalty for wideness
    89 
    90     int optOrder = psMetadataLookupS32(&mdok, recipe, "OPTIMUM.ORDER"); // Order for search
    91     float optThresh = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.TOL"); // Tolerance for search
    92 
    9382    int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations
    9483    float rej = psMetadataLookupF32(NULL, recipe, "REJ"); // Rejection threshold
    … …  
    9887    float poorFrac = psMetadataLookupF32(&mdok, recipe, "POOR.FRACTION"); // Fraction for "poor"
    9988
    100     // XXX EAM : do we need to / want to define different values for BAD and POOR subtraction vs BAD and POOR warp?
    101     psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config); // Bits to mask going in to pmSubtractionMatch
     89    // Options which control the generation of optimal parameters
     90    bool optimum = psMetadataLookupBool(&mdok, recipe, "OPTIMUM"); // Derive optimum parameters?
     91    float optMin = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.MIN"); // Minimum width for search
     92    float optMax = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.MAX"); // Maximum width for search
     93    float optStep = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.STEP"); // Step for search
     94    psVector *optWidths = NULL;         // Vector with FWHMs for optimum search
     95    if (optimum) {
     96        optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32);
     97    }
     98    int optOrder = psMetadataLookupS32(&mdok, recipe, "OPTIMUM.ORDER"); // Order for search
     99    float optThresh = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.TOL"); // Tolerance for search
     100
     101    // XXX Do we need/want to define different values for BAD and POOR subtraction vs BAD and POOR warp?
     102    psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config); // Bits to mask in inputs
    102103    psImageMaskType maskPoor = pmConfigMaskGet("POOR.WARP", config); // Bits to mask for poor pixels
    103104    psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask for bad pixels
    … …  
    124125    pmSubtractionThreadsFinalize(inRO, refRO);
    125126
     127    psImageCovarianceTransfer(inConv->variance, inConv->covariance);
     128    psImageCovarianceTransfer(refConv->variance, refConv->covariance);
     129
    126130    // XXX drop the pixels associated with inRO and refRO (now that we have inConv and refConf)
    127131

  • trunk/ppSub/src/ppSubReadout.c

    r21424 r21524  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.112 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-09 21:26:05 $
     8 *  @version $Revision: 1.113 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-18 00:31:20 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
     12
     13#ifdef HAVE_CONFIG_H
     14#include <config.h>
     15#endif
     16
     17#include <stdio.h>
     18#include <pslib.h>
     19#include <psmodules.h>
    1220
    1321#include "ppSub.h"
    … …  
    1725    psTimerStart("PPSUB_MATCH");
    1826
    19     if (!ppSubSetMasks (config, view)) {
    20         psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
     27    if (!ppSubSetMasks(config, view)) {
     28        psError(PS_ERR_UNKNOWN, false, "Unable to set masks.");
    2129        return false;
    2230    }
    2331
    24     if (!ppSubMatchPSFs (config, view)) {
    25         psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
     32    if (!ppSubMatchPSFs(config, view)) {
     33        psError(PS_ERR_UNKNOWN, false, "Unable to match PSFs.");
    2634        return false;
    2735    }
    2836
    29     if (!ppSubDefineOutput (config, view)) {
    30         psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
    31         return false;
    32     }
    33 
    34     if (!ppSubVarianceFactors (config, stats, view)) {
    35         psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
     37    if (!ppSubDefineOutput(config, view)) {
     38        psError(PS_ERR_UNKNOWN, false, "Unable to define output.");
    3639        return false;
    3740    }
    3841
    3942    if (!ppSubMakePSF(config, view)) {
    40         psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
     43        psError(PS_ERR_UNKNOWN, false, "Unable to generate PSF.");
    4144        return false;
    4245    }
    4346
    44     if (!ppSubReadoutSubtract (config, view)) {
    45         psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
     47    if (!ppSubReadoutSubtract(config, view)) {
     48        psError(PS_ERR_UNKNOWN, false, "Unable to subtract images.");
    4649        return false;
    4750    }
    … …  
    5356    }
    5457
    55     if (!ppSubReadoutPhotometry (config, stats, view)) {
    56         psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
     58    if (!ppSubReadoutPhotometry(config, stats, view)) {
     59        psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry.");
    5760        return false;
    5861    }
    5962
    60     if (!ppSubReadoutUpdate (config, view)) {
    61         psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
     63    if (!ppSubReadoutUpdate(config, stats, view)) {
     64        psError(PS_ERR_UNKNOWN, false, "Unable to update.");
    6265        return false;
    6366    }

  • trunk/ppSub/src/ppSubReadoutPhotometry.c

    r21374 r21524  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-06 03:20:02 $
     8 *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-18 00:31:20 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
    1212
     13#ifdef HAVE_CONFIG_H
     14#include <config.h>
     15#endif
     16
     17#include <stdio.h>
     18#include <pslib.h>
     19#include <psmodules.h>
     20#include <psphot.h>
     21
    1322#include "ppSub.h"
    1423
    15 /**
    16  * Photometry stage 2: find and measure sources on the subtracted image
    17  */
    18 bool ppSubReadoutPhotometry (pmConfig *config, psMetadata *stats, const pmFPAview *view) {
    19 
    20     bool mdok = false;
     24bool ppSubReadoutPhotometry (pmConfig *config, psMetadata *stats, const pmFPAview *view)
     25{
     26    psAssert(config, "Require configuration");
     27    psAssert(view, "Require view");
    2128
    2229    // Look up recipe values
    … …  
    3239    // The PSF (measured in ppSubMakePSF) is stored on the chip->analysis of PSPHOT.INPUT
    3340    // In order to use an incoming PSF, it must be stored on the chip->analysis of PSPHOT.PSF.LOAD
    34     pmChip *psfInputChip = pmFPAfileThisChip(config->files, view, "PSPHOT.INPUT");
     41    pmChip *psfInputChip = pmFPAfileThisChip(config->files, view, "PSPHOT.INPUT"); // Chip with PSF
    3542    psAssert (psfInputChip, "should have been generated for ppSubMakePSF");
    36 
    37     pmChip *psfLoadChip = pmFPAfileThisChip(config->files, view, "PSPHOT.PSF.LOAD");
     43    pmChip *psfLoadChip = pmFPAfileThisChip(config->files, view, "PSPHOT.PSF.LOAD"); // Chip to have PSF
    3844    psAssert (psfLoadChip, "PSPHOT.PSF.LOAD should have been defined in ppSubCamera");
    39 
    40     pmPSF *psf = psMetadataLookupPtr(NULL, psfInputChip->analysis, "PSPHOT.PSF");
     45    pmPSF *psf = psMetadataLookupPtr(NULL, psfInputChip->analysis, "PSPHOT.PSF"); // PSF for photometry
    4146    if (!psf) {
    4247        psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find PSF from psphot");
    4348        return false;
    4449    }
    45     psMetadataAddPtr(psfLoadChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN | PS_META_REPLACE, "PSF from matched addition", psf);
     50    psMetadataAddPtr(psfLoadChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN | PS_META_REPLACE,
     51                     "PSF from ppSubMakePSF", psf);
     52
     53    bool mdok = false;
    4654
    4755    // psphotReadoutMinimal performs the photometry analysis on PSPHOT.INPUT; we need to move
    4856    // around the pointers so PSPHOT.INPUT corresponds to the output image; previously, it was
    4957    // equivalent to the minuend image.
    50     pmFPAfile *outputFile = psMetadataLookupPtr(&mdok, config->files, "PPSUB.OUTPUT");
     58    pmFPAfile *outputFile = psMetadataLookupPtr(&mdok, config->files, "PPSUB.OUTPUT"); // Output file
    5159    pmReadout *outRO = pmFPAviewThisReadout(view, outputFile->fpa); // Readout with the sources
    5260
    5361    // XXX possibly rename this to PPSUB.RESID?
    54     pmFPAfile *photFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.RESID");
    55     pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout with the sources
     62    pmFPAfile *photFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.INPUT"); // Photometry file
     63    pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout to photometer
    5664    if (!photRO) {
    57         pmCell *cell = pmFPAfileThisCell(config->files, view, "PSPHOT.RESID"); // Output cell
     65        pmCell *cell = pmFPAviewThisCell(view, photFile->fpa); // Cell to photometer
    5866        photRO = pmReadoutAlloc(cell); // Output readout: subtraction
    5967    }
    … …  
    6270        photRO->variance = psImageCopy(photRO->variance, outRO->variance, PS_TYPE_F32);
    6371    } else {
    64         psFree (photRO->variance);
     72        psFree(photRO->variance);
    6573        photRO->variance = NULL;
    6674    }
    … …  
    6876        photRO->mask = psImageCopy(photRO->mask, outRO->mask, PS_TYPE_IMAGE_MASK);
    6977    } else {
    70         psFree (photRO->mask);
     78        psFree(photRO->mask);
    7179        photRO->mask = NULL;
    7280    }
    7381
    74     // pmFPAfile *photFile = outputFile;
    75     psMetadataAddPtr (config->files, PS_LIST_TAIL, "PSPHOT.INPUT", PS_DATA_UNKNOWN | PS_META_REPLACE, "psphot input : view on another pmFPAfile", photFile);
    76 
    77     // old-style variance renormalization
    78     if (!ppSubReadoutRenormPhot (config, recipe, photRO)) {
    79         psError(PS_ERR_UNKNOWN, false, "failure in renormalization");
    80         return false;
    81     }
     82#if 0
     83    psMetadataAddPtr(config->files, PS_LIST_TAIL, "PSPHOT.INPUT", PS_DATA_UNKNOWN | PS_META_REPLACE,
     84                     "psphot input: view on another pmFPAfile", photFile);
     85#endif
    8286
    8387    if (!psphotReadoutMinimal(config, view)) {
    … …  
    8690        psErrorClear();
    8791    }
     92#if 1
     93    photRO->data_exists = true;
     94    photRO->parent->data_exists = true;
     95    photRO->parent->parent->data_exists = true;
     96#endif
    8897
    8998    if (stats) {
    … …  
    101110}
    102111
    103     // XXX not sure that this is still needed (only if psphotReadoutMinimal measures the background)
    104     // if (!pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL") ||
    105     //  !pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV") ||
    106     //  !pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND")) {
    107     //  psError(PS_ERR_UNKNOWN, false, "Unable to drop PSPHOT internal files.");
    108     //  return false;
    109     // }
    110112
    111     // pmFPAfileActivate(config->files, false, "PSPHOT.INPUT");
    112     // pmFPAfileActivate(config->files, false, "PSPHOT.LOAD.PSF");
     113
    113114
    114115

  • trunk/ppSub/src/ppSubReadoutSubtract.c

    r21424 r21524  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-09 21:26:05 $
     8 *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-18 00:31:20 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
    1212
     13#ifdef HAVE_CONFIG_H
     14#include <config.h>
     15#endif
     16
     17#include <stdio.h>
     18#include <pslib.h>
     19#include <psmodules.h>
     20
    1321#include "ppSub.h"
     22
    1423#define WCS_TOLERANCE 0.001             // Tolerance for WCS
    1524
    16 bool ppSubReadoutSubtract (pmConfig *config, const pmFPAview *view) {
     25bool ppSubReadoutSubtract(pmConfig *config, const pmFPAview *view)
     26{
     27    psAssert(config, "Require configuration");
     28    psAssert(view, "Require view");
    1729
    18     bool mdok = false;
    1930
    2031    // Look up recipe values
    21     psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSUB_RECIPE); // Recipe for ppSim
     32    bool mdok = false;                  // Status of MD lookup
     33    psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSUB_RECIPE); // Recipe for ppSub
    2234    psAssert(recipe, "We checked this earlier, so it should be here.");
    2335
    … …  
    2537
    2638    // Subtraction is: minuend - subtrahend
    27     pmReadout *minuend = NULL;
    28     pmReadout *subtrahend = NULL;
     39    pmReadout *minuend = NULL;          // Positive image
     40    pmReadout *subtrahend = NULL;       // Negative image
    2941    if (reverse) {
    3042        minuend = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV");
    … …  
    5971    outRO->mask = (psImage*)psBinaryOp(outRO->mask, minuend->mask, "|", subtrahend->mask);
    6072
     73    psArray *covars = psArrayAlloc(2);  // Covariance pseudo-matrices
     74    covars->data[0] = psMemIncrRefCounter(minuend->covariance);
     75    covars->data[1] = psMemIncrRefCounter(subtrahend->covariance);
     76    outRO->covariance = psImageCovarianceSum(covars);
     77    psFree(covars);
     78
    6179    outRO->data_exists = true;
    6280    outRO->parent->data_exists = true;
    … …  
    6583    pmSubtractionVisualShowSubtraction(minuend->image, subtrahend->image, outRO->image);
    6684
    67     // copy concepts from the input to the output (XXX should this always use minuend?)
    68     pmFPAfile *inFile = psMetadataLookupPtr (&mdok, config->files, "PPSUB.INPUT");
    69     pmFPA *inFPA = inFile->fpa;
    70 
    71     pmFPAfile *outFile = psMetadataLookupPtr (&mdok, config->files, "PPSUB.OUTPUT");
    72     pmFPA *outFPA = outFile->fpa;
    73 
     85    // Copy concepts from the input to the output
     86    pmFPAfile *inFile = psMetadataLookupPtr(&mdok, config->files, "PPSUB.INPUT"); // Input file
     87    pmFPA *inFPA = inFile->fpa;         // Input FPA
     88    pmFPAfile *outFile = psMetadataLookupPtr(&mdok, config->files, "PPSUB.OUTPUT"); // Output file
     89    pmFPA *outFPA = outFile->fpa;       // Output FPA
    7490    if (!pmFPACopyConcepts(outFPA, inFPA)) {
    7591        psError(PS_ERR_UNKNOWN, false, "Unable to copy concepts from input to output.");
    … …  
    7894    }
    7995
    80     // get the HDUs and output Chip to copy the astrometry
    81     pmHDU *inHDU = inFPA->hdu;        // input HDU
    82     pmHDU *outHDU = outFPA->hdu;
    83     pmChip *outChip = pmFPAfileThisChip(config->files, view, "PPSUB.OUTPUT");
    84 
     96    // Copy astrometry over
     97    // It should find its way into the output images and photometry
     98    pmHDU *inHDU = inFPA->hdu;          // Input HDU
     99    pmHDU *outHDU = outFPA->hdu;        // Output HDU
     100    pmChip *outChip = pmFPAfileThisChip(config->files, view, "PPSUB.OUTPUT"); // Output chip
    85101    if (!outHDU || !inHDU) {
    86102        psError(PS_ERR_UNKNOWN, false, "Unable to find HDU at FPA level to copy astrometry.");
    87103        return false;
    88104    }
    89 
    90     // Copy astrometry over
    91     // It should get into the output images and photometry
    92105    if (!pmAstromReadWCS(outFPA, outChip, inHDU->header, 1.0)) {
    93106        psError(PS_ERR_UNKNOWN, false, "Unable to read WCS astrometry from input FPA.");
    94107        return false;
    95108    }
    96 
    97109    if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_TOLERANCE)) {
    98110        psError(PS_ERR_UNKNOWN, false, "Unable to write WCS astrometry to output FPA.");
    … …  
    102114    return true;
    103115}
    104 
    105 // XXX this test code was in place to check for and squash unexpected NANs
    106 
    107 #if 0
    108     pmReadoutMaskApply(outRO, maskBad);
    109 
    110     for (int y = 0; y < outRO->image->numRows; y++) {
    111         for (int x = 0; x < outRO->image->numCols; x++) {
    112             if (isnan(outRO->image->data.F32[y][x]) && !(outRO->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) {
    113                 printf("Unmasked NAN at %d %d --> %d\n", x, y, outRO->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x]);
    114             }
    115         }
    116     }
    117 #endif
    118 

  • trunk/ppSub/src/ppSubReadoutUpdate.c

    r21396 r21524  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-07 00:15:00 $
     8 *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-18 00:31:20 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
    1212
     13#ifdef HAVE_CONFIG_H
     14#include <config.h>
     15#endif
     16
     17#include <stdio.h>
     18#include <pslib.h>
     19#include <psmodules.h>
     20
    1321#include "ppSub.h"
    1422
    15 /**
    16  * Renormalize, update headers and generate JPEGs
    17  */
    18 bool ppSubReadoutUpdate (pmConfig *config, const pmFPAview *view) {
     23bool ppSubReadoutUpdate(pmConfig *config, psMetadata *stats, const pmFPAview *view)
     24{
     25    psAssert(config, "Require configuration");
     26    psAssert(view, "Require view");
    1927
    20     bool mdok = false;
     28    bool mdok = false;                  // Status of MD lookup
    2129
    2230    // Look up recipe values
    … …  
    2432    psAssert(recipe, "We checked this earlier, so it should be here.");
    2533
    26     // select the output readout
    27     pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT");
    28 
    29     // Renormalising for pixels, because that's what magic desires
    30     if (!ppSubReadoutRenormPixels (config, recipe, outRO)) {
    31         psError(PS_ERR_BAD_PARAMETER_VALUE, false, "failure in renormalization");
    32         return false;
    33     }
    34 
    35     // select the output FPA and HDU to get the output header
    36     pmFPAfile *outFile = psMetadataLookupPtr (&mdok, config->files, "PPSUB.OUTPUT");
    37     pmFPA *outFPA = outFile->fpa;
    38     pmHDU *outHDU = outFPA->hdu;
     34    pmFPAfile *outFile = psMetadataLookupPtr(&mdok, config->files, "PPSUB.OUTPUT"); // Output file
     35    pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT"); // Output image
     36    pmFPA *outFPA = outFile->fpa;       // Output FPA
     37    pmHDU *outHDU = outFPA->hdu;        // Output HDU
    3938
    4039    // Add additional data to the header
    … …  
    4645                     "Subtraction input", inFile->filename);
    4746
     47    // Statistics on the matching
     48    if (stats) {
     49        psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MODE);
     50        psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_STAMPS);
     51        psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_DEV_MEAN);
     52        psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_DEV_RMS);
     53        psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_NORM);
     54        psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_BGDIFF);
     55        psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MX);
     56        psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MY);
     57        psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MXX);
     58        psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MXY);
     59        psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MYY);
     60
     61        psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_MATCH", 0, "Time to match PSFs",
     62                         psTimerClear("PPSUB_MATCH"));
     63    }
     64
    4865    // Generate binned JPEGs
    4966    {
    … …  
    5471
    5572        // Target cells
    56         pmCell *cell1 = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT.JPEG1");
    57         pmCell *cell2 = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT.JPEG2");
     73        pmCell *cell1 = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT.JPEG1"); // Rebinned cell once
     74        pmCell *cell2 = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT.JPEG2"); // Rebinned cell twice
    5875
    59         pmReadout *ro1 = pmReadoutAlloc(cell1);
    60         pmReadout *ro2 = pmReadoutAlloc(cell2); // Binned readouts
     76        pmReadout *ro1 = pmReadoutAlloc(cell1), *ro2 = pmReadoutAlloc(cell2); // Binned readouts
    6177        if (!pmReadoutRebin(ro1, outRO, maskBad, bin1, bin1)) {
    6278            psError(PS_ERR_UNKNOWN, false, "Unable to bin output (1st binning)");

  • trunk/ppSub/src/ppSubSetMasks.c

    r21360 r21524  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-06 01:37:17 $
     8 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-18 00:31:20 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
    1212
     13#ifdef HAVE_CONFIG_H
     14#include <config.h>
     15#endif
     16
     17#include <stdio.h>
     18#include <pslib.h>
     19#include <psmodules.h>
     20#include <psphot.h>
     21
    1322#include "ppSub.h"
    1423
    15 /** this function generates (if needed) and sets or updates the masks for both the input and
    16  * reference images.  this function also has the code for interpolation over bad pixels, but it
    17  * is currently if-def-ed out
    18  */
     24bool ppSubSetMasks(pmConfig *config, const pmFPAview *view)
     25{
     26    psAssert(config, "Require configuration");
     27    psAssert(view, "Require view");
    1928
    20 bool ppSubSetMasks (pmConfig *config, const pmFPAview *view) {
    21 
    22     psImageMaskType maskValue;
    23     psImageMaskType markValue;
    24     bool mdok = false;
    25 
     29    psImageMaskType maskValue, markValue; // Mask values
    2630    if (!pmConfigMaskSetBits(&maskValue, &markValue, config)) {
    2731        psError(PS_ERR_UNKNOWN, false, "Unable to determine mask value.");
    … …  
    2933    }
    3034
    31     // set the mask bits needed by psphot (in psphot recipe)
     35    // Set the mask bits needed by psphot (in psphot recipe)
    3236    psphotSetMaskRecipe (config, maskValue, markValue);
    3337
    3438    // Look up recipe values
    35     psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSUB_RECIPE); // Recipe for ppSim
     39    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim
    3640    psAssert(recipe, "We checked this earlier, so it should be here.");
    3741
    3842    psImageMaskType satValue = pmConfigMaskGet("SAT", config);
    39     psAssert (satValue, "SAT must be non-zero");
     43    psAssert(satValue, "SAT must be non-zero");
    4044
    4145    psImageMaskType badValue = pmConfigMaskGet("BAD", config);
    42     psAssert (badValue, "BAD must be non-zero");
     46    psAssert(badValue, "BAD must be non-zero");
    4347
    4448    // input images
    … …  
    4953    psImage *reference = refRO->image;  // Reference image
    5054    PS_ASSERT_IMAGES_SIZE_EQUAL(input, reference, false);
    51 
    52     // XXX use psImageCopy below to avoid caring about image dimensions
    53     int numCols = input->numCols;
    54     int numRows = input->numRows;
     55    int numCols = input->numCols, numRows = input->numRows; // Size of image
    5556
    5657    // Generate masks if they don't exist
    … …  
    7879    }
    7980    if (!pmReadoutMaskNonfinite(refRO, satValue)) {
    80         psError(PS_ERR_UNKNOWN, false, "Unable to mask non-finite pixels in reference.");
    81         return false;
     81        psError(PS_ERR_UNKNOWN, false, "Unable to mask non-finite pixels in reference.");
     82        return false;
    8283    }
    8384
    84 #if (0)
    85     // Look up recipe values
    86     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim
    87     psAssert(recipe, "We checked this earlier, so it should be here.");
    8885
     86#if 0
     87    // Interpolation over bad pixels: this takes a while!
    8988    bool mdok = false;
    90 
    9189    psString interpModeStr = psMetadataLookupStr(&mdok, recipe, "INTERPOLATION"); // Interpolation mode
    9290    psImageInterpolateMode interpMode = psImageInterpolateModeFromString(interpModeStr); // Interp
    … …  
    102100    // Interpolate over bad pixels, so the bad pixels don't explode
    103101    if (!pmReadoutInterpolateBadPixels(inRO, maskVal, interpMode, poorFrac, maskPoor, maskBad)) {
    104         psError(PS_ERR_UNKNOWN, false, "Unable to interpolate bad pixels for input image.");
    105         return false;
     102        psError(PS_ERR_UNKNOWN, false, "Unable to interpolate bad pixels for input image.");
     103        return false;
    106104    }
    107105    if (!pmReadoutInterpolateBadPixels(refRO, maskVal, interpMode, poorFrac, maskPoor, maskBad)) {
    108         psError(PS_ERR_UNKNOWN, false, "Unable to interpolate bad pixels for reference image.");
    109         return false;
     106        psError(PS_ERR_UNKNOWN, false, "Unable to interpolate bad pixels for reference image.");
     107        return false;
    110108    }
    111109    maskVal |= maskBad;

  • trunk/ppSub/src/ppSubVarianceFactors.c

    r21374 r21524  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-06 03:20:02 $
     8 *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-18 00:31:20 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
    … …  
    9696    vfItem->data.F32 = vf;
    9797
    98     // Statistics on the matching
    99     if (stats) {
    100         psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MODE);
    101         psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_STAMPS);
    102         psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_DEV_MEAN);
    103         psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_DEV_RMS);
    104         psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_NORM);
    105         psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_BGDIFF);
    106         psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MX);
    107         psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MY);
    108         psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MXX);
    109         psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MXY);
    110         psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MYY);
    111 
    112         psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_MATCH", 0, "Time to match PSFs",
    113                          psTimerClear("PPSUB_MATCH"));
    114     }
    115 
    11698    return true;
    11799}

  • trunk/ppSub/src/ppSubVersion.c

    r21360 r21524  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-06 01:37:17 $
     8 *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-18 00:31:20 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
     12
     13#ifdef HAVE_CONFIG_H
     14#include <config.h>
     15#endif
     16
     17#include <stdio.h>
     18#include <pslib.h>
     19#include <psmodules.h>
     20#include <ppStats.h>
    1221
    1322#include "ppSub.h"
Note: See TracChangeset for help on using the changeset viewer.