IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5717


Ignore:
Timestamp:
Dec 6, 2005, 6:43:52 PM (21 years ago)
Author:
Paul Price
Message:

Importing current working PAP version (many bug fixes since this branch) straight on top of this version, and tidying up a bit

Location:
trunk/pois
Files:
17 edited

Legend:

Unmodified
Added
Removed

  • trunk/pois/Makefile.am

    r3792 r5717  
    1 SUBDIRS = src swig
     1SUBDIRS = src #swig
    22
    33CLEANFILES = *.pyc *~ core core.*

  • trunk/pois/configure.ac

    r3867 r5717  
    1515AC_PROG_LIBTOOL
    1616
    17 PKG_CHECK_MODULES(PSLIB, pslib >= 0.5.0)
     17PKG_CHECK_MODULES(PSLIB, pslib >= 0.8.0)
    1818
    19 dnl is this the best was to setup recursive CFLAGS?
    20 pois_CFLAGS="-Wall -g -O3 -std=c99"
     19dnl Set CFLAGS for build
     20pois_CFLAGS="-Wall -Werror -g -O3 -std=c99"
    2121AC_SUBST([pois_CFLAGS])
    2222
    2323AC_C_BIGENDIAN([PS_BIGENDIAN=1],[PS_BIGENDIAN=0])
    24                
     24
    2525AC_PROG_SWIG(1.3.24)
    2626SWIG_MULTI_MODULE_SUPPORT
    … …  
    2929dnl lifted from autoconf 2.59's configure.ac
    3030AC_PATH_PROG([PERL], perl, no)
    31 AC_SUBST([PERL])dnl
     31AC_SUBST([PERL])
    3232if test "$PERL" = no; then
    3333  AC_MSG_ERROR([perl is not found])

  • trunk/pois/src/Makefile.am

    r3830 r5717  
    1 AM_CFLAGS = $(pois_CFLAGS) $(PSLIB_CFLAGS)
     1AM_CFLAGS = $(PSLIB_CFLAGS) $(pois_CFLAGS)
    22AM_LDFLAGS = $(PSLIB_LIBS)
     3
     4pois_CPPFLAGS = $(PSLIB_CFLAGS) $(pois_CFLAGS)
     5pois_LDFLAGS = $(PSLIB_LIBS)
    36
    47include_HEADERS = pois.h poisErrorCodes.h
    … …  
    1922libpois_la_SOURCES = \
    2023        pois.h \
    21         poisParseConfig.c \
    22         poisKernelBasisFunctions.c \
     24        poisAddPenalty.c \
     25        poisCalculateDeviations.c \
    2326        poisCalculateEquation.c \
     27        poisCheckKernel.c \
    2428        poisConvolveImage.c \
    2529        poisExtractKernel.c \
     30        poisFindStamps.c \
     31        poisKernelBasisFunctions.c \
    2632        poisMakeMask.c \
    27         poisFindStamps.c \
    28         poisCalculateDeviations.c \
    29         poisAddPenalty.c \
     33        poisReadStamps.c \
     34        poisRejectStamps.c \
     35        poisParseConfig.c \
     36        poisSolveEquation.c \
    3037        poisStamp.c \
    31         poisRejectStamps.c \
    32         poisSolveEquation.c \
    33         poisMaskOps.c \
    3438        errorCodes.c
     39###     poisMaskOps.c
    3540
    3641#
    … …  
    4954        -$(RM) testout.fits stamp*.fits kernel*.fits
    5055        ./pois -v test1.fits test2.fits testout.fits
     56
     57# Tags for emacs
     58tags:
     59        etags `find . -name \*.[ch] -print`

  • trunk/pois/src/pois.c

    r3813 r5717  
    77
    88#include <stdio.h>
     9#include <string.h>
    910#include <sys/time.h>
    1011#include "pslib.h"
    … …  
    2526int main(int argc, char **argv)
    2627{
    27     poisErrorRegister();
    28 
    2928    double startTime = getTime();
    3029
    3130    // Set trace levels
    3231    psTraceSetLevel(".", 0);
    33     psTraceSetLevel("pois", 10);
    34     psTraceSetLevel("pois.config", 10);
    35     psTraceSetLevel("pois.time", 10);
    36     psTraceSetLevel("pois.solution", 0);
    37     psTraceSetLevel("pois.kernelBasisFunctions", 10);
    38     psTraceSetLevel("pois.calculateEquation", 10);
    39     psTraceSetLevel("pois.convolveImage", 10);
    40     psTraceSetLevel("pois.findStamps", 10);
    41     psTraceSetLevel("pois.parseConfig", 10);
    42     psTraceSetLevel("pois.makeMask", 10);
    43     psTraceSetLevel("pois.extractKernel", 10);
    44     psTraceSetLevel("pois.calculateDeviations", 10);
    45 
    46     // Set logging level
    47     psLogSetLevel(9);
    4832
    4933    // Parse the command line
    5034    poisConfig *config = poisParseConfig(argc, argv); // Configuration values
    5135
     36    if (config->verbose) {
     37        psTraceSetLevel("pois", 10);
     38        psTraceSetLevel("pois.config", 10);
     39        psTraceSetLevel("pois.time", 10);
     40        psTraceSetLevel("pois.solution", 0);
     41        psTraceSetLevel("pois.kernelBasisFunctions", 10);
     42        psTraceSetLevel("pois.calculateEquation", 10);
     43        psTraceSetLevel("pois.convolveImage", 10);
     44        psTraceSetLevel("pois.findStamps", 10);
     45        psTraceSetLevel("pois.parseConfig", 10);
     46        psTraceSetLevel("pois.makeMask", 10);
     47        psTraceSetLevel("pois.extractKernel", 10);
     48        psTraceSetLevel("pois.calculateDeviations", 10);
     49        psTraceSetLevel("pois.checkKernel", 10);
     50        if (config->stampFile) {
     51            psTraceSetLevel("pois.checkStamp", 10);
     52        } else {
     53            psTraceSetLevel("pois.checkStamp", 0);
     54        }
     55        // Set logging level
     56        psLogSetLevel(9);
     57    }
     58
    5259    // Read inputs
    53     psRegion *imageRegion = psRegionAlloc(0, 0, 0, 0);
     60    psMetadata *header = NULL;          // Header for output image
     61
     62    psRegion imageRegion = {0, 0, 0, 0};
    5463    psFits *reference = psFitsAlloc(config->refFile);
    55     //psMetadata *refHeader = psFitsReadHeader(NULL, reference);
    56     psImage *refImage = psFitsReadImage(NULL, reference, *imageRegion, 0);
     64    if (config->reverse) {
     65        header = psFitsReadHeader(NULL, reference);
     66    }
     67    psImage *refImage = psFitsReadImage(NULL, reference, imageRegion, 0);
    5768    if (refImage == NULL) {
    5869        psErrorStackPrint(stderr, "Fatal error: unable to read %s\n", config->refFile);
    … …  
    6172    psTrace("pois", 3, "Read %s: %dx%d\n", config->refFile, refImage->numCols, refImage->numRows);
    6273    psFree(reference);
     74    // Convert to 32-bit floating point, in necessary
     75    if (refImage->type.type != PS_TYPE_F32) {
     76        psTrace("pois", 3, "Converting %s to floating point in memory....\n", config->refFile);
     77        psImage *temp = psImageCopy(NULL, refImage, PS_TYPE_F32);
     78        psFree(refImage);
     79        refImage = temp;
     80    }
    6381
    6482    psFits *input = psFitsAlloc(config->inFile);
    65     //psMetadata *inHeader = psFitsReadHeader(NULL, input);
    66     psImage *inImage = psFitsReadImage(NULL, input, *imageRegion, 0);
     83    if (! config->reverse) {
     84        header = psFitsReadHeader(NULL, input);
     85    }
     86    psImage *inImage = psFitsReadImage(NULL, input, imageRegion, 0);
    6787    if (inImage == NULL) {
    6888        psErrorStackPrint(stderr, "Fatal error: unable to read %s\n", config->inFile);
    … …  
    7191    psTrace("pois", 3, "Read %s: %dx%d\n", config->inFile, inImage->numCols, inImage->numRows);
    7292    psFree(input);
    73     psFree(imageRegion);
     93    // Convert to 32-bit floating point, in necessary
     94    if (inImage->type.type != PS_TYPE_F32) {
     95        psTrace("pois", 3, "Converting %s to floating point in memory....\n", config->inFile);
     96        psImage *temp = psImageCopy(NULL, inImage, PS_TYPE_F32);
     97        psFree(inImage);
     98        inImage = temp;
     99    }
    74100
    75101    if ((refImage->numCols != inImage->numCols) || (refImage->numRows != inImage->numRows)) {
    … …  
    102128    // Write the mask out
    103129    psFits *maskFile = psFitsAlloc("mask.fits");
    104     if (!psFitsWriteImage(maskFile, NULL, mask, 0, NULL)) {
     130    if (!psFitsWriteImage(maskFile, NULL, mask, 0)) {
    105131        psErrorStackPrint(stderr, "Unable to write mask: mask.fits\n");
    106132    }
    … …  
    110136
    111137    psArray *stamps = NULL;             // Array of stamps
    112     psVector *stampMask = NULL;         // Mask for stamps
    113138    psVector *solution = NULL;          // Solution vector
    114139
     140    FILE *stampsFP = NULL;              // File pointer for stamps
     141    if (config->stampFile) {
     142        psTrace("pois", 5, "Opening stamp file, %s\n", config->stampFile);
     143        stampsFP = fopen(config->stampFile, "r");
     144        if (! stampsFP) {
     145            psError(PS_ERR_IO, true, "Unable to open stamps file, %s!\n", config->stampFile);
     146            exit(EXIT_FAILURE);
     147        }
     148    }
     149
    115150    // Iterate for a good solution
    116     psList *badStamps = NULL;           // Do we have bad stamps, such that we need to continue to iterate?
     151    bool badStamps = true;              // Do we have bad stamps, such that we need to continue to iterate?
    117152    for (int iterNum = 0; iterNum < config->numIter && badStamps; iterNum++) {
    118153        psTrace("pois", 1, "Iteration %d...\n", iterNum);
    119154
    120         psFree(badStamps);              // left over from last iteration
    121        
    122155        // Find stamps
    123         stamps = poisFindStamps(stamps, refImage, mask, config);
     156        if (config->stampFile) {
     157            stamps = poisReadStamps(stamps, &stampsFP, refImage, mask, config);
     158        } else {
     159            stamps = poisFindStamps(stamps, refImage, mask, config);
     160        }
    124161        int numStamps = 0;
    125162        for (int s = 0; s < stamps->n; s++) {
    … …  
    159196        psImage *kernelImage = poisExtractKernel(solution, kernelBasisFunctions, 0.0, 0.0, config);
    160197        psFits *kernelFile = psFitsAlloc(kernelName);
    161         if (!psFitsWriteImage(kernelFile, NULL, kernelImage, 0, NULL)) {
     198        if (!psFitsWriteImage(kernelFile, NULL, kernelImage, 0)) {
    162199            psErrorStackPrint(stderr, "Unable to write kernel: %s\n", kernelName);
    163200        }
    … …  
    177214    // If there was rejection on the last round, re-solve the equation using only the good stamps
    178215    if (badStamps) {
    179         psFree(badStamps);
    180216        solution = poisSolveEquation(solution, stamps, config);
     217    }
     218
     219    if (stampsFP) {
     220        fclose(stampsFP);
    181221    }
    182222
    … …  
    185225        for (int i = 0; i < mask->numCols; i++) {
    186226            if (mask->data.U8[j][i] & POIS_MASK_STAMP) {
    187                 mask->data.U8[j][i] = POIS_MASK_OK;
    188             }
    189         }
    190     }
     227                mask->data.U8[j][i] &= ~POIS_MASK_STAMP;
     228            }
     229        }
     230    }
     231
     232    // Check the result
     233    (void)poisCheckKernel(solution, kernelBasisFunctions, config);
    191234
    192235    // Convolve the input image
    … …  
    199242    snprintf(convName, MAXCHAR, "%s.conv", config->outFile);
    200243    psFits *convFile = psFitsAlloc(convName);
    201     if (!psFitsWriteImage(convFile, NULL, convImage, 0, NULL)) {
     244    if (!psFitsWriteImage(convFile, NULL, convImage, 0)) {
    202245        psErrorStackPrint(stderr, "Unable to write convolved image: %s\n", convName);
    203246    }
    … …  
    208251    // Do the subtraction
    209252    psImage *subImage = psImageAlloc(config->xImage, config->yImage, PS_TYPE_F32);
    210     (void)psBinaryOp(subImage, inImage, "-", convImage);
     253    if (config->reverse) {
     254        (void)psBinaryOp(subImage, convImage, "-", inImage);
     255    } else {
     256        (void)psBinaryOp(subImage, inImage, "-", convImage);
     257    }
    211258    psTrace("pois.time", 1, "Subtraction completed at %f sec\n", getTime() - startTime);
    212259
    … …  
    231278            }
    232279        }
    233        
    234         for (int y = 0; y < config->yKernel; y++) {
    235             for (int x = 0; x < subImage->numCols; x++) {
    236                 subImage->data.F32[y][x] = 0.0;
    237             }
    238         }
    239         for (int y = subImage->numRows - config->yKernel; y < subImage->numRows; y++) {
    240             for (int x = 0; x < subImage->numCols; x++) {
    241                 subImage->data.F32[y][x] = 0.0;
    242             }
    243         }
    244     }
     280    }
     281    for (int y = 0; y < config->yKernel; y++) {
     282        for (int x = 0; x < subImage->numCols; x++) {
     283            subImage->data.F32[y][x] = 0.0;
     284        }
     285    }
     286    for (int y = subImage->numRows - config->yKernel; y < subImage->numRows; y++) {
     287        for (int x = 0; x < subImage->numCols; x++) {
     288            subImage->data.F32[y][x] = 0.0;
     289        }
     290    }
     291
    245292
    246293    // Mask out bad pixels
    … …  
    253300    }
    254301
     302    if (config->mask) {
     303        char maskName[80];              // Name of mask image
     304        sprintf(maskName, "%s.mask", config->outFile);
     305        psFits *maskFile = psFitsAlloc(maskName);
     306        if (!psFitsWriteImage(maskFile, NULL, mask, 0)) {
     307            psErrorStackPrint(stderr, "Unable to write image: %s\n", maskName);
     308        }
     309        psTrace("pois", 1, "Mask image written to %s\n", maskName);
     310        psFree(maskFile);
     311    }
     312
     313
     314    int numNaN = psImageClipNaN(subImage, 0.0);
     315    if (numNaN > 0) {
     316        printf("Masked %d NAN pixels to zero.\n", numNaN);
     317    }
     318
    255319    psFits *subFile = psFitsAlloc(config->outFile);
    256     if (!psFitsWriteImage(subFile, NULL, subImage, 0, NULL)) {
     320    if (!psFitsWriteImage(subFile, header, subImage, 0)) {
    257321        psErrorStackPrint(stderr, "Unable to write subtracted image: %s\n", config->outFile);
    258322    }
    … …  
    261325    psTrace("pois.time", 1, "I/O completed at %f sec\n", getTime() - startTime);
    262326
    263 
    264327    // Clean up
    265     (void)psFree(kernelBasisFunctions);
    266     (void)psFree(refImage);
    267     (void)psFree(inImage);
    268     (void)psFree(convImage);
    269     (void)psFree(subImage);
    270     (void)psFree(config);
     328    psFree(header);
     329    psFree(kernelBasisFunctions);
     330    psFree(refImage);
     331    psFree(inImage);
     332    psFree(convImage);
     333    psFree(subImage);
     334    psFree(config);
    271335}
     336

  • trunk/pois/src/pois.h

    r3813 r5717  
    33
    44#include "pslib.h"
    5 #include "poisErrorCodes.h"
    65
    76#define abs(x) ((x) >= 0 ? (x) : (-(x)))
    … …  
    2120    int nsx, nsy;                       // Number of stamps in x and y
    2221    float refSat, inSat;                // Saturation values for reference and input images
    23     float bad;                          // Value below which pixels are bad
     22    float refBad, inBad;                // Bad values for reference and input images
    2423    float background;                   // Value to add to the background in images to get above zero
    2524    int spatialOrder;                   // Order of spatial variations in the kernel
    … …  
    2726    float sigmaRej;                     // Limit (in std dev) for rejection
    2827    float penalty;                      // Penalty value
     28    char *stampFile;                    // Filename for file with stamps
     29    bool reverse;                       // Reverse the subtraction?
     30    bool mask;                          // Output the mask?
    2931} poisConfig;
    3032
    … …  
    6466void poisStampFree(poisStamp *stamp);
    6567
     68// Return true if the stamp is OK
     69bool poisCheckStamp(const psImage *image, // Image to check for threshold
     70                    const psImage *mask, // Mask to check for bad pixels
     71                    int x, int y,       // Pixel coordinates
     72                    const poisConfig *config // Configuration
     73    );
     74
    6675/************************************************************************************************************/
    6776
    … …  
    7887
    7988// Print usage information
    80 //void help(void);
     89void help(void);
    8190
    8291// Parse the command line arguments
    83 poisConfig *restrict poisConfigAlloc(void);
    84 
    8592poisConfig *poisParseConfig(int argc,   // Number of command-line arguments
    8693                            char **argv // Command-line arguments
    … …  
    111118
    112119// Find stamps
    113 psArray *poisFindStamps(psArray *stampsIO, // Existing list of stamps, or NULL
     120psArray *poisFindStamps(psArray *stamps, // Existing list of stamps, or NULL
    114121                        const psImage *image, // Image for which to find stamps
    115122                        const psImage *mask, // Mask image
    … …  
    158165
    159166// Calculate deviations from the best fit for the stamps
    160 psVector *poisCalculateDeviations(psVector *deviationsIO, // Output array of deviations, or NULL
     167psVector *poisCalculateDeviations(psVector *deviations, // Output array of deviations, or NULL
    161168                                  psArray *stamps, // Array of stamps
    162169                                  psImage *refImage, // Reference image
    … …  
    173180
    174181// Reject stamps, based on the deviations
    175 psList *restrict poisRejectStamps(psArray *stamps,      // Array of stamps
    176                                   psImage *mask,        // Mask image
    177                                   const psVector *deviations, // Vector of deviations for the stamps
    178                                   const poisConfig *config // Configuration
     182bool poisRejectStamps(psArray *stamps,  // Array of stamps
     183                      psImage *mask,    // Mask image
     184                      const psVector *deviations, // Vector of deviations for the stamps
     185                      const poisConfig *config // Configuration
    179186    );
    180187
    … …  
    184191
    185192// Solve the matrix equation
    186 psVector *poisSolveEquation(psVector *solutionIO,       // Solution vector, or NULL
     193psVector *poisSolveEquation(psVector *solution, // Solution vector, or NULL
    187194                            const psArray *stamps, // Array of stamps
    188195                            const poisConfig *config
    … …  
    203210/************************************************************************************************************/
    204211
    205 psImage *poisImageSetVal(psImage *restrict IO,       // Input image
    206                          const psC64 val             // set to this value
    207     );
    208 
    209 psImage *poisImageSetValInMask(psImage *IO,                     // Image to update, or NULL
    210                                const psImage *restrict in,      // Input image
    211                                const psImage *mask,             // mask for image
    212                                const psC64 val,                 // set to this value
    213                                const unsigned long bits         // set pixels where mask & bits != 0
    214     );
     212// poisReadStamps.c
     213
     214// Read stamps from a file
     215psArray *poisReadStamps(psArray *stamps,// Stamps
     216                        FILE **stampFP, // File pointer for stamps file
     217                        const psImage *image, // Image for which to find stamps
     218                        const psImage *mask, // Mask image
     219                        const poisConfig *config // Configuration
     220    );
     221
     222/************************************************************************************************************/
     223
     224// poisCheckKernel.c
     225
     226// Check kernel for power at large radii
     227bool poisCheckKernel(const psVector *solution,// Kernel solution
     228                     const psArray *kernels, // Kernel basis functions
     229                     const poisConfig *config // Configuration
     230    );
     231
     232/************************************************************************************************************/
    215233
    216234#endif

  • trunk/pois/src/poisCalculateDeviations.c

    r3798 r5717  
    4646        int y = stamp->y;               // Stamp y coord
    4747        if (stamp->status == POIS_STAMP_USED) {
    48             psImage *refStamp = psImageSubset(refImage, x - xSize, y - ySize, x + xSize, y + ySize);
    49             psImage *inStamp = psImageSubset(inImage, x - xSize, y - ySize, x + xSize, y + ySize);
    50             psImage *maskStamp = psImageSubset(mask, x - xSize, y - ySize, x + xSize, y + ySize);
     48            psRegion stampRegion = {x - xSize, x + xSize, y - ySize, y + ySize};
     49            psImage *refStamp = psImageSubset(refImage, stampRegion);
     50            psImage *inStamp = psImageSubset(inImage, stampRegion);
     51            psImage *maskStamp = psImageSubset(mask, stampRegion);
    5152            psImage *convRefStamp = poisConvolveImage(refStamp, maskStamp, solution, kernelParams, config);
    5253            // Calculate chi^2
    5354            (void)psBinaryOp(subStamp, inStamp, "-", convRefStamp);
     55            (void)psBinaryOp(subStamp, subStamp, "*", subStamp);
    5456            (void)psBinaryOp(subStamp, subStamp, "/", inStamp);
    55             (void)psBinaryOp(subStamp, subStamp, "*", subStamp);
    56             psImage *subStampTrim = psImageSubset(subStamp, config->xKernel, config->yKernel,
    57                                                   config->xKernel + 2 * footprint,
    58                                                   config->yKernel + 2 * footprint);
    59             psImage *maskStampTrim = psImageSubset(maskStamp, config->xKernel, config->yKernel,
    60                                                    config->xKernel + 2 * footprint,
    61                                                    config->yKernel + 2 * footprint);
     57            psRegion stampTrim = { config->xKernel, config->xKernel + 2 * footprint, config->yKernel,
     58                                   config->yKernel + 2 * footprint };
     59            psImage *subStampTrim = psImageSubset(subStamp, stampTrim);
     60            psImage *maskStampTrim = psImageSubset(maskStamp, stampTrim);
    6261            // Copy image to workaround bug 305
    6362            psImage *tempImage = psImageCopy(NULL, subStampTrim, PS_TYPE_F32);
    … …  
    6968            psFree(tempMask);
    7069           
    71             deviations->data.F32[s] = stats->sampleMean * (float)footprint * (float)footprint * 4.0;
     70            deviations->data.F32[s] = sqrtf(stats->sampleMean / 2.0);
    7271            psTrace("pois.calculateDeviations", 5, "Deviation of stamp %d (%d,%d) is %f\n", s, x, y,
    7372                    deviations->data.F32[s]);
    … …  
    7776            snprintf(stampName, MAXCHAR, "stamp%d.fits", s);
    7877            psFits *stampFile = psFitsAlloc(stampName);
    79             if (!psFitsWriteImage(stampFile, NULL, subStampTrim, 0, NULL)) {
     78            if (!psFitsWriteImage(stampFile, NULL, subStampTrim, 0)) {
    8079                psErrorStackPrint(stderr, "Unable to write stamp: %s\n", stampName);
    8180            }
    … …  
    8382#endif
    8483
     84#if 0
    8585            psFree(convRefStamp);
    8686            psFree(maskStampTrim);
    … …  
    8989            psFree(refStamp);
    9090            psFree(inStamp);
     91#endif
    9192        }
    9293    }

  • trunk/pois/src/poisCalculateEquation.c

    r3799 r5717  
    77/* Calculate the matrix and vector for the matrix equation */
    88int poisCalculateEquation(psArray *stamps, // The stamps
    9                           const psImage *refImage, // The reference image
    10                           const psImage *inImage, // The input image
    11                           const psArray *kernelParams, // Kernel parameters
    12                           const poisConfig *config // Configuration
     9                          const psImage *refImage, // The reference image
     10                          const psImage *inImage, // The input image
     11                          const psArray *kernelParams, // Kernel parameters
     12                          const poisConfig *config // Configuration
    1313    )
    1414{
    … …  
    2323    // in the background levels.
    2424    int numSolveParams = kernelParams->n + 1;
    25     int bgIndex = kernelParams->n;      // Index in matrix for the background (it's the last one, after the
    26                                         // kernel parameters)
     25    int bgIndex = kernelParams->n;      // Index in matrix for the background (it's the last one, after the
     26                                        // kernel parameters)
    2727    float **reference = refImage->data.F32; // The reference pixels
    28     float **input = inImage->data.F32;  // The input pixels
    29     int nPix = 0;                       // Number of pixels used in calculation (for interest)
     28    float **input = inImage->data.F32;  // The input pixels
     29    int nPix = 0;                       // Number of pixels used in calculation (for interest)
    3030    float nx = 0.5 * (float)refImage->numCols; // Half-size of input image in x
    3131    float ny = 0.5 * (float)refImage->numRows; // Half-size of input image in y
    3232    int numKernelParams = kernelParams->n; // Number of kernel parameters
    3333
    34     psDPolynomial2D *poly = psDPolynomial2DAlloc(config->spatialOrder + 1, config->spatialOrder + 1,
    35                                                 PS_POLYNOMIAL_ORD); // Polynomial for evaluation
     34    psPolynomial2D *poly = psPolynomial2DAlloc(config->spatialOrder + 1, config->spatialOrder + 1,
     35                                              PS_POLYNOMIAL_ORD); // Polynomial for evaluation
    3636    for (int j = 0; j < config->spatialOrder + 1; j++) {
    37         for (int i = 0; i < config->spatialOrder + 1; i++) {
    38             poly->coeff[j][i] = 1.0;
    39             poly->mask[j][i] = 1;       // Mask all coefficients; unmask to evaluate
    40         }
     37        for (int i = 0; i < config->spatialOrder + 1; i++) {
     38            poly->coeff[j][i] = 1.0;
     39            poly->mask[j][i] = 1;       // Mask all coefficients; unmask to evaluate
     40        }
    4141    }
    4242
    4343    psImage *polyValues = psImageAlloc(config->spatialOrder + 1, config->spatialOrder + 1,
    44                                        PS_TYPE_F64); // Evaluations for each of the polynomial orders
     44                                       PS_TYPE_F64); // Evaluations for each of the polynomial orders
    4545
    4646    /* Iterate over the stamps */
    4747    for (int s = 0; s < stamps->n; s++) {
    48         poisStamp *stamp = stamps->data[s]; // The stamp
     48        poisStamp *stamp = stamps->data[s]; // The stamp
    4949
    50         if (stamp->status == POIS_STAMP_RECALC) {
    51             psTrace("pois.calculateEquation", 5, "Calculating for stamp %d: %d,%d\n", s, stamp->x, stamp->y);
     50        if (stamp->status == POIS_STAMP_RECALC) {
     51            psTrace("pois.calculateEquation", 5, "Calculating for stamp %d: %d,%d\n", s, stamp->x, stamp->y);
    5252
    53             // Initialise the matrix and vector if required
    54             if (! stamp->matrix) {
    55                 psTrace("pois.calculateEquation", 9, "Allocating matrix: %dx%d\n", numSolveParams,
    56                         numSolveParams);
    57                 psImage *matrix = psImageAlloc(numSolveParams, numSolveParams, PS_TYPE_F64);
    58                 for (int j = 0; j < numSolveParams; j++) {
    59                     for (int i = 0; i < numSolveParams; i++) {
    60                         matrix->data.F64[j][i] = 0.0;
    61                     }
    62                 }
    63                 stamp->matrix = matrix;
    64             }
    65             if (! stamp->vector) {
    66                 psTrace("pois.calculateEquation", 9, "Allocating vector: %d\n", numSolveParams);
    67                 psVector *vector = psVectorAlloc(numSolveParams, PS_TYPE_F64);
    68                 for (int i = 0; i < numSolveParams; i++) {
    69                     vector->data.F64[i] = 0.0;
    70                 }
    71                 stamp->vector = vector;
    72             }
    73            
    74             // Dereference, for convenience
    75             psImage *matrix = stamp->matrix;
    76             psVector *vector = stamp->vector;
    77            
    78             // Will assert on the type, since can't do much about it if it's wrong.
    79             assert(matrix->type.type == PS_TYPE_F64);
    80             assert(vector->type.type == PS_TYPE_F64);
    81            
    82             // Evaluate the polynomial for each order
    83             for (int j = 0; j < config->spatialOrder + 1; j++) {
    84                 for (int i = 0; i < config->spatialOrder + 1 - j; i++) {
    85                     poly->mask[j][i] = 0;
    86                     polyValues->data.F64[j][i] = psDPolynomial2DEval(poly, ((double)stamp->x - nx)/nx,
    87                                                                      ((double)stamp->y - ny)/ny);
    88                     poly->mask[j][i] = 1;
    89                 }
    90             }
    91            
    92             // Iterate over the pixels
    93             for (int y = stamp->y - config->footprint; y < stamp->y + config->footprint; y++) {
    94                 for (int x = stamp->x - config->footprint; x < stamp->x + config->footprint; x++) {
    95                     float invNoise2 = 1.0/reference[y][x]; // The inverse of the noise, squared.
    96                     nPix++;
    97                    
    98                     /* Iterate over the first convolution */
    99                     for (int k1 = 0; k1 < numKernelParams; k1++) {
    100                         poisKernelBasis *kernel1 = kernelParams->data[k1]; // Kernel basis #1
    101                        
    102                         int u1 = kernel1->u; // Offset in x
    103                         int v1 = kernel1->v; // Offset in y
    104                         int i1 = kernel1->xOrder; // Polynomial order in x
    105                         int j1 = kernel1->yOrder; // Polynomial order in y
    106                         // First convolution
    107                         float conv1 = polyValues->data.F64[j1][i1] * reference[y - v1][x - u1];
     53            // Initialise the matrix and vector if required
     54            if (! stamp->matrix) {
     55                psTrace("pois.calculateEquation", 9, "Allocating matrix: %dx%d\n", numSolveParams,
     56                        numSolveParams);
     57                stamp->matrix = psImageAlloc(numSolveParams, numSolveParams, PS_TYPE_F64);
     58            }
     59            if (! stamp->vector) {
     60                psTrace("pois.calculateEquation", 3, "Allocating vector: %d\n", numSolveParams);
     61                stamp->vector = psVectorAlloc(numSolveParams, PS_TYPE_F64);
     62            }
    10863
    109                         // Assuming that the first kernel component is 0 order in x and y, and 0 offset
    110                         if (k1 != 0) {
    111                             conv1 -= reference[y][x];
    112                         }
    113                        
    114                         /* Iterate over the second convolution */
    115                         for (int k2 = k1; k2 < numKernelParams; k2++) {
    116                             poisKernelBasis *kernel2 = kernelParams->data[k2]; // Kernel basis #2
    117                            
    118                             int u2 = kernel2->u; // Offset in x
    119                             int v2 = kernel2->v; // Offset in y
    120                             int i2 = kernel2->xOrder; // Polynomial order in x
    121                             int j2 = kernel2->yOrder; // Polynomial order in y
    122                             // Second convolution
    123                             float conv2 = polyValues->data.F64[j2][i2] * reference[y - v2][x - u2];
    124                            
    125                             // Assuming that the first kernel component is 0 order in x and y, and 0 offset
    126                             if (k2 != 0) {
    127                                 conv2 -= reference[y][x];
    128                             }
     64            // Dereference, for convenience
     65            psImage *matrix = stamp->matrix;
     66            psVector *vector = stamp->vector;
    12967
    130                             // Add into the matrix element
    131                             matrix->data.F64[k1][k2] += conv1 * // First convolution
    132                                 conv2 * // Second convolution
    133                                 invNoise2; // Divide by sigma^2, ~ poisson.
    134                         } // Second convolution
    135                        
    136                         // Add into the vector element
    137                         vector->data.F64[k1] += input[y][x] * // Input image, no convolution
    138                             conv1 *     // Convolved reference image
    139                             invNoise2; // Divide by sigma^2, ~ poisson.
    140                        
    141                         /* Background term */
    142                         matrix->data.F64[k1][bgIndex] += conv1 * // Convolved reference image
    143                             invNoise2; // Divide by sigma^2, ~ poisson.
    144                        
    145                     } // First convolution
    146                    
    147                     /* Background only terms */
    148                     matrix->data.F64[bgIndex][bgIndex] += invNoise2;
    149                     vector->data.F64[bgIndex] += input[y][x] * invNoise2;
    150                 }
    151             } // Iterating over pixels in stamp
     68            // Will assert on the type, since can't do much about it if it's wrong.
     69            assert(matrix->type.type == PS_TYPE_F64);
     70            assert(vector->type.type == PS_TYPE_F64);
    15271
    153             // Fill in other side of symmetric matrix
    154             for (int k1 = 0; k1 < kernelParams->n; k1++) {
    155                 for (int k2 = 0; k2 < k1; k2++) {
    156                     matrix->data.F64[k1][k2] = matrix->data.F64[k2][k1];
    157                 }
    158                 matrix->data.F64[bgIndex][k1] = matrix->data.F64[k1][bgIndex];
    159             }
     72            // Initialise the matrix and vector
     73            for (int j = 0; j < numSolveParams; j++) {
     74                for (int i = 0; i < numSolveParams; i++) {
     75                    matrix->data.F64[j][i] = 0.0;
     76                }
     77                vector->data.F64[j] = 0.0;
     78            }
    16079
    161             for (int k = 0; k < kernelParams->n; k++) {
    162                 poisKernelBasis *kernel = kernelParams->data[k]; // Kernel basis #2
    163                 int u = kernel->u; // Offset in x
    164                 int v = kernel->v; // Offset in y
    165                 matrix->data.F64[k][k] += config->penalty * (float)(u*u + v*v);
    166             }
     80            // Evaluate the polynomial for each order
     81            for (int j = 0; j < config->spatialOrder + 1; j++) {
     82                for (int i = 0; i < config->spatialOrder + 1 - j; i++) {
     83                    poly->mask[j][i] = 0;
     84                    polyValues->data.F64[j][i] = psPolynomial2DEval(poly, ((double)stamp->x - nx)/nx,
     85                                                                    ((double)stamp->y - ny)/ny);
     86                    poly->mask[j][i] = 1;
     87                }
     88            }
    16789
    168             stamp->status = POIS_STAMP_USED; // We've calculated it now, so don't need to do so again.
     90            // Iterate over the pixels
     91            for (int y = stamp->y - config->footprint; y < stamp->y + config->footprint; y++) {
     92                for (int x = stamp->x - config->footprint; x < stamp->x + config->footprint; x++) {
     93                    float invNoise2 = 1.0/reference[y][x]; // The inverse of the noise, squared.
     94                    nPix++;
    16995
    170         } // Legitimate stamps
     96                    /* Iterate over the first convolution */
     97                    for (int k1 = 0; k1 < numKernelParams; k1++) {
     98                        poisKernelBasis *kernel1 = kernelParams->data[k1]; // Kernel basis #1
     99
     100                        int u1 = kernel1->u; // Offset in x
     101                        int v1 = kernel1->v; // Offset in y
     102                        int i1 = kernel1->xOrder; // Polynomial order in x
     103                        int j1 = kernel1->yOrder; // Polynomial order in y
     104                        // First convolution
     105                        float conv1 = polyValues->data.F64[j1][i1] * reference[y - v1][x - u1];
     106
     107                        // Assuming that the first kernel component is 0 order in x and y, and 0 offset
     108                        if (k1 != 0) {
     109                            conv1 -= reference[y][x];
     110                        }
     111
     112                        /* Iterate over the second convolution */
     113                        for (int k2 = k1; k2 < numKernelParams; k2++) {
     114                            poisKernelBasis *kernel2 = kernelParams->data[k2]; // Kernel basis #2
     115
     116                            int u2 = kernel2->u; // Offset in x
     117                            int v2 = kernel2->v; // Offset in y
     118                            int i2 = kernel2->xOrder; // Polynomial order in x
     119                            int j2 = kernel2->yOrder; // Polynomial order in y
     120                            // Second convolution
     121                            float conv2 = polyValues->data.F64[j2][i2] * reference[y - v2][x - u2];
     122
     123                            // Assuming that the first kernel component is 0 order in x and y, and 0 offset
     124                            if (k2 != 0) {
     125                                conv2 -= reference[y][x];
     126                            }
     127
     128                            // Add into the matrix element
     129                            matrix->data.F64[k1][k2] += conv1 * // First convolution
     130                                conv2 * // Second convolution
     131                                invNoise2; // Divide by sigma^2, ~ poisson.
     132                        } // Second convolution
     133
     134                        // Add into the vector element
     135                        vector->data.F64[k1] += input[y][x] * // Input image, no convolution
     136                            conv1 *     // Convolved reference image
     137                            invNoise2; // Divide by sigma^2, ~ poisson.
     138
     139                        /* Background term */
     140                        matrix->data.F64[k1][bgIndex] += conv1 * // Convolved reference image
     141                            invNoise2; // Divide by sigma^2, ~ poisson.
     142
     143                    } // First convolution
     144
     145                    /* Background only terms */
     146                    matrix->data.F64[bgIndex][bgIndex] += invNoise2;
     147                    vector->data.F64[bgIndex] += input[y][x] * invNoise2;
     148                }
     149            } // Iterating over pixels in stamp
     150
     151            // Fill in other side of symmetric matrix
     152            for (int k1 = 0; k1 < kernelParams->n; k1++) {
     153                for (int k2 = 0; k2 < k1; k2++) {
     154                    matrix->data.F64[k1][k2] = matrix->data.F64[k2][k1];
     155                }
     156                matrix->data.F64[bgIndex][k1] = matrix->data.F64[k1][bgIndex];
     157            }
     158
     159            for (int k = 0; k < kernelParams->n; k++) {
     160                poisKernelBasis *kernel = kernelParams->data[k]; // Kernel basis #2
     161                int u = kernel->u; // Offset in x
     162                int v = kernel->v; // Offset in y
     163                matrix->data.F64[k][k] += config->penalty * (float)(u*u + v*v);
     164            }
     165
     166            stamp->status = POIS_STAMP_USED; // We've calculated it now, so don't need to do so again.
     167
     168        } // Legitimate stamps
    171169
    172170    } // Iterating over stamps

  • trunk/pois/src/poisConvolveImage.c

    r3801 r5717  
    77/* Convolve an image by the kernel */
    88psImage *poisConvolveImage(const psImage *input, // Input image, to be convolved
    9                            const psImage *mask, // Mask image
    10                            const psVector *solution, // The solution vector, containing the coefficients
    11                            const psArray *kernelParams, // The kernel parameters
    12                            const poisConfig *config // The configuration
     9                           const psImage *mask, // Mask image
     10                           const psVector *solution, // The solution vector, containing the coefficients
     11                           const psArray *kernelParams, // The kernel parameters
     12                           const poisConfig *config // The configuration
    1313    )
    1414{
    15     assert(solution->n == kernelParams->n + 1); // Extra parameter for background
     15    assert(solution->n == kernelParams->n + 1); // Extra parameter for background
    1616    assert(solution->type.type == PS_TYPE_F64);
    1717    assert(input->numCols == mask->numCols && input->numRows == mask->numRows);
    1818    assert(mask->type.type == PS_TYPE_U8);
    1919    assert(input->type.type == PS_TYPE_F32);
    20    
    21     int nx = input->numCols;            // The size of the image in x
    22     int ny = input->numRows;            // The size of the image in y
    23     int xKernel = config->xKernel;      // The half-size of the kernel in x
    24     int yKernel = config->yKernel;      // The half-size of the kernel in y
    25     float background = solution->data.F64[solution->n - 1]; // The difference in background
    26     int numParams = kernelParams->n;    // Number of kernel parameters
     20
     21    int nx = input->numCols;            // The size of the image in x
     22    int ny = input->numRows;            // The size of the image in y
     23    float nxHalf = 0.5 * (float)nx;     // Half the size of the image
     24    float nyHalf = 0.5 * (float)ny;     // Half the size of the image
     25    int xKernel = config->xKernel;      // The half-size of the kernel in x
     26    int yKernel = config->yKernel;      // The half-size of the kernel in y
     27    float background = solution->data.F64[solution->n-1]; // The difference in background
     28    int numParams = kernelParams->n;    // Number of kernel parameters
    2729
    2830    psImage *convolved = psImageAlloc(nx, ny, PS_TYPE_F32); // The convolved image, to be returned
    2931
    30     psDPolynomial2D *poly = psDPolynomial2DAlloc(config->spatialOrder + 1, config->spatialOrder + 1,
    31                                                  PS_POLYNOMIAL_ORD); // Polynomial
     32    psPolynomial2D *poly = psPolynomial2DAlloc(config->spatialOrder + 1, config->spatialOrder + 1,
     33                                               PS_POLYNOMIAL_ORD ); // Polynomial
    3234    for (int i = 0; i < config->spatialOrder + 1; i++) {
    33         for (int j = 0; j < config->spatialOrder + 1; j++) {
    34             poly->coeff[i][j] = 1.0;
    35             poly->mask[i][j] = 1;       // Mask all coefficients; unmask to evaluate
    36         }
     35        for (int j = 0; j < config->spatialOrder + 1; j++) {
     36            poly->coeff[i][j] = 1.0;
     37            poly->mask[i][j] = 1;       // Mask all coefficients; unmask to evaluate
     38        }
    3739    }
    3840
    3941    // Convolve only the middle part
    4042    for (int y = yKernel; y < ny - yKernel; y++) {
    41         for (int x = xKernel; x < nx - xKernel; x++) {
    42             if (! mask->data.U8[y][x] & (POIS_MASK_BAD | POIS_MASK_NEAR_BAD)) {
    43                 double conv = background; // Convolved value
     43        for (int x = xKernel; x < nx - xKernel; x++) {
     44            if (! mask->data.U8[y][x] & (POIS_MASK_BAD | POIS_MASK_NEAR_BAD)) {
     45                double conv = background; // Convolved value
    4446
    45 #if 0
    46                 float imageX = (x - nxHalf) / nxHalf; // Normalised position in x
    47                 float imageY = (y - nyHalf) / nyHalf; // Normalised position in y
    48 #endif
    49                
    50                 // Iterate over the kernel basis functions
    51                 for (int k = 0; k < numParams; k++) {
    52                     poisKernelBasis *kernel = kernelParams->data[k]; // The kernel basis function
    53                     int u = kernel->u;
    54                     int v = kernel->v;
    55                    
    56 #if 0
    57                     int xOrder = kernel->xOrder;
    58                     int yOrder = kernel->yOrder;
    59                     // Evaluate the polynomial
    60                     poly->mask[xOrder][yOrder] = 0;
    61                     double polyVal = psDPolynomial2DEval(poly, imageX, imageY);
    62                     poly->mask[xOrder][yOrder] = 1;
    63 #else
    64                     double polyVal = 1.0;
    65 #endif
     47                float imageX = (x - nxHalf) / nxHalf; // Normalised position in x
     48                float imageY = (y - nyHalf) / nyHalf; // Normalised position in y
    6649
    67                     if (k == 0) {
    68                         conv += solution->data.F64[k] * input->data.F32[y - v][x - u] * polyVal;
    69                     } else {
    70                         conv += solution->data.F64[k] *
    71                             (input->data.F32[y - v][x - u] * polyVal - input->data.F32[y][x]);
    72                     }
    73                 }
    74                 convolved->data.F32[y][x] = (float)conv;
    75             }
    76         }
    77            
    78         /* Pad the rest with zeros */
    79         for (int x = 0; x < xKernel; x++) {
    80             convolved->data.F32[y][x] = 0.0;
    81         }
    82         for (int x = nx - xKernel; x < nx; x++) {
    83             convolved->data.F32[y][x] = 0.0;
    84         }
     50                // Iterate over the kernel basis functions
     51                for (int k = 0; k < numParams; k++) {
     52                    poisKernelBasis *kernel = kernelParams->data[k]; // The kernel basis function
     53                    int u = kernel->u;
     54                    int v = kernel->v;
     55                    int xOrder = kernel->xOrder;
     56                    int yOrder = kernel->yOrder;
     57
     58                    double polyVal = 1.0;
     59                    if (xOrder != 0 || yOrder != 0) {
     60                        // Evaluate the polynomial
     61                        poly->mask[xOrder][yOrder] = 0;
     62                        polyVal = psPolynomial2DEval(poly, imageX, imageY);
     63                        poly->mask[xOrder][yOrder] = 1;
     64                    }
     65
     66                    if (k == 0) {
     67                        conv += solution->data.F64[k] * input->data.F32[y - v][x - u] * polyVal;
     68                    } else {
     69                        conv += solution->data.F64[k] *
     70                            (input->data.F32[y - v][x - u] * polyVal - input->data.F32[y][x]);
     71                    }
     72                }
     73                convolved->data.F32[y][x] = (float)conv;
     74            }
     75        }
     76
     77        /* Pad the rest with zeros */
     78        for (int x = 0; x < xKernel; x++) {
     79            convolved->data.F32[y][x] = 0.0;
     80        }
     81        for (int x = nx - xKernel; x < nx; x++) {
     82            convolved->data.F32[y][x] = 0.0;
     83        }
    8584    }
    8685
    8786    for (int y = 0; y < yKernel; y++) {
    88         for (int x = 0; x < nx; x++) {
    89             convolved->data.F32[y][x] = 0.0;
    90         }
     87        for (int x = 0; x < nx; x++) {
     88            convolved->data.F32[y][x] = 0.0;
     89        }
    9190    }
    9291    for (int y = ny - yKernel; y < ny; y++) {
    93         for (int x = 0; x < nx; x++) {
    94             convolved->data.F32[y][x] = 0.0;
    95         }
     92        for (int x = 0; x < nx; x++) {
     93            convolved->data.F32[y][x] = 0.0;
     94        }
    9695    }
    9796

  • trunk/pois/src/poisExtractKernel.c

    r3792 r5717  
    3030    }
    3131
    32     psDPolynomial2D *poly = psDPolynomial2DAlloc(config->spatialOrder + 1, config->spatialOrder + 1,
    33                                                 PS_POLYNOMIAL_ORD); // Polynomial for evaluation
     32    psPolynomial2D *poly = psPolynomial2DAlloc(config->spatialOrder + 1, config->spatialOrder + 1,
     33                                              PS_POLYNOMIAL_ORD); // Polynomial for evaluation
    3434    for (int j = 0; j < config->spatialOrder + 1; j++) {
    3535        for (int i = 0; i < config->spatialOrder + 1; i++) {
    … …  
    5151        // Evaluate polynomial
    5252        poly->mask[xOrder][yOrder] = 0;
    53         double evaluation = solution->data.F64[k] * psDPolynomial2DEval(poly, x, y);
     53        double evaluation = solution->data.F64[k] * psPolynomial2DEval(poly, x, y);
    5454        poly->mask[xOrder][yOrder] = 1;
    5555

  • trunk/pois/src/poisFindStamps.c

    r3813 r5717  
    77
    88psArray *poisFindStamps(psArray *stamps, // Existing list of stamps, or NULL
    9                         const psImage *image, // Image for which to find stamps
    10                         const psImage *mask, // Mask image
    11                         const poisConfig *config // Configuration values
     9                        const psImage *image, // Image for which to find stamps
     10                        const psImage *mask, // Mask image
     11                        const poisConfig *config // Configuration values
    1212    )
    1313{
    14     int nx = config->xImage;            // Size of image in x
    15     int ny = config->yImage;            // Size of image in y
     14    int nx = config->xImage;            // Size of image in x
     15    int ny = config->yImage;            // Size of image in y
    1616
    1717    assert(image->numCols == nx && image->numRows == ny);
    … …  
    2121    assert(!stamps || stamps->n == config->nsx * config->nsy);
    2222
    23     float threshold = config->threshold;// Stamp threshold
    24     int footprint = config->footprint;  // Footprint for stamps
     23    int footprint = config->footprint;  // Footprint for stamps
    2524    int borderx = footprint + config->xKernel; // Border around image in x
    2625    int bordery = footprint + config->yKernel; // Border around image in y
    27     int nsx = config->nsx;              // Number of stamps in x
    28     int nsy = config->nsy;              // Number of stamps in y
     26    int nsx = config->nsx;              // Number of stamps in x
     27    int nsy = config->nsy;              // Number of stamps in y
    2928
    30     psF32 **pixels = image->data.F32;   // The image pixels
    31     psU8 **maskpix = mask->data.U8; // The mask pixels
     29    psF32 **pixels = image->data.F32;   // The image pixels
    3230
    3331    if (stamps == NULL) {
    34         stamps = psArrayAlloc(nsx * nsy);
    35                                         // Initialise; should be done by psArrayAlloc
    36         for (int i = 0; i < nsx * nsy; i++) {
    37             stamps->data[i] = NULL;
    38         }
     32        stamps = psArrayAlloc(nsx * nsy);
     33        // Initialise
     34        for (int i = 0; i < nsx * nsy; i++) {
     35            stamps->data[i] = poisStampAlloc();
     36        }
    3937    }
    4038
    … …  
    4240    int num = 0;
    4341    for (int i = 0; i < nsx; i++) {
    44         for (int j = 0; j < nsy; j++) {
    45             poisStamp *stamp = stamps->data[num]; // The stamp
    46             if (stamp == NULL) {
    47                 stamp = stamps->data[num] = poisStampAlloc();
    48             }
     42        for (int j = 0; j < nsy; j++) {
     43            poisStamp *stamp = stamps->data[num]; // The stamp
    4944
    50             // Only find a new stamp if we need to
    51             if (stamp->status == POIS_STAMP_RESET) {
    52                 // Find maximum non-masked value in the image section, but don't include a footprint around
    53                 // the edge
    54                 float max = - INFINITY;         // Negative infinity
    55                 int bestx = 0, besty = 0;       // Position of maximum
    56                
    57                 for (int y = bordery + j * (ny - 2.0 * bordery) / nsy;
    58                      y < bordery + (j + 1) * (ny - 2.0 * bordery) / nsy; y++) {
    59                     for (int x = borderx + i * (nx - 2.0 * borderx) / nsx;
    60                          x < borderx + (i + 1) * (nx - 2.0 * borderx) / nsx; x++) {
    61                         if (pixels[y][x] > max && pixels[y][x] >= threshold) {
    62                             // Check the footprint
    63                             int ok = 1;
    64                             for (int v = y - footprint; v <= y + footprint && ok; v++) {
    65                                 for (int u = x - footprint; u <= x + footprint && ok; u++) {
    66                                     if (maskpix[v][u] &
    67                                         (POIS_MASK_BAD | POIS_MASK_NEAR_BAD | POIS_MASK_STAMP)) {
    68                                         ok = 0;
    69                                     }
    70                                 }
    71                             }
    72                             if (ok) {
    73                                 max = pixels[y][x];
    74                                 bestx = x;
    75                                 besty = y;
    76                             }
    77                         } // Done checking the candidate pixel
    78                     }
    79                 } // Done iterating over this image section
    80                
    81                 // Store the stamp
    82                 if (max != - INFINITY) {
    83                     stamp->x = bestx;
    84                     stamp->y = besty;
    85                     stamp->status = POIS_STAMP_RECALC;
    86                 } else {
    87                     // No stamp in this section
    88                     stamp->status = POIS_STAMP_BAD;
    89                 }
    90             } // Finding a new stamp
     45            // Only find a new stamp if we need to
     46            if ((stamp->x == 0 && stamp->y == 0) || stamp->status == POIS_STAMP_RESET) {
     47                // Find maximum non-masked value in the image section, but don't include a footprint around
     48                // the edge
     49                float max = - INFINITY;         // Negative infinity
     50                int bestx = 0, besty = 0;       // Position of maximum
    9151
    92             psTrace("pois.findStamps", 7, "Stamp %d: %d,%d\n", num, stamp->x, stamp->y);
    93             num++;
    94         }
     52                for (int y = bordery + j * (ny - 2.0 * bordery) / nsy;
     53                     y < bordery + (j + 1) * (ny - 2.0 * bordery) / nsy; y++) {
     54                    for (int x = borderx + i * (nx - 2.0 * borderx) / nsx;
     55                         x < borderx + (i + 1) * (nx - 2.0 * borderx) / nsx; x++) {
     56                        if (pixels[y][x] > max && poisCheckStamp(image, mask, x, y, config)) {
     57                            max = pixels[y][x];
     58                            bestx = x;
     59                            besty = y;
     60                        }
     61                    }
     62                } // Done iterating over this image section
     63
     64                // Store the stamp
     65                if (max != - INFINITY) {
     66                    stamp->x = bestx;
     67                    stamp->y = besty;
     68                    stamp->status = POIS_STAMP_RECALC;
     69                } else {
     70                    // No stamp in this section
     71                    stamp->status = POIS_STAMP_BAD;
     72                }
     73            } // Finding a new stamp
     74
     75            psTrace("pois.findStamps", 7, "Stamp %d: %d,%d\n", num, stamp->x, stamp->y);
     76            num++;
     77        }
    9578    } // Done iterating over sections
    9679
    9780    return stamps;
    98 }                                   
     81}

  • trunk/pois/src/poisKernelBasisFunctions.c

    r3800 r5717  
    3636            array->data[num] = kernel;
    3737
    38             psTrace("pois.kernelBasisFunctions", 8, "--> %d: (%d,%d) x^%d y^%d\n", i, kernel->u,
     38            psTrace("pois.kernelBasisFunctions", 3, "--> %d: (%d,%d) x^%d y^%d\n", i, kernel->u,
    3939                    kernel->v, kernel->xOrder, kernel->yOrder);
    4040            num++;
    … …  
    6464                        array->data[num] = kernel;
    6565
    66                         psTrace("pois.kernelBasisFunctions", 8, "--> %d: (%d,%d) x^%d y^%d\n", num,
     66                        psTrace("pois.kernelBasisFunctions", 3, "--> %d: (%d,%d) x^%d y^%d\n", num,
    6767                                kernel->u, kernel->v, kernel->xOrder, kernel->yOrder);
    6868                        num++;

  • trunk/pois/src/poisMakeMask.c

    r3792 r5717  
    4141            // Mask all around something that's saturated in the reference, since it gets convolved
    4242            // But only mask single pixel in the input image, since it doesn't get convolved
    43             if (refImage->data.F32[y][x] >= config->refSat || refImage->data.F32[y][x] <= config->bad ||
     43            if (refImage->data.F32[y][x] >= config->refSat || refImage->data.F32[y][x] <= config->refBad ||
    4444                isnan(refImage->data.F32[y][x])) {
    4545                mask->data.U8[y][x] |= POIS_MASK_BAD;
    … …  
    5353                    }
    5454                }
    55             } else if (inImage->data.F32[y][x] >= config->inSat || inImage->data.F32[y][x] <= config->bad ||
     55            } else if (inImage->data.F32[y][x] >= config->inSat || inImage->data.F32[y][x] <= config->inBad ||
    5656                isnan(inImage->data.F32[y][x])) {
    5757                mask->data.U8[y][x] |= POIS_MASK_BAD;

  • trunk/pois/src/poisParseConfig.c

    r3805 r5717  
    44#include "pois.h"
    55
    6 static void help(void)
     6void help(void)
    77{
    88    fprintf (stderr, "POIS: Pan-STARRS (or Pricey's) Optimal Image Subtraction\n"
    9              "Usage: pois [-h] [-v] [-k X Y] [-t T] [-s S1 S2] [-b BAD] [-B BG] [-o N] [-n NSTAMPS] [-i ITER] [-r REJ] REF IN OUT\n"
     9             "Usage: pois [-h] [-v] [-k X Y] [-t T] [-s S1 S2] [-b B1 B2] [-B BG] [-o N] [-n NSTAMPS] [-S STAMPS] [-i ITER] [-r REJ] [-m] REF IN OUT\n"
    1010             "where\n"
    1111             "\t-h         Help (this info)\n"
    … …  
    1414             "\t-t T       Threshold for stamps on reference image (0)\n"
    1515             "\t-s S1 S2   Saturation level for reference (S1) and input (S2) (50000)\n"
    16              "\t-b BAD     Bad level for both images (0)\n"
     16             "\t-b B1 B2   Bad level for both images (0)\n"
    1717             "\t-B BG      Level to add to background to get above zero (0)\n"
    1818             "\t-o ORDER   Order of spatial variations in the kernel (0)\n"
    1919             "\t-n NSTAMPS Number of stamps in each dimension (5)\n"
     20             "\t-S STAMPS  File from which to read stamps\n"
     21             "\t-R         Reverse the sense of the subtraction (REF - IN)\n"
    2022             "\t-i ITER    Number of rejection iterations (10)\n"
    2123             "\t-r REJ     Rejection level in standard deviations (2.5)\n"
     24             "\t-m         Output mask frame (OUT.mask)\n"
    2225             "\tREF        Reference image\n"
    2326             "\tIN         Input image (to be matched to REF)\n"
    … …  
    2629}
    2730
    28 poisConfig *restrict poisConfigAlloc(void)
     31poisConfig *poisParseConfig(int argc, // Number of command-line arguments
     32                            char **argv // Command-line arguments
     33    )
    2934{
    30     poisConfig *config = (poisConfig *) psAlloc(sizeof(poisConfig));
     35    poisConfig *config;                 // Configuration values
     36
     37    /* Variables for getopt */
     38    int opt;   /* Option, from getopt */
     39    extern char *optarg;   /* Argument accompanying switch */
     40    extern int optind;   /* getopt variables */
     41
     42    /* Initialise configuration */
     43    config = (poisConfig *) psAlloc(sizeof(poisConfig));
    3144    config->verbose = 0;
    3245    config->refFile = NULL;
    … …  
    3952    config->refSat = 60000.0;
    4053    config->inSat = 60000.0;
    41     config->bad = 0.0;
     54    config->refBad = 0.0;
     55    config->inBad = 0.0;
    4256    config->background = 0.0;
    4357    config->spatialOrder = 0;
    … …  
    4761    config->sigmaRej = 2.5;
    4862    config->penalty = 0.0;
    49 
    50     return config;
    51 }
    52 
    53 
    54 poisConfig *poisParseConfig(int argc, // Number of command-line arguments
    55                             char **argv // Command-line arguments
    56     )
    57 {
    58     poisConfig *config  = poisConfigAlloc(); // Configuration values
    59 
    60     /* Variables for getopt */
    61     int opt;   /* Option, from getopt */
    62     extern char *optarg;   /* Argument accompanying switch */
    63     extern int optind;   /* getopt variables */
    64 
    65     /* Initialise configuration */
     63    config->stampFile = NULL;
     64    config->reverse = false;
     65    config->mask = false;
    6666
    6767    /* Parse command-line arguments using getopt */
    68     while ((opt = getopt(argc, argv, "hvk:f:q:t:s:o:b:B:n:i:r:p:")) != -1) {
     68    while ((opt = getopt(argc, argv, "hvRmk:f:q:t:s:o:b:B:n:i:r:p:S:")) != -1) {
    6969        switch (opt) {
    7070          case 'h':
    … …  
    115115            break;
    116116          case 'b':
    117             if (sscanf(optarg, "%f", &config->bad) != 1) {
    118                 help();
     117            if ((sscanf(argv[optind-1], "%f", &config->refBad) != 1) ||
     118                (sscanf(argv[optind++], "%f", &config->inBad) != 1)) {
     119                /*
     120                  Note: incrementing optind, so I can read more than
     121                  one parameter.
     122                */
     123                help();
    119124                exit(EXIT_FAILURE);
    120             }
     125            }
    121126            break;
    122127          case 'B':
    … …  
    151156            }
    152157            break;
     158          case 'S':
     159            config->stampFile = argv[optind-1];
     160            config->threshold = -INFINITY;
     161            break;
     162          case 'R':
     163            config->reverse = true;
     164            break;
     165          case 'm':
     166            config->mask = true;
     167            break;
    153168          default:
    154169            help();

  • trunk/pois/src/poisReadStamps.c

    r4644 r5717  
    33#include "pois.h"
    44
    5 #define MAXLINE 128                     // Maximum length of a line in the stamps file
     5#define MAXLINE 128                     // Maximum length of a line in the stamps file
    66
    77
    88psArray *poisReadStamps(psArray *stamps,// Stamps
    9                         FILE **stampFP, // File pointer for stamps file
    10                         const psImage *image, // Image for which to find stamps
    11                         const psImage *mask, // Mask image
    12                         const poisConfig *config // Configuration
     9                        FILE **stampFP, // File pointer for stamps file
     10                        const psImage *image, // Image for which to find stamps
     11                        const psImage *mask, // Mask image
     12                        const poisConfig *config // Configuration
    1313    )
    1414{
    15     FILE *fp = *stampFP;                // File pointer
    16 
    17     int num = 0;                        // Number of stamps
     15    FILE *fp = *stampFP;                // File pointer
    1816
    1917    if (stamps == NULL) {
    20         stamps = psArrayAlloc(config->nsx);
    21         // Initialise
    22         for (int i = 0; i < config->nsx; i++) {
    23             stamps->data[i] = poisStampAlloc();
    24         }
     18        stamps = psArrayAlloc(config->nsx);
     19        // Initialise
     20        for (int i = 0; i < config->nsx; i++) {
     21            stamps->data[i] = poisStampAlloc();
     22        }
    2523    }
    2624
    2725    for (int i = 0; i < stamps->n; i++) {
    28         poisStamp *stamp = stamps->data[i];
     26        poisStamp *stamp = stamps->data[i];
    2927
    30         if (stamp->status == POIS_STAMP_RESET) {
    31             if (fp) {
    32                 float x, y;             // Position of stamp
    33                 int xPixel, yPixel;     // Integer position
    34                 int result = 0;         // Result of fscanf
    35                 bool badPixels = false; // Result of check on stamp footprint
    36                 char line[MAXLINE];     // The line from the stamps file
    37                 char *lineBack;         // The line back from fgets
    38                 do {
    39                     lineBack = fgets(line, MAXLINE, fp);
    40                     result = sscanf(line, "%f %f", &x, &y);
     28        if (stamp->status == POIS_STAMP_RESET) {
     29            if (fp) {
     30                float x, y;             // Position of stamp
     31                int xPixel, yPixel;     // Integer position
     32                int result = 0;         // Result of fscanf
     33                char line[MAXLINE];     // The line from the stamps file
     34                char *lineBack;         // The line back from fgets
     35                do {
     36                    lineBack = fgets(line, MAXLINE, fp);
     37                    result = sscanf(line, "%f %f", &x, &y);
    4138
    42                     // Add 0.5 for rounding, and subtract 1.0 for unit-offset to zero-offset
    43                     xPixel = (int)(x - 0.5);
    44                     yPixel = (int)(y - 0.5);
    45                 } while (result == 2 && lineBack && ! poisCheckStamp(image, mask, xPixel, yPixel, config));
    46                 if (result == 2 && lineBack) {
    47                     stamp->x = xPixel;
    48                     stamp->y = yPixel;
    49                     stamp->status = POIS_STAMP_RECALC;
    50                     psTrace("pois.readStamps", 7, "Stamp %d: %d,%d\n", i, stamp->x, stamp->y);
    51                 } else {
    52                     // We've run out of stamps in the file --- probably got EOF
    53                     stamp->status = POIS_STAMP_BAD;
    54                 }
    55             } else {
    56                 // Stamp file isn't open!
    57                 stamp->status = POIS_STAMP_BAD;
    58             }
    59         }
     39                    // Add 0.5 for rounding, and subtract 1.0 for unit-offset to zero-offset
     40                    xPixel = (int)(x - 0.5);
     41                    yPixel = (int)(y - 0.5);
     42                } while (result == 2 && lineBack && ! poisCheckStamp(image, mask, xPixel, yPixel, config));
     43                if (result == 2 && lineBack) {
     44                    stamp->x = xPixel;
     45                    stamp->y = yPixel;
     46                    stamp->status = POIS_STAMP_RECALC;
     47                    psTrace("pois.readStamps", 7, "Stamp %d: %d,%d\n", i, stamp->x, stamp->y);
     48                } else {
     49                    // We've run out of stamps in the file --- probably got EOF
     50                    stamp->status = POIS_STAMP_BAD;
     51                }
     52            } else {
     53                // Stamp file isn't open!
     54                stamp->status = POIS_STAMP_BAD;
     55            }
     56        }
    6057    }
    6158

  • trunk/pois/src/poisRejectStamps.c

    r3875 r5717  
    66#define SQUARE(x) ((x)*(x))
    77
    8 psList *restrict poisRejectStamps(psArray *stamps, // Array of stamps
    9                                   psImage *mask,        // Mask image
    10                                   const psVector *deviations, // Vector of deviations for the stamps
    11                                   const poisConfig *config // Configuration
     8bool poisRejectStamps(psArray *stamps, // Array of stamps
     9                      psImage *mask,    // Mask image
     10                      const psVector *deviations, // Vector of deviations for the stamps
     11                      const poisConfig *config // Configuration
    1212    )
    1313{
    … …  
    2121    assert(deviations->type.type == PS_TYPE_F32);
    2222
    23     psVector *devMask = psVectorAlloc(stamps->n, PS_TYPE_U8); // Mask for statistics
     23    bool masked = false;                // Have we masked any stamps?
    2424
    2525    // Don't want the standard deviation, but rather the deviation from zero
    26     double meanDev = 0.0;
     26    double sum = 0.0;
    2727    int numDev = 0;
    2828    for (int i = 0; i < deviations->n; i++) {
    29         poisStamp *stamp = stamps->data[i];     // The stamp
    30         // Only interested in the stamps that we're actually using
    31         if (stamp != NULL && stamp->status == POIS_STAMP_USED) {
    32             meanDev += SQUARE(deviations->data.F32[i]);
    33             numDev++;
    34         }
     29        poisStamp *stamp = stamps->data[i];     // The stamp
     30        // Only interested in the stamps that we're actually using
     31        if (stamp->status == POIS_STAMP_USED) {
     32            sum += SQUARE(deviations->data.F32[i]);
     33            numDev++;
     34        }
    3535    }
    36     float rmsDev = sqrt(meanDev/numDev);
    37     float limit = rmsDev * config->sigmaRej;
    38     psTrace("pois.rejectStamps", 2, "RMS deviation: %f\n", rmsDev);
     36    float meanDev = sqrtf(sum / (float)numDev);
     37    float limit = meanDev * config->sigmaRej;
     38    psTrace("pois.rejectStamps", 2, "Mean RMS deviation: %f\n", meanDev);
    3939    psTrace("pois.rejectStamps", 2, "Rejection limit: %f\n", limit);
    4040
    4141    // Reject stamps
    42     psList *rejected = NULL;            // rejected stamps
    4342    for (int s = 0; s < deviations->n; s++) {
    44         poisStamp *stamp = stamps->data[s];
    45         if (stamp == NULL) {
    46             continue;
    47         }
    48        
    49         if (stamp->status == POIS_STAMP_USED && fabsf(deviations->data.F32[s]) > limit) {
    50             psTrace("pois.rejectStamps", 1, "Rejecting stamp %d (%d,%d): %f\n", s, stamp->x, stamp->y,
    51                     deviations->data.F32[s]);
    52            
    53             // Mask out the stamp in the image so you don't find it again
    54             for (int y = stamp->y - config->footprint; y < stamp->y + config->footprint; y++) {
    55                 for (int x = stamp->x - config->footprint; x < stamp->x + config->footprint; x++) {
    56                     mask->data.U8[y][x] |= POIS_MASK_STAMP;
    57                 }
    58             }
    59            
    60             // Mark stamp for replacement
    61             stamps->data[s] = NULL;
    62            
    63             if (rejected == NULL) {
    64                 rejected = psListAlloc(stamp);
    65             } else {
    66                 bool success = psListAdd(rejected, PS_LIST_TAIL, stamp);
    67                 assert(success);
    68             }
    69             psFree(stamp);
    70         } // Bad stamps
     43        poisStamp *stamp = stamps->data[s];
     44        if (stamp->status == POIS_STAMP_USED &&
     45            (fabsf(deviations->data.F32[s]) > limit || ! isfinite(deviations->data.F32[s]))) {
     46            masked = true;
     47            psTrace("pois.rejectStamps", 1, "Rejecting stamp %d (%d,%d): %f\n", s, stamp->x, stamp->y,
     48                    deviations->data.F32[s]);
     49
     50            // Mask out the stamp in the image so you don't find it again
     51            for (int y = stamp->y - config->footprint; y < stamp->y + config->footprint; y++) {
     52                for (int x = stamp->x - config->footprint; x < stamp->x + config->footprint; x++) {
     53                    mask->data.U8[y][x] |= POIS_MASK_STAMP;
     54                }
     55            }
     56
     57            // Set stamp for replacement
     58            stamp->x = 0;
     59            stamp->y = 0;
     60            stamp->status = POIS_STAMP_RESET;
     61
     62        } // Bad stamps
    7163    } // Iterating over stamps
    7264
    73     psFree(devMask);
    74 
    75     return rejected;
     65    return masked;
    7666}

  • trunk/pois/src/poisSolveEquation.c

    r3792 r5717  
    3939
    4040        if (stamp->status == POIS_STAMP_USED) {
    41 
    42 #ifdef TESTING
    43             printf("Matrix %d:\n", s);
    44             for (int i = 0; i < matrix->numCols; i++) {
    45                 for (int j = 0; j < matrix->numRows; j++) {
    46                     printf("%f ", matrix->data.F64[i][j]);
    47                 }
    48                 printf("\n");
    49             }
    50 #endif
    51 
    5241            (void)psBinaryOp(matrix, matrix, "+", stampMatrix);
    5342            (void)psBinaryOp(vector, vector, "+", stampVector);
    … …  
    6150    psTrace("pois.solveEquation", 1, "Background difference is %f\n", solution->data.F64[solution->n - 1]);
    6251
     52#if 0
     53    printf("Matrix:\n");
     54    for (int i = 0; i < matrix->numCols; i++) {
     55        for (int j = 0; j < matrix->numRows; j++) {
     56            printf("%f ", matrix->data.F64[i][j]);
     57        }
     58        printf("\n");
     59    }
     60#endif
     61
    6362    // Clean up
    6463    psFree(luMatrix);

  • trunk/pois/src/poisStamp.c

    r3836 r5717  
    1111    stamp->matrix = NULL;
    1212    stamp->vector = NULL;
    13     psMemSetDeallocator(stamp, (psFreeFcn)poisStampFree);
     13    psMemSetDeallocator(stamp, (psFreeFunc)poisStampFree);
    1414
    1515    return stamp;
    … …  
    1818void poisStampFree(poisStamp *stamp)
    1919{
    20     psFree(stamp->matrix);
    21     psFree(stamp->vector);
     20    if (stamp->matrix) {
     21        psFree(stamp->matrix);
     22    }
     23    if (stamp->vector) {
     24        psFree(stamp->vector);
     25    }
     26    psFree(stamp);
    2227}
     28
     29
     30// Return true if the stamp is OK
     31bool poisCheckStamp(const psImage *image, // Image to check for threshold
     32                    const psImage *mask, // Mask to check for bad pixels
     33                    int x, int y,       // Pixel coordinates
     34                    const poisConfig *config // Configuration
     35    )
     36{
     37    if (x < config->footprint + config->xKernel ||
     38        y < config->footprint + config->yKernel ||
     39        x > config->xImage - config->footprint - config->xKernel ||
     40        y > config->yImage - config->footprint - config->yKernel ||
     41        image->data.F32[y][x] < config->threshold) {
     42        psTrace("pois.checkStamp", 9, "Rejecting stamp at %d,%d (border/threshold)\n", x, y);
     43        return false;
     44    }
     45
     46    // Check the footprint
     47    bool ok = true;                     // Is the footprint OK?
     48    int footprint = config->footprint;  // Footprint size
     49    for (int v = y - footprint; v <= y + footprint && ok; v++) {
     50        for (int u = x - footprint; u <= x + footprint && ok; u++) {
     51            if (mask->data.U8[v][u] &
     52                (POIS_MASK_BAD | POIS_MASK_NEAR_BAD | POIS_MASK_STAMP)) {
     53                psTrace("pois.checkStamp", 9, "Rejecting stamp at %d,%d (bad footprint)\n", x, y);
     54                ok = false;
     55            }
     56        }
     57    }
     58   
     59    return ok;
     60}
     61
Note: See TracChangeset for help on using the changeset viewer.