IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 21257


Ignore:
Timestamp:
Feb 1, 2009, 11:55:34 AM (17 years ago)
Author:
eugene
Message:

re-factored ppSubReadout into a sequence of logical steps; defined convolved images as their own pmFPAfiles; switched to new psphot functions (psphotReadoutFindPSF and psphotReadoutMinimal); simplified mask and psphot visualization handling

Location:
trunk/ppSub/src
Files:
10 added
9 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppSub/src/Makefile.am

    r20611 r21257  
    44ppSub_LDFLAGS  = $(PSLIB_LIBS)   $(PSMODULE_LIBS)   $(PPSTATS_LIBS)   $(PSPHOT_LIBS)   $(PPSUB_LIBS)
    55
    6 ppSub_SOURCES =                 \
    7         ppSub.c                 \
    8         ppSubArguments.c        \
    9         ppSubBackground.c       \
    10         ppSubCamera.c           \
    11         ppSubLoop.c             \
    12         ppSubReadout.c          \
    13         ppSubVersion.c
     6ppSub_SOURCES =                  \
     7        ppSub.c                  \
     8        ppSubArguments.c         \
     9        ppSubVersion.c           \
     10        ppSubBackground.c        \
     11        ppSubCamera.c            \
     12        ppSubLoop.c              \
     13        ppSubReadout.c           \
     14        ppSubDefineOutput.c      \
     15        ppSubExtras.c            \
     16        ppSubMakePSF.c           \
     17        ppSubMatchPSFs.c         \
     18        ppSubReadoutPhotometry.c \
     19        ppSubReadoutSubtract.c   \
     20        ppSubReadoutUpdate.c     \
     21        ppSubSetMasks.c          \
     22        ppSubReadoutRenorm.c     \
     23        ppSubVarianceFactors.c
    1424
    1525ppSubKernel_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PPSUB_CFLAGS)

  • trunk/ppSub/src/ppSub.c

    r20413 r21257  
    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>
    9 
    101#include "ppSub.h"
    112

  • trunk/ppSub/src/ppSub.h

    r20611 r21257  
    11#ifndef PP_SUB_H
    22#define PP_SUB_H
     3
     4#ifdef HAVE_CONFIG_H
     5#include <config.h>
     6#endif
     7
     8#include <stdio.h>
     9#include <string.h>
     10#include <pslib.h>
     11#include <psmodules.h>
     12#include <psphot.h>
     13#include <ppStats.h>
    314
    415#define PPSUB_RECIPE "PPSUB"            /// Name of the recipe to use
    … …  
    617/// Setup the arguments parsing
    718bool ppSubArgumentsSetup(int argc, char *argv[], ///< Command-line arguments
    8                     pmConfig *config    ///< Configuration
     19                        pmConfig *config    ///< Configuration
    920    );
    1021
    … …  
    2738    );
    2839
     40/// generate (if needed) and set or update the masks for input and reference images
     41bool ppSubSetMasks (pmConfig *config,     ///< Configuration
     42                    const pmFPAview *view ///< View of active readout
     43    );
     44
     45/// Generate the PSF-matching kernel and convolve the images as needed.  Most of this function
     46/// involves looking up the parameters in the recipe and supplying them to the function
     47/// pmSubtractionMatch()
     48bool ppSubMatchPSFs (pmConfig *config,    ///< Configuration
     49                     const pmFPAview *view ///< View of active readout
     50    );
     51
     52/// generate the output readout and pass the kernel info to the header
     53bool ppSubDefineOutput(pmConfig *config,          ///< Configuration
     54                       const pmFPAview *view ///< View of active readout
     55    );
     56
     57/// Calculate the variance factor for the output image based on the input images
     58bool ppSubVarianceFactors(pmConfig *config,       ///< Configuration
     59                          psMetadata *stats,    ///< Statistics, for output
     60                          const pmFPAview *view ///< View of active readout
     61    );
     62
     63/// Photometry stage 1: measure the PSF from the minuend image
     64bool ppSubMakePSF(pmConfig *config,       ///< Configuration
     65                  const pmFPAview *view ///< View of active readout
     66    );
     67
     68/// Perform the actual image subtraction, update output concepts
     69bool ppSubReadoutSubtract(pmConfig *config,       ///< Configuration
     70                          const pmFPAview *view ///< View of active readout
     71    );
     72
     73
     74/// Photometry stage 2: find and measure sources on the subtracted image
     75bool ppSubReadoutPhotometry(pmConfig *config,     ///< Configuration
     76                            psMetadata *stats,    ///< Statistics, for output
     77                            const pmFPAview *view ///< View of active readout
     78    );
     79
     80/// Renormalize, update headers and generate JPEGs
     81bool ppSubReadoutUpdate(pmConfig *config,         ///< Configuration
     82                        const pmFPAview *view ///< View of active readout
     83    );
     84
    2985/// Higher-order background subtraction
    30 bool ppSubBackground(pmReadout *ro,     ///< Readout of interest
    31                      pmConfig *config,  ///< Configuration
     86bool ppSubBackground(pmConfig *config,  ///< Configuration
    3287                     const pmFPAview *view ///< View to readout
    3388    );
    … …  
    3792    );
    3893
     94/// Generate and Set the masks if needed
     95bool ppSubSetMasks (
     96    pmConfig *config,                   ///< Configuration
     97    const pmFPAview *view               ///< view to readout
     98    );
     99
     100/// Renormalize readout for peak pixels
     101bool ppSubReadoutRenormPixels (
     102    pmConfig *config,                   ///< Configuration
     103    psMetadata *recipe,                         ///< Recipe
     104    pmReadout *readout                  ///< Readout
     105    );
     106
     107/// Renormalize readout for photometry analysis
     108bool ppSubReadoutRenormPhot (
     109    pmConfig *config,                   ///< Configuration
     110    psMetadata *recipe,                         ///< Recipe
     111    pmReadout *readout                  ///< Readout
     112    );
     113
     114// Copy every instance of a single keyword from one metadata to another
     115bool psMetadataCopySingle(psMetadata *target, psMetadata *source, const char *name);
    39116
    40117#endif

  • trunk/ppSub/src/ppSubArguments.c

    r21183 r21257  
    1 #ifdef HAVE_CONFIG_H
    2 #include <config.h>
    3 #endif
    4 
    5 #include <stdio.h>
    6 #include <string.h>
    7 #include <pslib.h>
    8 #include <psmodules.h>
    9 #include <psphot.h>
    10 
    111#include "ppSub.h"
    122
    … …  
    172162bool ppSubArgumentsSetup(int argc, char *argv[], pmConfig *config)
    173163{
     164    int argnum;                         // Argument Number
    174165    assert(config);
    175166
    … …  
    180171    }
    181172
    182     {
    183         int arg;                        // Argument Number
    184         if ((arg = psArgumentGet(argc, argv, "-psphot-visual"))) {
    185             psArgumentRemove(arg, &argc, argv);
    186             psMetadataAddBool(recipe, PS_LIST_TAIL, "PSPHOT.VISUAL", 0, "Visual guide to psphot?", true);
    187         }
     173    if ((argnum = psArgumentGet(argc, argv, "-psphot-visual"))) {
     174        psArgumentRemove(argnum, &argc, argv);
     175        psphotSetVisual(true);
    188176    }
    189177
    … …  
    381369    }
    382370
     371    // XXX move this to ppSubArguments
     372    int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads
     373    if (threads > 0) {
     374        if (!psThreadPoolInit(threads)) {
     375            psError(PS_ERR_UNKNOWN, false, "Unable to setup %d threads", threads);
     376            return false;
     377        }
     378    }
     379
    383380    return true;
    384381

  • trunk/ppSub/src/ppSubBackground.c

    r21183 r21257  
    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>
    91#include "ppSub.h"
    102
    11 // Copied from ppImageSubtractBackground()
    12 bool ppSubBackground(pmReadout *ro, pmConfig *config, const pmFPAview *view)
     3// Based on ppImageSubtractBackground()
     4bool ppSubBackground(pmConfig *config, const pmFPAview *view)
    135{
    146    psAssert(config, "Need configuration");
    157    psAssert(view, "Need view to chip");
    168
    17     // The background model file may be defined, though the model hasn't been built
    18     // If so, we will build it below.
    19     bool status;                        // Status of lookup
    20     pmFPAfile *modelFile = psMetadataLookupPtr(&status, config->files, "PSPHOT.BACKMDL"); // Background model
    21     if (!status) {
    22         psError(PS_ERR_UNEXPECTED_NULL, true, "PSPHOT.BACKMDL file is not defined");
    23         return false;
    24     }
     9    bool status; // Status of metadata lookups
    2510
    2611    psMetadata *ppSubRecipe = psMetadataLookupPtr(NULL, config->recipes, PPSUB_RECIPE);
    2712    psAssert(ppSubRecipe, "Need PPSUB recipe");
     13
    2814    psMetadata *psphotRecipe = psMetadataLookupPtr(NULL, config->recipes, PSPHOT_RECIPE);
    29     psAssert(psphotRecipe, "Need PSPHOT recipe");
     15    psAssert(psphotRecipe, "Need PSPHOT recipe for binning");
    3016
    31     psString maskBadStr = psMetadataLookupStr(NULL, ppSubRecipe, "MASK.BAD"); // Name of bits to mask for bad
    32     psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
     17    psImageMaskType maskBad = pmConfigMaskGet("BLANK", config);
    3318
    34     // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    35     psMetadataAddImageMask(psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "user-defined mask", maskBad);
     19    // select the output readout
     20    pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT");
    3621
    37     psImage *image = ro->image, *mask = ro->mask; // Image and mask of interest
     22    // select the model readout, if it exist already; if not, generate it
     23    pmReadout *modelRO = pmFPAfileThisReadout(config->files, view, "PSPHOT.BACKMDL");
    3824
    39     pmReadout *modelRO = pmFPAviewThisReadout(view, modelFile->fpa); // Background model
     25    // if necessary, generate the background model
    4026    if (!modelRO) {
    4127        // Create the background model
    … …  
    4430            return false;
    4531        }
    46         modelRO = (modelFile->mode == PM_FPA_MODE_INTERNAL) ? modelFile->readout :
    47                    pmFPAviewThisReadout(view, modelFile->fpa);
     32        // select the model readout (should now exist)
     33        modelRO = pmFPAfileThisReadout(config->files, view, "PSPHOT.BACKMDL");
    4834        if (!modelRO) {
    4935            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find background model");
    … …  
    5440                                                  "PSPHOT.BACKGROUND.BINNING"); // Binning for model
    5541    psImage *modelImage = modelRO->image; // Background model
     42
     43    psImage *image = outRO->image; // Image of interest
     44    psImage *mask = outRO->mask; // Mask of interest
    5645
    5746    // Do the background subtraction

  • trunk/ppSub/src/ppSubCamera.c

    r20614 r21257  
    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>
    9 
    101#include "ppSub.h"
    112
    123bool ppSubCamera(pmConfig *config)
    134{
    14     bool status = false;                // Status of definition
     5    bool status = false;                // result from pmFPAfileDefine operations
    156
    167    // Input image
    … …  
    119110    }
    120111
     112    // Convolved input image
     113    pmFPAfile *inConv = pmFPAfileDefineFromFile(config, input, 1, 1, "PPSUB.INPUT.CONV");
     114    if (!inConv) {
     115        psError(PS_ERR_IO, false, _("Unable to generate output file for PPSUB.INPUT.CONV"));
     116        return false;
     117    }
     118    if (output->type != PM_FPA_FILE_IMAGE) {
     119        psError(PS_ERR_IO, true, "PPSUB.INPUT.CONV is not of type IMAGE");
     120        return false;
     121    }
     122    // XXX should be based on recipe : inConv->save = true;
     123
     124    // Convolved input mask
     125    pmFPAfile *inConvMask = pmFPAfileDefineOutput(config, inConv->fpa, "PPSUB.INPUT.CONV.MASK");
     126    if (!inConvMask) {
     127        psError(PS_ERR_IO, false, _("Unable to generate output file from PPSUB.INPUT.CONV.MASK"));
     128        return false;
     129    }
     130    if (inConvMask->type != PM_FPA_FILE_MASK) {
     131        psError(PS_ERR_IO, true, "PPSUB.INPUT.CONV.MASK is not of type MASK");
     132        return false;
     133    }
     134    // XXX should be based on recipe : inConvMask->save = true;
     135
     136    // Convolved input weight
     137    if (inputWeight) {
     138        pmFPAfile *inConvWeight = pmFPAfileDefineOutput(config, inConv->fpa, "PPSUB.INPUT.CONV.WEIGHT");
     139        if (!inConvWeight) {
     140            psError(PS_ERR_IO, false, _("Unable to generate output file from PPSUB.OUTPUT.WEIGHT"));
     141            return false;
     142        }
     143        if (inConvWeight->type != PM_FPA_FILE_WEIGHT) {
     144            psError(PS_ERR_IO, true, "PPSUB.INPUT.CONV.WEIGHT is not of type WEIGHT");
     145            return false;
     146        }
     147        // XXX should be based on recipe : inConvWeight->save = true;
     148    }
     149
     150    // Convolved ref image
     151    pmFPAfile *refConv = pmFPAfileDefineFromFile(config, ref, 1, 1, "PPSUB.REF.CONV");
     152    if (!refConv) {
     153        psError(PS_ERR_IO, false, _("Unable to generate output file for PPSUB.REF.CONV"));
     154        return false;
     155    }
     156    if (output->type != PM_FPA_FILE_IMAGE) {
     157        psError(PS_ERR_IO, true, "PPSUB.REF.CONV is not of type IMAGE");
     158        return false;
     159    }
     160    // XXX should be based on recipe : refConv->save = true;
     161
     162    // Convolved ref mask
     163    pmFPAfile *refConvMask = pmFPAfileDefineOutput(config, refConv->fpa, "PPSUB.REF.CONV.MASK");
     164    if (!refConvMask) {
     165        psError(PS_ERR_IO, false, _("Unable to generate output file from PPSUB.REF.CONV.MASK"));
     166        return false;
     167    }
     168    if (refConvMask->type != PM_FPA_FILE_MASK) {
     169        psError(PS_ERR_IO, true, "PPSUB.REF.CONV.MASK is not of type MASK");
     170        return false;
     171    }
     172    // XXX should be based on recipe : refConvMask->save = true;
     173
     174    // Convolved ref weight
     175    if (refWeight) {
     176        pmFPAfile *refConvWeight = pmFPAfileDefineOutput(config, refConv->fpa, "PPSUB.REF.CONV.WEIGHT");
     177        if (!refConvWeight) {
     178            psError(PS_ERR_IO, false, _("Unable to generate output file from PPSUB.OUTPUT.WEIGHT"));
     179            return false;
     180        }
     181        if (refConvWeight->type != PM_FPA_FILE_WEIGHT) {
     182            psError(PS_ERR_IO, true, "PPSUB.REF.CONV.WEIGHT is not of type WEIGHT");
     183            return false;
     184        }
     185        // XXX should be based on recipe : refConvWeight->save = true;
     186    }
     187
    121188    // Output JPEGs
    122189    pmFPAfile *jpeg1 = pmFPAfileDefineOutput(config, NULL, "PPSUB.OUTPUT.JPEG1");
    … …  
    192259        }
    193260        pmFPAfileActivate(config->files, false, "PSPHOT.PSF.LOAD");
     261
     262        pmFPAfile *psphotResid = pmFPAfileDefineOutputForFormat (config, NULL, "PSPHOT.RESID", output->cameraName, output->formatName);
     263        if (!psphotResid) {
     264            psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.RESID");
     265            return false;
     266        }
     267        // make this optional: psphotResid->save = true;
    194268    }
    195269
    196270    // Always do this, since we're using psphot for the background model.
    197     {
     271    // XXX Not certain that I need to generate a separate pmFPAfile for PSPHOT.INPUT
     272    if (0) {
    198273        pmFPAfile *psphot = pmFPAfileDefineFromFPA(config, output->fpa, 1, 1, "PSPHOT.INPUT");
    199274        if (!psphot) {

  • trunk/ppSub/src/ppSubLoop.c

    r20117 r21257  
    1 #ifdef HAVE_CONFIG_H
    2 #include <config.h>
    3 #endif
    4 
    5 #include <stdio.h>
    6 #include <string.h>
    7 #include <pslib.h>
    8 #include <psmodules.h>
    9 #include <ppStats.h>
    10 
    111#include "ppSub.h"
    122
    … …  
    2818        }
    2919        psFree(resolved);
    30     }
    31 
    32     if (!pmConfigMaskSetBits(NULL, NULL, config)) {
    33         psError(PS_ERR_UNKNOWN, false, "Unable to determine mask value.");
    34         return false;
    3520    }
    3621
    … …  
    141126                    return false;
    142127                }
    143                 ppStatsFPA(stats, output->fpa, view, pmConfigMaskGet("MASK.VALUE", config), config);
     128                psImageMaskType maskValue = pmConfigMaskGet("MASK.VALUE", config);
     129                ppStatsFPA(stats, output->fpa, view, maskValue, config);
    144130            }
    145131

  • trunk/ppSub/src/ppSubReadout.c

    r21183 r21257  
    1 #ifdef HAVE_CONFIG_H
    2 #include <config.h>
    3 #endif
    4 
    5 #include <stdio.h>
    6 #include <string.h>
    7 #include <pslib.h>
    8 #include <psmodules.h>
    9 #include <psphot.h>
    10 
    111#include "ppSub.h"
    12 
    13 #define WCS_TOLERANCE 0.001             // Tolerance for WCS
    14 //#define TESTING                         // For test output
    15 
    16 #define SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_SATSTAR | PM_SOURCE_MODE_BLEND | \
    17                      PM_SOURCE_MODE_BADPSF | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \
    18                      PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to input sources
    19 
    20 
    21 // Copy every instance of a single keyword from one metadata to another
    22 static bool metadataCopySingle(psMetadata *target, psMetadata *source, const char *name)
    23 {
    24     PS_ASSERT_METADATA_NON_NULL(target, false);
    25     PS_ASSERT_METADATA_NON_NULL(source, false);
    26     PS_ASSERT_STRING_NON_EMPTY(name, false);
    27 
    28     psString regex = NULL;      // Regular expression
    29     psStringAppend(&regex, "^%s$", name);
    30     psMetadataIterator *iter = psMetadataIteratorAlloc(source, PS_LIST_HEAD, regex); // Iterator
    31     psFree(regex);
    32     psMetadataItem *item;       // Item from iteration
    33     while ((item = psMetadataGetAndIncrement(iter))) {
    34         psMetadataAddItem(target, item, PS_LIST_TAIL, PS_META_DUPLICATE_OK);
    35     }
    36     psFree(iter);
    37 
    38     return true;
    39 }
    40 
    412
    423bool ppSubReadout(pmConfig *config, psMetadata *stats, const pmFPAview *view)
    434{
    445    psTimerStart("PPSUB_MATCH");
    45     pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); // Input readout
    46     pmReadout *refRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF"); // Reference readout
    47     pmReadout *sourcesRO = pmFPAfileThisReadout(config->files, view, "PPSUB.SOURCES"); // Readout with sources
    48     pmReadout *inConv = pmReadoutAlloc(NULL); // Convolved version of input
    49     pmReadout *refConv = pmReadoutAlloc(NULL); // Convolved version of reference
    50     pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT"); // Output cell
    51     pmReadout *outRO = pmReadoutAlloc(outCell); // Output readout: subtraction
    52     pmFPA *outFPA = outCell->parent->parent; // Output FPA
    53     pmHDU *outHDU = outFPA->hdu; // Output HDU
    54     if (!outHDU->header) {
    55         outHDU->header = psMetadataAlloc();
    56     }
    576
    58     psImage *input = inRO->image;       // Input image
    59     psImage *reference = refRO->image;  // Reference image
    60     PS_ASSERT_IMAGES_SIZE_EQUAL(input, reference, false);
    61 
    62     int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads
    63     if (threads > 0) {
    64         if (!psThreadPoolInit(threads)) {
    65             psError(PS_ERR_UNKNOWN, false, "Unable to setup %d threads", threads);
    66             return false;
    67         }
    68         pmSubtractionThreadsInit(inRO, refRO);
    69     }
    70 
    71     // Look up recipe values
    72     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim
    73     psAssert(recipe, "We checked this earlier, so it should be here.");
    74 
    75     const char *typeStr = psMetadataLookupStr(NULL, recipe, "KERNEL.TYPE"); // Kernel type
    76     psAssert(typeStr, "We put it here in ppSubArguments.c");
    77     pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Type of kernel
    78     if (type == PM_SUBTRACTION_KERNEL_NONE) {
    79         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unrecognised kernel type: %s", typeStr);
    80         psFree(inConv);
    81         psFree(refConv);
    82         psFree(outRO);
     7    if (!ppSubSetMasks (config, view)) {
     8        psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
    839        return false;
    8410    }
    8511
    86     bool mdok;                          // Status of MD lookup
    87     int size = psMetadataLookupS32(NULL, recipe, "KERNEL.SIZE"); // Kernel half-size
    88     int order = psMetadataLookupS32(NULL, recipe, "SPATIAL.ORDER"); // Spatial polynomial order
    89     float regionSize = psMetadataLookupF32(NULL, recipe, "REGION.SIZE"); // Size of iso-kernel regs
    90     float spacing = psMetadataLookupF32(NULL, recipe, "STAMP.SPACING"); // Typical stamp spacing
    91     int footprint = psMetadataLookupS32(NULL, recipe, "STAMP.FOOTPRINT"); // Stamp half-size
    92     float threshold = psMetadataLookupF32(NULL, recipe, "STAMP.THRESHOLD"); // Threshold for stmps
    93     int stride = psMetadataLookupS32(NULL, recipe, "STRIDE"); // Size of convolution patches
    94     int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations
    95     float rej = psMetadataLookupF32(NULL, recipe, "REJ"); // Rejection threshold
    96     float sys = psMetadataLookupF32(NULL, recipe, "SYS"); // Relative systematic error in kernel
    97     bool reverse = psMetadataLookupBool(NULL, config->arguments, "REVERSE"); // Reverse sense of subtraction?
    98     psVector *widths = psMetadataLookupPtr(NULL, recipe, "ISIS.WIDTHS"); // ISIS Gaussian widths
    99     psVector *orders = psMetadataLookupPtr(NULL, recipe, "ISIS.ORDERS"); // ISIS Polynomial orders
    100     int inner = psMetadataLookupS32(NULL, recipe, "INNER"); // Inner radius
    101     int ringsOrder = psMetadataLookupS32(NULL, recipe, "RINGS.ORDER"); // RINGS polynomial order
    102     int binning = psMetadataLookupS32(NULL, recipe, "SPAM.BINNING"); // Binning for SPAM kernel
    103     float penalty = psMetadataLookupF32(NULL, recipe, "PENALTY"); // Penalty for wideness
    104     psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.IN"); // Name of bits to mask going in
    105     psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
    106     psString maskPoorStr = psMetadataLookupStr(NULL, recipe, "MASK.POOR"); // Name of bits to mask for poor
    107     psImageMaskType maskPoor = pmConfigMaskGet(maskPoorStr, config); // Bits to mask for poor pixels
    108     psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
    109     psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
    110     float badFrac = psMetadataLookupF32(NULL, recipe, "BADFRAC"); // Maximum bad fraction
    111     const char *stampsName = psMetadataLookupStr(&mdok, config->arguments, "STAMPS"); // Filename for stamps
    112 
    113     bool optimum = psMetadataLookupBool(&mdok, recipe, "OPTIMUM"); // Derive optimum parameters?
    114     float optMin = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.MIN"); // Minimum width for search
    115     float optMax = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.MAX"); // Maximum width for search
    116     float optStep = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.STEP"); // Step for search
    117     float optThresh = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.TOL"); // Tolerance for search
    118     int optOrder = psMetadataLookupS32(&mdok, recipe, "OPTIMUM.ORDER"); // Order for search
    119     bool dual = psMetadataLookupBool(&mdok, recipe, "DUAL"); // Dual convolution?
    120 
    121     // Statistics for renormalisation
    122     psStatsOptions renormMean = psStatsOptionFromString(psMetadataLookupStr(NULL, recipe, "RENORM.MEAN"));
    123     psStatsOptions renormStdev = psStatsOptionFromString(psMetadataLookupStr(NULL, recipe, "RENORM.STDEV"));
    124     if (renormMean == PS_STAT_NONE || renormStdev == PS_STAT_NONE) {
    125         psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    126                 "Unable to parse renormalisation statistics from recipe.");
    127         return false;
    128     }
    129     int renormNum = psMetadataLookupS32(&mdok, recipe, "RENORM.NUM"); // Number of samples for renormalisation
    130     float renormWidth = psMetadataLookupS32(&mdok, recipe, "RENORM.WIDTH"); // Width of Gaussian phot
    131 
    132     psString interpModeStr = psMetadataLookupStr(&mdok, recipe, "INTERPOLATION"); // Interpolation mode
    133     psImageInterpolateMode interpMode = psImageInterpolateModeFromString(interpModeStr); // Interp
    134     if (interpMode == PS_INTERPOLATE_NONE) {
    135         psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unknown interpolation mode: %s", interpModeStr);
    136         return false;
    137     }
    138     float poorFrac = psMetadataLookupF32(&mdok, recipe, "POOR.FRACTION"); // Fraction for "poor"
    139 
    140     pmSubtractionMode subMode = dual ? PM_SUBTRACTION_MODE_DUAL : PM_SUBTRACTION_MODE_UNSURE; // Subtracn mode
    141 
    142     // Generate masks if they don't exist
    143     int numCols = input->numCols, numRows = input->numRows; // Image dimensions
    144     if (!inRO->mask) {
    145         if (psMetadataLookupBool(NULL, recipe, "MASK.GENERATE")) {
    146             pmReadoutSetMask(inRO, pmConfigMaskGet("SAT", config), pmConfigMaskGet("BAD", config));
    147         } else {
    148             inRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
    149             psImageInit(inRO->mask, 0);
    150         }
    151     }
    152     if (!refRO->mask) {
    153         if (psMetadataLookupBool(NULL, recipe, "MASK.GENERATE")) {
    154             pmReadoutSetMask(refRO, pmConfigMaskGet("SAT", config), pmConfigMaskGet("BAD", config));
    155         } else {
    156             refRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
    157             psImageInit(refRO->mask, 0);
    158         }
    159     }
    160 
    161     if (!pmReadoutMaskNonfinite(inRO, pmConfigMaskGet("SAT", config))) {
    162         psError(PS_ERR_UNKNOWN, false, "Unable to mask non-finite pixels in input.");
    163         return false;
    164     }
    165     if (!pmReadoutMaskNonfinite(refRO, pmConfigMaskGet("SAT", config))) {
    166         psError(PS_ERR_UNKNOWN, false, "Unable to mask non-finite pixels in reference.");
     12    if (!ppSubMatchPSFs (config, view)) {
     13        psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
    16714        return false;
    16815    }
    16916
    170     psVector *optWidths = NULL;         // Vector with FWHMs for optimum search
    171     if (optimum) {
    172         optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32);
    173     }
    174 
    175     psArray *sources = NULL;            // Sources in image
    176     if (sourcesRO) {
    177         sources = psMetadataLookupPtr(&mdok, sourcesRO->analysis, "PSPHOT.SOURCES");
    178     }
    179 
    180 #if 0
    181     // Interpolate over bad pixels, so the bad pixels don't explode
    182     if (!pmReadoutInterpolateBadPixels(inRO, maskVal, interpMode, poorFrac, maskPoor, maskBad)) {
    183         psError(PS_ERR_UNKNOWN, false, "Unable to interpolate bad pixels for input image.");
    184         return false;
    185     }
    186     if (!pmReadoutInterpolateBadPixels(refRO, maskVal, interpMode, poorFrac, maskPoor, maskBad)) {
    187         psError(PS_ERR_UNKNOWN, false, "Unable to interpolate bad pixels for reference image.");
    188         return false;
    189     }
    190     maskVal |= maskBad;
    191 #endif
    192 
    193     // Match the PSFs
    194     if (!pmSubtractionMatch(inConv, refConv, inRO, refRO, footprint, stride, regionSize, spacing, threshold,
    195                             sources, stampsName, type, size, order, widths, orders, inner, ringsOrder,
    196                             binning, penalty, optimum, optWidths, optOrder, optThresh, iter, rej, sys,
    197                             maskVal, maskBad, maskPoor, poorFrac, badFrac, subMode)) {
    198         psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
    199         psFree(inConv);
    200         psFree(refConv);
    201         psFree(outRO);
    202         return false;
    203     }
    204     psFree(optWidths);
    205 
    206     pmSubtractionThreadsFinalize(inRO, refRO);
    207 
    208     // Add kernel descrption to header
    209     pmSubtractionKernels *kernels = psMetadataLookupPtr(&mdok, inConv->analysis,
    210                                                         PM_SUBTRACTION_ANALYSIS_KERNEL); // Subtraction kernel
    211     if (!kernels) {
    212         kernels = psMetadataLookupPtr(&mdok, refConv->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL);
    213     }
    214     if (!kernels) {
    215         psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find SUBTRACTION.KERNEL");
    216         psFree(inConv);
    217         psFree(refConv);
    218         psFree(outRO);
    219         return false;
    220     }
    221     psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.KERNEL", 0,
    222                      "Subtraction kernel", kernels->description);
    223 
    224     {
    225         if (psMetadataLookup(inConv->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL)) {
    226             outRO->analysis = psMetadataCopy(outRO->analysis, inConv->analysis);
    227         } else if (psMetadataLookup(refConv->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL)) {
    228             outRO->analysis = psMetadataCopy(outRO->analysis, refConv->analysis);
    229         } else {
    230             psWarning("Unable to find subtraction kernel in analysis metadata.");
    231         }
    232 
    233         // Update variance factors
    234         // It's not possible to do this perfectly, since we're combining different images:
    235         // vf_sum sigma_sum^2 = vf_1 * sigma_1^2 + vf_2 * sigma_2^2
    236         // It's not possible to write vf_sum as a function of vf_1 and vf_2 with no reference to the sigmas.
    237         // Instead, we're going to cheat.
    238         bool mdok;                      // Status of MD lookup
    239         pmSubtractionKernels *kernels = psMetadataLookupPtr(&mdok, outRO->analysis,
    240                                                             PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels
    241         if (!mdok) {
    242             psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find subtraction kernels.");
    243             psFree(inConv);
    244             psFree(refConv);
    245             psFree(outRO);
    246             return false;
    247         }
    248         float vfIn = psMetadataLookupF32(NULL, inRO->parent->concepts,
    249                                          "CELL.VARFACTOR"); // Variance factor for input
    250         if (!isfinite(vfIn)) {
    251             vfIn = 1.0;
    252         }
    253         float vfRef = psMetadataLookupF32(NULL, refRO->parent->concepts,
    254                                           "CELL.VARFACTOR"); // Variance factor for ref
    255         if (!isfinite(vfRef)) {
    256             vfRef = 1.0;
    257         }
    258         if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
    259             vfIn *= pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false);
    260         }
    261         if (kernels->mode == PM_SUBTRACTION_MODE_2) {
    262             vfRef *= pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false);
    263         } else if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
    264             vfRef *= pmSubtractionVarianceFactor(kernels, 0.0, 0.0, true);
    265         }
    266 
    267         psStats *vfStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics
    268         psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
    269         if (!psImageBackground(vfStats, NULL, inConv->weight, inConv->mask, maskVal | maskBad, rng)) {
    270             psError(PS_ERR_UNKNOWN, false, "Unable to measure mean variance for convolved input");
    271             psFree(inConv);
    272             psFree(refConv);
    273             psFree(outRO);
    274             psFree(vfStats);
    275             psFree(rng);
    276             return false;
    277         }
    278         float inMeanVar = vfStats->robustMedian; // Mean variance of input
    279         if (!psImageBackground(vfStats, NULL, refConv->weight, refConv->mask, maskVal | maskBad, rng)) {
    280             psError(PS_ERR_UNKNOWN, false, "Unable to measure mean variance for convolved reference");
    281             psFree(inConv);
    282             psFree(refConv);
    283             psFree(outRO);
    284             psFree(vfStats);
    285             psFree(rng);
    286             return false;
    287         }
    288         float refMeanVar = vfStats->robustMedian; // Mean variance of reference
    289         psFree(rng);
    290         psFree(vfStats);
    291 
    292         float vf = (vfIn * inMeanVar + vfRef * refMeanVar) / (inMeanVar + refMeanVar); // Resulting var factor
    293         psMetadataItem *vfItem = psMetadataLookup(outRO->parent->concepts, "CELL.VARFACTOR");
    294         vfItem->data.F32 = vf;
    295 
    296         // Statistics on the matching
    297         if (stats) {
    298             metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MODE);
    299             metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_STAMPS);
    300             metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_DEV_MEAN);
    301             metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_DEV_RMS);
    302             metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_NORM);
    303             metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_BGDIFF);
    304             metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MX);
    305             metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MY);
    306             metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MXX);
    307             metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MXY);
    308             metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MYY);
    309 
    310             psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_MATCH", 0, "Time to match PSFs",
    311                              psTimerClear("PPSUB_MATCH"));
    312         }
    313     }
    314 
    315 
    316 #ifdef TESTING
    317     psImage *kernelImage = psMetadataLookupPtr(&mdok, inConv->analysis,
    318                                                "SUBTRACTION.KERNEL.IMAGE"); // Image of the kernels
    319     if (!kernelImage) {
    320         kernelImage = psMetadataLookupPtr(&mdok, refConv->analysis, "SUBTRACTION.KERNEL.IMAGE");
    321     }
    322     psFits *fits = psFitsOpen("kernel.fits", "w");
    323     psFitsWriteImage(fits, NULL, kernelImage, 0, NULL);
    324     psFitsClose(fits);
    325 #endif
    326 
    327     // Subtraction is: minuend - subtrahend
    328     pmReadout *minuend = inConv;
    329     pmReadout *subtrahend = refConv;
    330 
    331     if (reverse) {
    332         pmReadout *temp = subtrahend;
    333         subtrahend = minuend;
    334         minuend = temp;
    335     }
    336 
    337 #ifdef TESTING
    338     {
    339         pmReadoutMaskApply(minuend, maskVal);
    340         psFits *fits = psFitsOpen("minuend.fits", "w");
    341         psFitsWriteImage(fits, NULL, minuend->image, 0, NULL);
    342         psFitsClose(fits);
    343     }
    344     {
    345         pmReadoutMaskApply(subtrahend, maskVal);
    346         psFits *fits = psFitsOpen("subtrahend.fits", "w");
    347         psFitsWriteImage(fits, NULL, subtrahend->image, 0, NULL);
    348             psFitsClose(fits);
    349     }
    350 #endif
    351 
    352     // Photometry is to be performed in two stages:
    353     // 1. Measure the PSF using the PSF-matched images
    354     // 2. Find and measure sources on the subtracted image
    355     psTimerStart("PPSUB_PHOT");
    356 
    357     // Photometry stage 1: measure the PSF
    358     pmPSF *psf = NULL;                  // PSF for photometry
    359     if (psMetadataLookupBool(NULL, recipe, "PHOTOMETRY")) {
    360         outRO->image = psImageCopy(outRO->image, minuend->image, PS_TYPE_F32);
    361         if (minuend->weight) {
    362             outRO->weight = psImageCopy(outRO->weight, minuend->weight, PS_TYPE_F32);
    363         }
    364         if (minuend->mask) {
    365             outRO->mask = psImageCopy(outRO->mask, minuend->mask, PS_TYPE_IMAGE_MASK);
    366         }
    367         outRO->data_exists = outCell->data_exists = outCell->parent->data_exists = true;
    368 
    369         if (psMetadataLookupBool(&mdok, recipe, "RENORM")) {
    370             psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask
    371             if (!pmReadoutWeightRenormPhot(outRO, maskValue, renormNum, renormWidth,
    372                                            renormMean, renormStdev, NULL)) {
    373                 psError(PS_ERR_UNKNOWN, false, "Unable to renormalise variances.");
    374                 psFree(outRO);
    375                 return false;
    376             }
    377         }
    378 
    379         pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT");
    380         pmFPACopy(photFile->fpa, outRO->parent->parent->parent);
    381 
    382         // We have a list of sources, so we could use those to redetermine the PSF model.
    383         // Until Gene makes the necessary adaptations to psphot, we will simply redetermine the PSF model from
    384         // scratch.
    385 
    386         if (psMetadataLookupBool(&mdok, recipe, "PSPHOT.VISUAL")) {
    387             psphotSetVisual(true);
    388         }
    389 
    390         // Need to ensure aperture residual is not calculated
    391         psMetadata *psphotRecipe = psMetadataLookupMetadata(NULL, config->recipes, PSPHOT_RECIPE); // Recipe
    392         if (!psphotRecipe) {
    393             psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find %s recipe.", PSPHOT_RECIPE);
    394             return false;
    395         }
    396         const char *breakpoint = psMetadataLookupStr(&mdok, psphotRecipe, "BREAK_POINT"); // Current break point
    397         psMetadataAddStr(psphotRecipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE,
    398                          "ALTERED break point for psphot operations", "PSFMODEL");
    399 
    400         // set maskValue and markValue in the psphot recipe
    401         psImageMaskType maskValue = maskVal;
    402         psImageMaskType markValue = pmConfigMaskGet("MARK.VALUE", config); // Bits to use for marking
    403         psMetadataAddImageMask(psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "Bits to mask", maskValue);
    404         psMetadataAddImageMask(psphotRecipe, PS_LIST_TAIL, "MARK.PSPHOT", PS_META_REPLACE, "Bits to use for marking",
    405                         markValue);
    406 
    407         if (!psphotReadout(config, view)) {
    408             psWarning("Unable to perform photometry on subtracted image.");
    409             psErrorStackPrint(stderr, "Error stack from photometry:");
    410             psErrorClear();
    411         }
    412 
    413         if (!pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL") ||
    414             !pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV") ||
    415             !pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND")) {
    416             psError(PS_ERR_UNKNOWN, false, "Unable to drop PSPHOT internal files.");
    417             return false;
    418         }
    419 
    420         if (breakpoint) {
    421             psMetadataAddStr(psphotRecipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE,
    422                              "RESTORED break point for psphot operations", breakpoint);
    423         }
    424 
    425         // Blow away the sources psphot found --- they're irrelevant for the subtraction
    426         pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout with sources
    427         psMetadataRemoveKey(photRO->analysis, "PSPHOT.SOURCES");
    428         psMetadataRemoveKey(photRO->analysis, "PSPHOT.HEADER");
    429 
    430         psf = psMemIncrRefCounter(psMetadataLookupPtr(NULL, photRO->parent->parent->analysis, "PSPHOT.PSF"));
    431         if (!psf) {
    432             psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find PSF from psphot");
    433             return false;
    434         }
    435 
    436         // Record the FWHM
    437         pmHDU *hdu = pmHDUFromCell(outRO->parent); // HDU with header
    438         psMetadataItemSupplement(hdu->header, psphotRecipe, "FWHM_MAJ");
    439         psMetadataItemSupplement(hdu->header, psphotRecipe, "FWHM_MIN");
    440 
    441 
    442         pmCell *photCell = pmFPAfileThisCell(config->files, view, "PSPHOT.INPUT");
    443         pmCellFreeReadouts(photCell);
    444     }
    445 
    446     // Do the actual subtraction
    447     outRO->image = (psImage*)psBinaryOp(outRO->image, minuend->image, "-", subtrahend->image);
    448     if (inConv->weight && refConv->weight) {
    449         outRO->weight = (psImage*)psBinaryOp(outRO->weight, inConv->weight, "+", refConv->weight);
    450     }
    451     outRO->mask = (psImage*)psBinaryOp(outRO->mask, inConv->mask, "|", refConv->mask);
    452 
    453     psFree(inConv);
    454     psFree(refConv);
    455 
    456     outRO->data_exists = outCell->data_exists = outCell->parent->data_exists = true;
    457 
    458     // Copy astrometry over
    459     // It should get into the output images and photometry
    460     pmFPA *refFPA = refRO->parent->parent->parent; // Reference FPA
    461     pmHDU *refHDU = refFPA->hdu;        // Reference HDU
    462     if (!outHDU || !refHDU) {
    463         psWarning("Unable to find HDU at FPA level to copy astrometry.");
    464     } else if (!pmAstromReadWCS(outFPA, outCell->parent, refHDU->header, 1.0)) {
    465         psWarning("Unable to read WCS astrometry from reference FPA.");
    466         psErrorClear();
    467     } else if (!pmAstromWriteWCS(outHDU->header, outFPA, outCell->parent, WCS_TOLERANCE)) {
    468         psWarning("Unable to write WCS astrometry to output FPA.");
    469         psErrorClear();
    470     }
    471 
    472 #if 0
    473     pmReadoutMaskApply(outRO, maskBad);
    474 #endif
    475 
    476 #ifdef TESTING
    477     for (int y = 0; y < outRO->image->numRows; y++) {
    478         for (int x = 0; x < outRO->image->numCols; x++) {
    479             if (isnan(outRO->image->data.F32[y][x]) && !(outRO->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) {
    480                 printf("Unmasked NAN at %d %d --> %d\n", x, y, outRO->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x]);
    481             }
    482         }
    483     }
    484 #endif
    485 
    486     if (!pmFPACopyConcepts(outCell->parent->parent, inRO->parent->parent->parent)) {
    487         psError(PS_ERR_UNKNOWN, false, "Unable to copy concepts from input to output.");
    488         psFree(outRO);
     17    if (!ppSubDefineOutput (config, view)) {
     18        psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
    48919        return false;
    49020    }
    49121
    492 #if 0
    493     // XXX Under development
     22    if (!ppSubVarianceFactors (config, stats, view)) {
     23        psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
     24        return false;
     25    }
    49426
    495     // QA: photometry of known sources with measured PSF
    496     if (sourcesRO && psMetadataLookupBool(NULL, recipe, "PHOTOMETRY")) {
    497         sources = psMetadataLookupPtr(&mdok, sourcesRO->analysis, "PSPHOT.SOURCES");
     27    if (!ppSubMakePSF(config, view)) {
     28        psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
     29        return false;
     30    }
    49831
    499 
    500 
    501         // For each source, need to:
    502         // * generate appropriate model
    503         // * select pixels
    504         // * define fit radius
    505         // Don't need to measure sky (psphotFitSourcesLinear does that)
    506 
    507         psArray *dummy = psArrayAlloc(1); // Dummy array of sources, for psphotFitSourcesLinear with 1 source
    508         for (long i = 0; i < sources->n; i++) {
    509             pmSource *source = sources->data[i];
    510             if (!source || (source->mode & SOURCE_MASK)) {
    511                 continue;
    512             }
    513 
    514             float origMag = source->psfMag; // Original magnitude
    515 
    516             // Coordinates of source
    517             float x = source->modelPSF->params->data.F32[PM_PAR_XPOS];
    518             float y = source->modelPSF->params->data.F32[PM_PAR_YPOS];
    519 
    520             model = pmModelFromPSFforXY(psf, x, y, NAN);
    521             if (!fakeModel) {
    522                 psWarning("Can't generate model");
    523                 continue;
    524             }
    525             if (!pmSourceDefinePixels(fakeSource, readout, x, y, fakeRadius)) {
    526                 psErrorClear();
    527                 continue;
    528             }
    529 pmModel *pmModelFromPSFforXY (
    530     const pmPSF *psf,
    531     float Xo,
    532     float Yo,
    533     float Io
    534     );
    535 
    536 bool pmSourceDefinePixels(
    537     pmSource *mySource,                 ///< source to be re-defined
    538     const pmReadout *readout,  ///< base the source on this readout
    539     psF32 x,                            ///< center coords of source
    540     psF32 y,                            ///< center coords of source
    541     psF32 Radius                        ///< size of box on source
    542 );
    543 
    544 
    545 bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) {
    546 
    547     }
    548 #endif
    549 
    550     // Photometry stage 2: find and measure sources on the subtracted image
    551     if (psMetadataLookupBool(NULL, recipe, "PHOTOMETRY")) {
    552         // The PSF should already be stored for the readout
    553         pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT");
    554         pmFPACopy(photFile->fpa, outFPA);
    555 
    556         pmReadout *psfRO = pmFPAfileThisReadout(config->files, view, "PSPHOT.PSF.LOAD");
    557         if (!psfRO) {
    558             psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find readout for file PSPHOT.PSF.LOAD");
    559             return false;
    560         }
    561         psMetadataAddPtr(psfRO->parent->parent->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN,
    562                          "PSF from matched addition", psf);
    563         psFree(psf);
    564 
    565         if (psMetadataLookupBool(&mdok, recipe, "PSPHOT.VISUAL")) {
    566             psphotSetVisual(true);
    567         }
    568 
    569         if (psMetadataLookupBool(&mdok, recipe, "RENORM")) {
    570             psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask
    571             if (!pmReadoutWeightRenormPhot(outRO, maskValue, renormNum, renormWidth,
    572                                            renormMean, renormStdev, NULL)) {
    573                 psError(PS_ERR_UNKNOWN, false, "Unable to renormalise variances.");
    574                 psFree(outRO);
    575                 return false;
    576             }
    577         }
    578 
    579         // Need to ensure aperture residual is not calculated
    580         psMetadata *psphotRecipe = psMetadataLookupMetadata(NULL, config->recipes, PSPHOT_RECIPE); // Recipe
    581         psMetadataItem *item = psMetadataLookup(psphotRecipe, "MEASURE.APTREND"); // Item determining aptrend
    582         if (!item) {
    583             psWarning("Unable to find MEASURE.APTREND in psphot recipe");
    584             psErrorClear();
    585         } else {
    586             item->data.B = false;
    587         }
    588 
    589         if (!psphotReadout(config, view)) {
    590             psWarning("Unable to perform photometry on subtracted image.");
    591             psErrorStackPrint(stderr, "Error stack from photometry:");
    592             psErrorClear();
    593         }
    594 
    595         if (!pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL") ||
    596             !pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV") ||
    597             !pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND")) {
    598             psError(PS_ERR_UNKNOWN, false, "Unable to drop PSPHOT internal files.");
    599             return false;
    600         }
    601 
    602         pmFPAfileActivate(config->files, false, "PSPHOT.INPUT");
    603         pmFPAfileActivate(config->files, false, "PSPHOT.LOAD.PSF");
    604 
    605         if (stats) {
    606             pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout with the sources
    607             psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources
    608             psMetadataAddS32(stats, PS_LIST_TAIL, "NUM_SOURCES", 0, "Number of sources detected",
    609                              sources ? sources->n : 0);
    610             psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_PHOT", 0, "Time to do photometry",
    611                              psTimerClear("PPSUB_PHOT"));
    612         }
    613 
    614 #ifdef TESTING
    615         // Record data about sources: not everything gets into the output CMF files
    616         {
    617             pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout with the sources
    618             psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources
    619             FILE *sourceFile = fopen("sources.dat", "w"); // File for sources
    620             fprintf(sourceFile,
    621                     "# x y mag mag_err psf_chisq cr_nsigma ext_nsigma psf_qf flags m_x m_y m_xx m_xy m_yy\n");
    622             for (int i = 0; i < sources->n; i++) {
    623                 pmSource *source = sources->data[i];
    624                 if (!source) {
    625                     continue;
    626                 }
    627 
    628                 float x, y;             // Position of source
    629                 float chi2;             // chi^2 for source
    630                 if (source->modelPSF) {
    631                     x = source->modelPSF->params->data.F32[PM_PAR_XPOS];
    632                     y = source->modelPSF->params->data.F32[PM_PAR_YPOS];
    633                     chi2 = source->modelPSF->chisq;
    634                 } else if (source->peak) {
    635                     x = source->peak->xf;
    636                     y = source->peak->yf;
    637                     chi2 = NAN;
    638                 } else {
    639                     psWarning("No position available for source.");
    640                     continue;
    641                 }
    642 
    643                 float xMoment = NAN, yMoment = NAN, xxMoment = NAN, xyMoment = NAN, yyMoment = NAN;
    644                 if (source->moments) {
    645                     xMoment = source->moments->Mx;
    646                     yMoment = source->moments->My;
    647                     xxMoment = source->moments->Mxx;
    648                     xyMoment = source->moments->Mxy;
    649                     yyMoment = source->moments->Myy;
    650                 }
    651 
    652                 fprintf(sourceFile, "%f %f %f %f %f %f %f %f %d %f %f %f %f %f\n",
    653                         x, y, source->psfMag, source->errMag, chi2, source->crNsigma, source->extNsigma,
    654                         source->pixWeight, source->mode, xMoment, yMoment, xxMoment, xyMoment, yyMoment);
    655             }
    656             fclose(sourceFile);
    657         }
    658 #endif
    659 
     32    if (!ppSubReadoutSubtract (config, view)) {
     33        psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
     34        return false;
    66035    }
    66136
    66237    // Higher order background subtraction using psphot
    663     if (!ppSubBackground(outRO, config, view)) {
     38    if (!ppSubBackground(config, view)) {
    66439        psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
    665         psFree(outRO);
     40        return false;
     41    }
     42 
     43    if (!ppSubReadoutPhotometry (config, stats, view)) {
     44        psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
    66645        return false;
    66746    }
    66847
    669     // Renormalising for pixels, because that's what magic desires
    670     if (psMetadataLookupBool(&mdok, recipe, "RENORM")) {
    671         psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask
    672         if (!pmReadoutWeightRenormPixels(outRO, maskValue, renormMean, renormStdev, NULL)) {
    673             psError(PS_ERR_UNKNOWN, false, "Unable to renormalise variances.");
    674             psFree(outRO);
    675             return false;
    676         }
     48    if (!ppSubReadoutUpdate (config, view)) {
     49        psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
     50        return false;
    67751    }
    678 
    679     // Add additional data to the header
    680     pmFPAfile *refFile = psMetadataLookupPtr(NULL, config->files, "PPSUB.REF"); // Reference file
    681     pmFPAfile *inFile = psMetadataLookupPtr(NULL, config->files, "PPSUB.INPUT"); // Input file
    682     psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.REFERENCE", 0,
    683                      "Subtraction reference", refFile->filename);
    684     psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.INPUT", 0,
    685                      "Subtraction input", inFile->filename);
    686 
    687 
    688     // Generate binned JPEGs
    689     {
    690         int bin1 = psMetadataLookupS32(NULL, recipe, "BIN1"); // First binning level
    691         int bin2 = psMetadataLookupS32(NULL, recipe, "BIN2"); // Second binning level
    692 
    693         // Target cells
    694         pmCell *cell1 = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT.JPEG1");
    695         pmCell *cell2 = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT.JPEG2");
    696 
    697         pmReadout *ro1 = pmReadoutAlloc(cell1), *ro2 = pmReadoutAlloc(cell2); // Binned readouts
    698         if (!pmReadoutRebin(ro1, outRO, maskBad, bin1, bin1) || !pmReadoutRebin(ro2, ro1, 0, bin2, bin2)) {
    699             psError(PS_ERR_UNKNOWN, false, "Unable to bin output.");
    700             psFree(ro1);
    701             psFree(ro2);
    702             psFree(outRO);
    703             return false;
    704         }
    705         psFree(ro1);
    706         psFree(ro2);
    707     }
    708 
    709 #ifdef TESTING
    710     // Significance image
    711     {
    712         psImage *sig = (psImage*)psBinaryOp(NULL, outRO->image, "*", outRO->image);
    713         psBinaryOp(sig, sig, "/", outRO->weight);
    714         psFits *fits = psFitsOpen("significance.fits", "w");
    715         psFitsWriteImage(fits, NULL, sig, 0, NULL);
    716         psFitsClose(fits);
    717         psFree(sig);
    718     }
    719 #endif
    720 
    721     psFree(outRO);
    72252
    72353    return true;

  • trunk/ppSub/src/ppSubVersion.c

    r13341 r21257  
    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 <ppStats.h>
    9 
    101#include "ppSub.h"
    112
Note: See TracChangeset for help on using the changeset viewer.