IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 16607


Ignore:
Timestamp:
Feb 22, 2008, 9:50:56 AM (18 years ago)
Author:
Paul Price
Message:

Previous merge didn't have the expected result, probably because of the revert-and-branch I did earlier. I think this is the expected result.

Location:
trunk/psModules/src/imcombine
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/imcombine/pmPSFEnvelope.c

    r16479 r16607  
    3131
    3232
    33 #define TESTING                         // Enable test output
     33//#define TESTING                         // Enable test output
    3434#define PEAK_FLUX 1.0e3                 // Peak flux for each source
    3535#define SKY_VALUE 0.0e0                 // Sky value for fake image
    … …  
    9494            fake->peak = pmPeakAlloc(xFake - dx, yFake - dy, PEAK_FLUX, PM_PEAK_LONE);
    9595            fake->type = PM_SOURCE_TYPE_STAR;
    96 
    97 #ifdef TESTING
    98             // Required to generate model image
    9996            fake->psfMag = -2.5 * log10(PEAK_FLUX);
    100 #endif
    10197
    10298            psTrace("psModules.imcombine", 5, "Source %d: %.2f,%.2f\n",

  • trunk/psModules/src/imcombine/pmStack.c

    r16604 r16607  
    88 *  @author GLG, MHPCC
    99 *
     10 *  @version $Revision: 1.19 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2008-02-22 19:50:56 $
    1012 *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
    1113 *
    … …  
    510512        for (int j = 0; j < pixels->n; j++) {
    511513            int x = pixels->data[j].x, y = pixels->data[j].y; // Coordinates of interest
     514            if (x < 0 || x >= numCols || y < 0 || y >= numRows) {
     515                psWarning("Bad pixel coordinate: %d,%d --- ignored.", x, y);
     516                continue;
     517            }
    512518            psArray *columns = map->data[y]; // The columns for that row
    513519            psVector *images = columns->data[x]; // The images for that column

  • trunk/psModules/src/imcombine/pmStackReject.c

    r16604 r16607  
    1111#define PIXEL_LIST_BUFFER 100           // Number of pixels to add to list at a time
    1212
    13 psPixels *pmStackReject(const psPixels *in, float threshold, const psArray *regions,
    14                         const psArray *solutions, const pmSubtractionKernels *kernels)
     13psPixels *pmStackReject(const psPixels *in, const psRegion *valid, float threshold,
     14                        const psArray *subRegions, const psArray *kernels)
    1515{
    1616    PS_ASSERT_PIXELS_NON_NULL(in, NULL);
    1717    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(threshold, 0.0, NULL);
    1818    PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(threshold, 1.0, NULL);
    19     PS_ASSERT_ARRAY_NON_NULL(regions, NULL);
    20     PS_ASSERT_ARRAY_NON_NULL(solutions, NULL);
    21     PS_ASSERT_ARRAYS_SIZE_EQUAL(regions, solutions, NULL);
    22     PS_ASSERT_PTR_NON_NULL(kernels, NULL);
     19    PS_ASSERT_ARRAY_NON_NULL(subRegions, NULL);
     20    PS_ASSERT_ARRAY_NON_NULL(kernels, NULL);
     21    PS_ASSERT_ARRAYS_SIZE_EQUAL(subRegions, kernels, NULL);
    2322
    2423    // Get the original image size
    25     int numRegions = regions->n;        // Number of regions
     24    int numRegions = subRegions->n;        // Number of regions
    2625    int numCols = 0, numRows = 0;       // Size of original image
    2726    int minCols = INT_MAX, minRows = INT_MAX; // Minimum coordinate for image
    2827    int size = 0;                       // Size of kernel
    2928    for (int i = 0; i < numRegions; i++) {
    30         psRegion *region = regions->data[i]; // Region of interest
    31         if (region->x0 < minCols) {
    32             minCols = region->x0;
     29        psRegion *subRegion = subRegions->data[i]; // Region of interest
     30        if (subRegion->x0 < minCols) {
     31            minCols = subRegion->x0;
    3332        }
    34         if (region->y0 < minRows) {
    35             minRows = region->y0;
     33        if (subRegion->y0 < minRows) {
     34            minRows = subRegion->y0;
    3635        }
    37         if (region->x1 > numCols) {
    38             numCols = region->x1;
     36        if (subRegion->x1 > numCols) {
     37            numCols = subRegion->x1;
    3938        }
    40         if (region->y1 > numRows) {
    41             numRows = region->y1;
     39        if (subRegion->y1 > numRows) {
     40            numRows = subRegion->y1;
     41        }
     42
     43        pmSubtractionKernels *kernel = kernels->data[i]; // Kernel of interest
     44        if (size == 0) {
     45            size = kernel->size;
     46        } else if (kernel->size != size) {
     47            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Kernel sizes are not identical: %d vs %d",
     48                    size, kernel->size);
     49            return NULL;
    4250        }
    4351    }
    44     if (minCols != 0 || minRows != 0) {
    45         psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    46                 "Some error with image regions --- minimum coordinate is not 0,0");
    47         return NULL;
     52
     53    // Adjust the size for the size of the subimage
     54    if (valid) {
     55        minCols = PS_MAX(valid->x0, minCols);
     56        minRows = PS_MAX(valid->y0, minRows);
     57        numCols = PS_MIN(valid->x1, numCols);
     58        numRows = PS_MIN(valid->y1, numRows);
    4859    }
    4960    psTrace("psModules.imcombine", 1, "Rejecting [%d:%d,%d:%d]\n", minCols, numCols, minRows, numRows);
    5061
    51     psImage *mask = psPixelsToMask(NULL, in, psRegionSet(0, numCols, 0, numRows), 0x01); // Mask image
     62    psImage *mask = psPixelsToMask(NULL, in, psRegionSet(minCols, numCols - 1, minRows, numRows - 1),
     63                                   0x01); // Mask
    5264    psImage *image = psImageCopy(NULL, mask, PS_TYPE_F32); // Floating-point version, so we can convolve
    5365    psFree(mask);
    … …  
    5769    pmReadout *inRO = pmReadoutAlloc(NULL); // Readout with input image
    5870    inRO->image = image;
     71    inRO->col0 = minCols;
     72    inRO->row0 = minRows;
    5973    for (int i = 0; i < numRegions; i++) {
    6074        psRegion *region = subRegions->data[i]; // Region of interest
    … …  
    7690
    7791        // Image of the kernel at the centre of the region
    78         float xNorm = (region->x0 + 0.5 * (region->x1 - region->x0) - numCols/2.0) / (float)numCols;
    79         float yNorm = (region->y0 + 0.5 * (region->y1 - region->y0) - numRows/2.0) / (float)numRows;
    80         psImage *kernel = pmSubtractionKernelImage(kernels, xNorm, yNorm, false);
    81         if (!kernel) {
     92        float xNorm = (region->x0 + 0.5 * (region->x1 - region->x0) - kernel->numCols/2.0) /
     93            (float)kernel->numCols;
     94        float yNorm = (region->y0 + 0.5 * (region->y1 - region->y0) - kernel->numRows/2.0) /
     95            (float)kernel->numRows;
     96        psImage *image = pmSubtractionKernelImage(kernel, xNorm, yNorm, false);
     97        if (!image) {
    8298            psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");
    8399            psFree(convRO);
    … …  
    86102        }
    87103        float sum = 0.0;
    88         for (int y = 0; y < kernel->numRows; y++) {
    89             for (int x = 0; x < kernel->numCols; x++) {
    90                 sum += kernel->data.F32[y][x];
     104        for (int y = 0; y < image->numRows; y++) {
     105            for (int x = 0; x < image->numCols; x++) {
     106                sum += image->data.F32[y][x];
    91107            }
    92108        }
    93         psFree(kernel);
     109        psFree(image);
    94110
    95         psImage *subConv = psImageSubset(convRO->image, *region); // Sub-image of convolved image
    96         psBinaryOp(subConv, subConv, "*", psScalarAlloc(1.0 / sum, PS_TYPE_F32));
     111        // Range for normalisation
     112        int yMin = PS_MAX(minRows, region->y0) - inRO->row0;
     113        int yMax = PS_MIN(numRows - 1, region->y1) - inRO->row0;
     114        int xMin = PS_MAX(minCols, region->x0) - inRO->col0;
     115        int xMax = PS_MIN(numCols - 1, region->x1) - inRO->col0;
     116        psTrace("psModules.imcombine", 2, "Normalising convolved mask image by %f over %d:%d,%d:%d\n",
     117                sum, xMin, xMax, yMin, yMax);
     118        for (int y = yMin; y <= yMax; y++) {
     119            for (int x = xMin; x <= xMax; x++) {
     120                convRO->image->data.F32[y][x] /= sum;
     121            }
     122        }
    97123    }
    98124    psFree(inRO);
    … …  
    124150        for (int y = yMin; y <= yMax; y++) {
    125151            for (int x = xMin; x <= xMax; x++) {
     152                assert(x < mask->numCols && y < mask->numRows);
    126153                mask->data.PS_TYPE_MASK_DATA[y][x] = 0xff;
    127154            }

  • trunk/psModules/src/imcombine/pmStackReject.h

    r16479 r16607  
    1010/// We apply a matched filter to the corresponding mask image, and threshold to find the original pixels
    1111psPixels *pmStackReject(const psPixels *in, ///< List of pixels in the convolved image
     12                        const psRegion *valid, ///< Valid region to consider
    1213                        float threshold, ///< Threshold on convolved image, 0..1
    1314                        const psArray *regions, ///< Array of image regions for image
    14                         const psArray *solutions, ///< Array of solution vectors for image
    15                         const pmSubtractionKernels *kernels ///< Kernel parameters
     15                        const psArray *kernels ///< Array of kernel parameters for each region
    1616    );
    1717

  • trunk/psModules/src/imcombine/pmSubtraction.c

    r16604 r16607  
    44 *  @author GLG, MHPCC
    55 *
     6 *  @version $Revision: 1.82 $ $Name: not supported by cvs2svn $
     7 *  @date $Date: 2008-02-22 19:50:56 $
     8 *
    69 *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
    710 *
    … …  
    4750
    4851    // Take the square of the normal kernel
    49     double sumNormal = 0.0, sumVariance = 0.0; // Sum of the normal and variance kernels
    5052    for (int v = yMin; v <= yMax; v++) {
    5153        for (int u = xMin; u <= xMax; u++) {
    52             float value = normalKernel->kernel[v][u]; // Value of interest
    53             float value2 = PS_SQR(value); // Value squared
    54             sumNormal += value;
    55             sumVariance += value2;
    56             out->kernel[v][u] = value2;
    57         }
    58     }
    59 
    60     // Normalise so that the sum of the variance kernel is the square of the sum of the normal kernel
    61     // This is required to keep the relative scaling between the image and the weight map
    62     psBinaryOp(out->image, out->image, "*", psScalarAlloc(PS_SQR(sumNormal) / sumVariance, PS_TYPE_F32));
     54            out->kernel[v][u] = PS_SQR(normalKernel->kernel[v][u]);
     55        }
     56    }
    6357
    6458    return out;
    … …  
    687681        PS_ASSERT_IMAGE_TYPE(out1->weight, PS_TYPE_F32, false);
    688682    }
    689     if (region) {
    690         if (psRegionIsNaN(*region)) {
    691             psString string = psRegionToString(*region);
    692             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) contains NAN values", string);
    693             psFree(string);
    694             return false;
    695         }
    696         if (region->x0 < 0 || region->x1 > ro1->image->numCols ||
    697             region->y0 < 0 || region->y1 > ro1->image->numRows) {
    698             psString string = psRegionToString(*region);
    699             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) does not fit in image (%dx%d)",
    700                     string, ro1->image->numCols, ro1->image->numRows);
    701             psFree(string);
    702             return false;
    703         }
     683    if (region && psRegionIsNaN(*region)) {
     684        psString string = psRegionToString(*region);
     685        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) contains NAN values", string);
     686        psFree(string);
     687        return false;
    704688    }
    705689
    … …  
    776760        yMax = PS_MIN(region->y1, yMax);
    777761    }
     762
     763    // Size to use when calculating normalised coordinates (different from actual size when convolving
     764    // subimage)
     765    int xNormSize = (kernels->numCols > 0 ? kernels->numCols : numCols);
     766    int yNormSize = (kernels->numRows > 0 ? kernels->numRows : numRows);
    778767
    779768    psMaskType maskSource;              // Mask these pixels when convolving
    … …  
    803792    for (int j = yMin; j < yMax; j += fullSize) {
    804793        int ySubMax = PS_MIN(j + fullSize, yMax); // Range for subregion of interest
    805         float yNorm = 2.0 * (float)(j + size + 1 - numRows/2.0) / (float)numRows; // Normalised coordinate
     794        float yNorm = 2.0 * (float)(j + y0 + size + 1 - yNormSize/2.0) /
     795            (float)yNormSize; // Normalised coordinate
    806796        for (int i = xMin; i < xMax; i += fullSize) {
    807797            int xSubMax = PS_MIN(i + fullSize, xMax); // Range for subregion of interest
    808             float xNorm = 2.0 * (float)(i + size + 1 - numCols/2.0) / (float)numCols; // Normalised coordinate
     798            float xNorm = 2.0 * (float)(i + x0 + size + 1 - xNormSize/2.0) /
     799                (float)xNormSize; // Normalised coordinate
    809800
    810801            // Only generate polynomial values every kernel footprint, since we have already assumed

  • trunk/psModules/src/imcombine/pmSubtraction.h

    r16604 r16607  
    66 * @author GLG, MHPCC
    77 *
    8  * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
    9  * @date $Date: 2008-02-22 19:24:42 $
     8 * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
     9 * @date $Date: 2008-02-22 19:50:56 $
    1010 * Copyright 2004-207 Institute for Astronomy, University of Hawaii
    1111 */

  • trunk/psModules/src/imcombine/pmSubtractionEquation.c

    r16479 r16607  
    1414
    1515
    16 #define TESTING
     16//#define TESTING
    1717
    1818//////////////////////////////////////////////////////////////////////////////////////////////////////////////

  • trunk/psModules/src/imcombine/pmSubtractionKernels.c

    r16479 r16607  
    180180    kernels->bgOrder = 0;
    181181    kernels->mode = mode;
     182    kernels->numCols = 0;
     183    kernels->numRows = 0;
    182184    kernels->solution1 = NULL;
    183185    kernels->solution2 = NULL;

  • trunk/psModules/src/imcombine/pmSubtractionKernels.h

    r16479 r16607  
    3939    int bgOrder;                        ///< The order for the background fitting
    4040    pmSubtractionMode mode;             ///< Mode for subtraction
     41    int numCols, numRows;               ///< Size of image (for normalisation), or zero to use image provided
    4142    psVector *solution1, *solution2;    ///< Solution for the PSF matching
    4243} pmSubtractionKernels;

  • trunk/psModules/src/imcombine/pmSubtractionMatch.c

    r16604 r16607  
    306306                                                       inner, binning, ringsOrder, mode);
    307307            }
    308             psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL",
    309                              PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels);
     308
     309            // Add analysis metadata
     310            {
     311                psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_KERNEL,
     312                                 PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels);
     313                if (conv2) {
     314                    psMetadataAddPtr(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_KERNEL,
     315                                     PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels);
     316                }
     317                psMetadataAddS32(conv1->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE,
     318                                 PS_META_DUPLICATE_OK, "Subtraction kernels", mode);
     319                if (conv2) {
     320                    psMetadataAddS32(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE,
     321                                     PS_META_DUPLICATE_OK, "Subtraction kernels", mode);
     322                }
     323                psRegion *subRegion;
     324                if (region) {
     325                    subRegion = psMemIncrRefCounter(region);
     326                } else {
     327                    subRegion = psRegionAlloc(0, numCols, 0, numRows);
     328                }
     329                psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION,
     330                                  PS_DATA_REGION | PS_META_DUPLICATE_OK,
     331                                 "Region over which subtraction was performed", subRegion);
     332                if (conv2) {
     333                    psMetadataAddPtr(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION,
     334                                     PS_DATA_REGION | PS_META_DUPLICATE_OK,
     335                                     "Region over which subtraction was performed", subRegion);
     336                }
     337                psFree(subRegion);
     338            }
    310339
    311340            memCheck("kernels");
    … …  
    442471            kernels = NULL;
    443472
    444 #if 0
    445             // Put the solution on the metadata
    446             {
    447                 psString comment = NULL; // Comment for metadata
    448                 psStringAppend(&comment, "Subtraction solution for region %s", regionString);
    449                 psMetadataAddVector(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.SOLUTION",
    450                                     PS_META_DUPLICATE_OK, comment, solution);
    451                 psFree(comment);
    452                 if (region) {
    453                     psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.REGION",
    454                                      PS_META_DUPLICATE_OK | PS_DATA_REGION, comment, region);
    455                 } else {
    456                     region = psRegionAlloc(0, numCols, 0, numRows);
    457                     psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.REGION",
    458                                      PS_META_DUPLICATE_OK | PS_DATA_REGION, comment, region);
    459                     psFree(region);
    460                     region = NULL;
    461                 }
    462             }
    463 #endif
    464 
    465473            // There is data in the readout now
    466474            conv1->data_exists = true;

  • trunk/psModules/src/imcombine/pmSubtractionMatch.h

    r16479 r16607  
    88#include <pmSubtractionKernels.h>
    99#include <pmSubtractionStamps.h>
     10
     11#define PM_SUBTRACTION_ANALYSIS_KERNEL "SUBTRACTION.KERNEL" // Name of kernel in the analysis metadata
     12#define PM_SUBTRACTION_ANALYSIS_MODE "SUBTRACTION.MODE" // Name of subtraction mode in the analysis metadata
     13#define PM_SUBTRACTION_ANALYSIS_REGION "SUBTRACTION.REGION" // Name of subtraction region in the analysis MD
     14
    1015
    1116/// Match two images
Note: See TracChangeset for help on using the changeset viewer.