IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 4417


Ignore:
Timestamp:
Jun 28, 2005, 1:23:55 PM (21 years ago)
Author:
gusciora
Message:

Significant modifications.

Location:
trunk/psModules/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/pmImageSubtract.c

    r4322 r4417  
    11/** @file  ImageSubtract.c
    2 *
    3 *  This file will ...
    4 *
    5 *  @author Paul Price, IfA (original prototype)
    6 *  @author GLG, MHPCC
    7 *
    8 *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
    9 *  @date $Date: 2005-06-20 23:21:05 $
    10 *
    11 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    12 *
    13 *   XXX: sync with iFa on this:
    14 *   The (x, y) (row, col) issue is becoming a problem.  In this file, and I
    15 *   think, the rest of psLib and psModules, the following conventions are used:
    16 *
    17 * 1) x will correspond to the column, and y will correspond to the row.
    18 * 2) When used in function prototypes, the column (and hence x) appears
    19 *    first.
    20 * 3) When used to index 2-D arrays, obviously, the row (and hence, y)
    21 *    appears first. (2 and 3 are the source of confusion).
    22 * 4) When (u, v) are used in certain structures.
    23 *  u corresponds to x
    24 *  v corresponds to y
    25 * 5) When element (a, b) is casually referred to in comments, or
    26 *    documentation it is unclear where a is the row, or the column.
    27 * 6) A convention on loop index variables (i, j) would be convenient.
    28 *    Currently, sometimes i corresponds to the column (x),
    29 *    usually it corresponds to the row (y).
    30 *
    31 *  XXX: The following variables are used an interpreted this way:
    32 * kernelSize: Means that the actual kernel is a square (1 + 2 * kernelSize) per side.
    33 * border: When accessing an image, a swath of pixels this wide is ignored.
    34 * footprint: When accessing a stample, a square of pixels, footprint pixels per side,
    35 *  are looked at.  We must ensure that (footprint+kernelSize) pixels exist
    36 *  around the center.
    37 */
     2 *
     3 *  This file will contain code which creates a set of kernel basis
     4 *  functions, solves for their solution, and applies them to an image.
     5 *
     6 *  @author Paul Price, IfA (original prototype)
     7 *  @author GLG, MHPCC
     8 *
     9 *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2005-06-28 23:23:55 $
     11 *
     12 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     13 *
     14 *   XXX: sync with IfA on this:
     15 *   The (x, y) (row, col) issue is becoming a problem.  In this file, and I
     16 *   think, the rest of psLib and psModules, the following conventions are used:
     17 *
     18 * 1) x will correspond to the column, and y will correspond to the row.
     19 * 2) When used in function prototypes, the column (and hence x) appears
     20 *    first.
     21 * 3) When used to index 2-D arrays, obviously, the row (and hence, y)
     22 *    appears first. (2 and 3 are the source of confusion).
     23 * 4) When (u, v) are used in certain structures.
     24 *  u corresponds to x
     25 *  v corresponds to y
     26 * 5) When element (a, b) is casually referred to in comments, or
     27 *    documentation it is unclear where a is the row, or the column.
     28 * 6) A convention on loop index variables (i, j) would be convenient.
     29 *    Currently, sometimes i corresponds to the column (x),
     30 *    usually it corresponds to the row (y).
     31 *
     32 *  XXX: The following variables are used an interpreted this way:
     33 * kernelSize: Means that the actual kernel is a square (1 + 2 * kernelSize) per side.
     34 * border: When accessing an image, a swath of pixels this wide is ignored.
     35 * footprint: When accessing a stample, a square of pixels, footprint pixels per side,
     36 *  are looked at.  We must ensure that (footprint+kernelSize) pixels exist
     37 *  around the center.
     38 */
    3839
    3940#include<stdio.h>
    … …  
    4344#include "pmImageSubtract.h"
    4445
    45 // XXX: Put this is psConstants.h
    46 #define PS_ASSERT_INT_LESS_THAN(VAR1, VAR2, RVAL) \
    47 if (!(VAR1 < VAR2)) { \
    48     psError(PS_ERR_UNKNOWN, true, \
    49             "Error: %s is not less than %s (%d, %d)", #VAR1, #VAR2, VAR1, VAR2); \
    50     return(RVAL); \
    51 }
    52 
    5346/*******************************************************************************
    5447    Private alloc/free functions.
    55  
    56 XXX: It's not clear if the SubtractionKernels alloc/free functions are doable.
     48XXX: It's not clear if the SubtractionKernels alloc/free functions are feasable.
    5749 ******************************************************************************/
    5850void p_pmStampFree(pmStamp *stamp)
    … …  
    9082}
    9183
    92 psSubtractionKernels *p_pmSubtractionKernelsAlloc(psS32 numBasisFunctions,
     84psSubtractionKernels *p_pmSubtractionKernelsAlloc(int numBasisFunctions,
    9385        pmSubtractionKernelsType type)
    9486{
    … …  
    10698        int spatialOrder)
    10799{
     100    psTrace("ImageSubtract.pmSubtractionKernelsAllocPOIS", 1,
     101            "Calling pmSubtractionKernelsAllocPOIS(%d, %d)\n", size, spatialOrder);
    108102    PS_ASSERT_INT_POSITIVE(size, NULL);
    109103    PS_ASSERT_INT_POSITIVE(spatialOrder, NULL);
    … …  
    152146
    153147    //
    154     // Iterate over (u,v).  Generate a set of kernels foreach (u, v).
     148    // Iterate over (u,v).  Generate a set of kernels for each (u, v).
    155149    //
    156150    for (psS32 v = -yKernelHalfSize; v <= yKernelHalfSize; v++) {
    … …  
    176170    }
    177171
     172    psTrace("ImageSubtract.pmSubtractionKernelsAllocPOIS", 1,
     173            "Exiting pmSubtractionKernelsAllocPOIS(%d, %d)\n", size, spatialOrder);
    178174    return(tmp);
    179175}
    … …  
    181177/*******************************************************************************
    182178XXX: Get the types correct (u, v, xOrder, yOrder).
    183  
    184 XXX: For the various spatial orders, verify that the preCalc data does not
    185 vary with the (x^i * y^j) since that comes from the image coordinates when the
    186 kernel is convolved with it.
    187179 
    188180XXX: Code review this.
    … …  
    193185        int spatialOrder)
    194186{
     187    PS_ASSERT_VECTOR_NON_NULL(sigmas, NULL);
     188    PS_ASSERT_VECTOR_NON_NULL(orders, NULL);
     189    psTrace("ImageSubtract.pmSubtractionKernelsAllocISIS", 1,
     190            "Calling pmSubtractionKernelsAllocISIS(%d, %d, %d, %d)\n",
     191            sigmas->n, orders->n, size, spatialOrder);
    195192    PS_ASSERT_INT_POSITIVE(size, NULL);
    196193    PS_ASSERT_INT_POSITIVE(spatialOrder, NULL);
    197     PS_ASSERT_VECTOR_NON_NULL(sigmas, NULL);
    198194    PS_ASSERT_VECTOR_TYPE(sigmas, PS_TYPE_F32, NULL);
    199     PS_ASSERT_VECTOR_NON_NULL(orders, NULL);
    200195    PS_ASSERT_VECTOR_TYPE(orders, PS_TYPE_S32, NULL);
    201196    //
    202     // Calculate the number of basis functions (nBF), and others.
     197    // Calculate the number of basis functions (nBF).
    203198    //
    204199    psS32 numSigmas = sigmas->n;
    … …  
    283278                                                         (2.0 * PS_SQR(sigmas->data.F32[sigPtr]));
    284279                                        currPreCalc->data.F32[v+size][u+size] =
    285                                             exp(exponent) * pow(uScaled, orderXTerm) * pow(vScaled, orderYTerm);
     280                                            exp(-exponent) *
     281                                            pow(uScaled, orderXTerm) *
     282                                            pow(vScaled, orderYTerm);
    286283                                    }
    287284                                }
    … …  
    298295    }
    299296
     297    psTrace("ImageSubtract.pmSubtractionKernelsAllocISIS", 1,
     298            "Exiting pmSubtractionKernelsAllocISIS(%d, %d, %d, %d)\n",
     299            sigmas->n, orders->n, size, spatialOrder);
    300300    return(tmp);
    301301}
    … …  
    310310 
    311311XXX: Do we need to care about case where yNum/xNum does not evenly divide the
    312 nnumber of rows/columns in the mage?
     312nnumber of rows/columns in the image?
    313313 ******************************************************************************/
    314314psArray *pmSubtractionFindStamps(psArray *stamps,        ///< Output stamps, or NULL
    … …  
    322322                                )
    323323{
     324    psTrace("ImageSubtract.pmSubtractionFindStamps", 1,
     325            "Calling pmSubtractionFindStamps(%d, %f, %d, %d, %d)\n",
     326            maskVal, threshold, xNum, yNum, border);
    324327    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
    325328    PS_ASSERT_IMAGE_NON_EMPTY(image, NULL);
    … …  
    360363    //
    361364    // XXX: Must handle cases where image size is not an even multiple of xNum or yNum
     365    // they are currently ignored.
    362366    //
    363367    psS32 num = 0;
    … …  
    440444        }
    441445    }
     446    psTrace("ImageSubtract.pmSubtractionFindStamps", 1,
     447            "Exiting pmSubtractionFindStamps(%d, %f, %d, %d, %d)\n",
     448            maskVal, threshold, xNum, yNum, border);
    442449    return(stamps);
    443450}
    … …  
    455462                                psF32 y)
    456463{
     464    psTrace("ImageSubtract.GenSpatialOrder", 5,
     465            "Calling GenSpatialOrder(%d, %f, %f)\n", spatialOrder, x, y);
     466
    457467    psImage *polyValues = psImageAlloc(spatialOrder+1, spatialOrder+1, PS_TYPE_F64);
    458468
    … …  
    468478    }
    469479
     480    psTrace("ImageSubtract.GenSpatialOrder", 5,
     481            "Exiting GenSpatialOrder(%d, %f, %f)\n", spatialOrder, x, y);
     482
    470483    return(polyValues);
    471484}
    … …  
    473486
    474487/*******************************************************************************
    475 GeneralKernelConvolve(input, kernels, kernelID, col, row): This routine
     488IsisKernelConvolve(input, kernels, kernelID, col, row): This routine
    476489convolves a single kernel basis function with a pixel in an image.
    477  
    478 XXX: merge this code with the other convolution code in this file.
    479  
    480 XXX: Must we ensure that the pixels being accessed are inside the image?  Is
    481 there any garantee that the kernelSize, footprint stuff will be set correctly
    482 to ensure no seg faults.
    483  ******************************************************************************/
    484 static psF32 GeneralKernelConvolve(const psImage *input,
    485                                    const psSubtractionKernels *kernels,
    486                                    psS32 kernelID,
    487                                    psS32 col,
    488                                    psS32 row)
    489 {
     490  ******************************************************************************/
     491static psF32 IsisKernelConvolve(const psImage *input,
     492                                const psSubtractionKernels *kernels,
     493                                psS32 kernelID,
     494                                psS32 col,
     495                                psS32 row)
     496{
     497
     498
     499    psTrace("ImageSubtract.IsisKernelConvolve", 5,
     500            "Calling IsisKernelConvolve(%d, %d, %d)\n", kernelID, col, row);
    490501    psS32 spatialOrder = kernels->p_spatialOrder;
    491502    psS32 kernelSize = kernels->p_size;
    … …  
    496507    psF32 imageX = (((psF32) col) - numColsHalf) / numColsHalf; // Normalised position
    497508    psF32 imageY = (((psF32) row) - numRowsHalf) / numRowsHalf; // Normalised position
     509
    498510    psImage *polyValues = GenSpatialOrder(spatialOrder, imageX, imageY);
    499     psF32 polyVal = polyValues->data.F32[yOrder][xOrder];
     511
     512    psF64 polyVal = polyValues->data.F64[yOrder][xOrder];
    500513
    501514    psImage *preCalc = (psImage *) kernels->preCalc->data[kernelID];
     515
     516    // XXX: Are the following asserts really necessary?
     517    PS_ASSERT_INT_LARGER_THAN_OR_EQUAL(row-kernelSize, 0, NAN);
     518    PS_ASSERT_INT_LESS_THAN(row+kernelSize, input->numRows, NAN);
     519    PS_ASSERT_INT_LARGER_THAN_OR_EQUAL(col-kernelSize, 0, NAN);
     520    PS_ASSERT_INT_LESS_THAN(col+kernelSize, input->numCols, NAN);
    502521    psF32 conv = 0.0;
    503522    for (psS32 yy = -kernelSize ; yy < kernelSize ; yy++) {
    … …  
    510529    psFree(polyValues);
    511530
     531    psTrace("ImageSubtract.IsisKernelConvolve", 5,
     532            "Exiting IsisKernelConvolve(%d, %d, %d)\n", kernelID, col, row);
    512533    return(conv);
    513534}
    514535
    515 
    516 
    517 
    518536/*******************************************************************************
    519 XXX: Should we assert that the footprint is equal to the kernel size, or does
    520 they have nothing in common.
     537ConvolvePixelPois(input, mask, badStampMaskVal, solution, kernels, col, row):
     538 
     539This routine takes a single pixel in the psImage input and convolves it with
     540the set of kernel basis functions and their appropriate weights in solution.
     541It returns the value of the convolved pixel.
     542 
     543XXX: Static structure for polyValues?
     544 ******************************************************************************/
     545static psF32 ConvolvePixelPois(const psImage *input,
     546                               const psImage *mask,
     547                               psU32 badStampMaskVal,
     548                               const psVector *solution,
     549                               const psSubtractionKernels *kernels,
     550                               psS32 col,
     551                               psS32 row)
     552{
     553    psTrace("ImageSubtract.ConvolvePixelPois", 5,
     554            "Calling ConvolvePixelPois(%d, %d)\n", col, row);
     555    psS32 nBF = kernels->u->n;
     556    psF32 numColsHalf = 0.5 * (psF32) input->numCols;
     557    psF32 numRowsHalf = 0.5 * (psF32) input->numRows;
     558    psF32 background = solution->data.F64[solution->n-1];
     559    psS32 spatialOrder = kernels->p_spatialOrder;
     560    psF32 conv = background; // Initial convolved value
     561
     562    if ((mask == NULL) || !(mask->data.U8[row][col] & badStampMaskVal)) {
     563        psF32 imageX = (((psF32) col) - numColsHalf) / numColsHalf; // Normalised position
     564        psF32 imageY = (((psF32) row) - numRowsHalf) / numRowsHalf; // Normalised position
     565        psImage *polyValues = GenSpatialOrder(spatialOrder, imageX, imageY);
     566
     567        // Iterate over the kernel basis functions
     568        for (psS32 k = 0; k < nBF; k++) {
     569            psS32 u = (psS32) kernels->u->data.F32[k];
     570            psS32 v = (psS32) kernels->v->data.F32[k];
     571
     572            // XXX: What's the story with this?
     573            #if 0
     574
     575            psS32 xOrder = (psS32) kernels->xOrder->data.F32[k];
     576            psS32 yOrder = (psS32) kernels->yOrder->data.F32[k];
     577            psF64 polyVal = polyValues->data.F64[yOrder][xOrder];
     578            #else
     579
     580            psF32 polyVal = 1.0;
     581            #endif
     582
     583            // XXX: Why this?
     584            if (k == 0) {
     585                conv += solution->data.F64[k] * input->data.F32[row - v][col - u] * polyVal;
     586            } else {
     587                conv += solution->data.F64[k] *
     588                        (input->data.F32[row - v][col - u] * polyVal - input->data.F32[row][col]);
     589            }
     590        }
     591        psFree(polyValues);
     592    }
     593
     594    psTrace("ImageSubtract.ConvolvePixelPois", 5,
     595            "Exiting ConvolvePixelPois(%d, %d)\n", col, row);
     596    return(conv);
     597}
     598
     599
     600
     601/*******************************************************************************
     602ConvolvePixelIsis(input, mask, badStampMaskVal, solution, kernels, col, row):
     603 
     604This routine takes a single pixel in the psImage input and convolves it with
     605the set of kernel basis functions and their appropriate weights in solution.
     606It returns the value of the convolved pixel.
     607 
     608XXX: Static structure for polyValues?
     609 ******************************************************************************/
     610static psF32 ConvolvePixelIsis(const psImage *input,
     611                               const psImage *mask,
     612                               psU32 badStampMaskVal,
     613                               const psVector *solution,
     614                               const psSubtractionKernels *kernels,
     615                               psS32 col,
     616                               psS32 row)
     617{
     618    psTrace("ImageSubtract.ConvolvePixelIsis", 5,
     619            "Calling ConvolvePixelIsis(%d, %d)\n", col, row);
     620    psF32 background = solution->data.F64[solution->n-1];
     621    psF32 conv = background; // Initial convolved value
     622
     623    if ((mask == NULL) || !(mask->data.U8[row][col] & badStampMaskVal)) {
     624        // Iterate over the kernel basis functions
     625        for (psS32 k = 0; k < kernels->u->n; k++) {
     626            conv += IsisKernelConvolve(input, kernels, k, col, row);
     627        }
     628    }
     629
     630    psTrace("ImageSubtract.ConvolvePixelIsis", 5,
     631            "Exiting ConvolvePixelIsis(%d, %d)\n", col, row);
     632    return(conv);
     633}
     634
     635/*******************************************************************************
     636ConvolveImage(input, mask, badStampMaskVal, solution, kernels): convolves an
     637arbitrary image with either an ISIS or POIS set of kernel basis functions.
     638 ******************************************************************************/
     639static psImage *ConvolveImage(const psImage *input,
     640                              const psImage *mask,
     641                              psU32 badStampMaskVal,
     642                              const psVector *solution,
     643                              const psSubtractionKernels *kernels)
     644{
     645    psTrace("ImageSubtract.ConvolveImage", 5, "Calling ConvolveImage()\n");
     646    PS_ASSERT_IMAGE_NON_NULL(input, NULL);
     647    PS_ASSERT_IMAGE_NON_EMPTY(input, NULL);
     648    PS_ASSERT_IMAGE_TYPE(input, PS_TYPE_F32, NULL);
     649    if (mask != NULL) {
     650        PS_ASSERT_IMAGES_SIZE_EQUAL(input, mask, NULL);
     651        PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_U8, NULL);
     652    }
     653    PS_ASSERT_VECTOR_NON_NULL(solution, NULL);
     654    PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, NULL);
     655    PS_ASSERT_PTR_NON_NULL(kernels, NULL);
     656    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, NULL);
     657    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, NULL);
     658    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, NULL);
     659    if (kernels->preCalc != NULL) {
     660        PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, NULL);
     661    } else {
     662        if (kernels->type != PM_SUBTRACTION_KERNEL_POIS) {
     663            psError(PS_ERR_BAD_PARAMETER_NULL, true,
     664                    "Unallowable operation: kernels->preCalc == NULL and kernels->type != PM_SUBTRACTION_KERNEL_POIS.\n");
     665            return(NULL);
     666        }
     667    }
     668    psS32 nBF = kernels->u->n;
     669    PS_ASSERT_VECTOR_SIZE(solution, nBF+1, NULL);
     670
     671    psS32 numCols = input->numCols;
     672    psS32 numRows = input->numRows;
     673    psS32 kernelSize = kernels->p_size;
     674
     675    psImage *convolved = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     676
     677    for (psS32 y = kernelSize; y < numRows - kernelSize; y++) {
     678        for (psS32 x = kernelSize; x < numCols - kernelSize; x++) {
     679            if (kernels->type == PM_SUBTRACTION_KERNEL_POIS) {
     680                convolved->data.F32[y][x] = ConvolvePixelPois(input, mask, badStampMaskVal,
     681                                            solution, kernels, x, y);
     682            } else if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) {
     683                convolved->data.F32[y][x] = ConvolvePixelIsis(input, mask, badStampMaskVal,
     684                                            solution, kernels, x, y);
     685            } else {
     686                psLogMsg(__func__, PS_LOG_WARN, "WARNING: unknown kernel type.  Returning NULL\n");
     687                return(NULL);
     688            }
     689        }
     690    }
     691
     692    //
     693    // Pad the rest of the convolved image with 0.0
     694    //
     695    for (psS32 y = kernelSize; y < numRows - kernelSize; y++) {
     696        for (psS32 x = 0; x < kernelSize; x++) {
     697            convolved->data.F32[y][x] = 0.0;
     698        }
     699        for (psS32 x = numCols - kernelSize; x < numCols; x++) {
     700            convolved->data.F32[y][x] = 0.0;
     701        }
     702    }
     703    for (psS32 y = 0; y < kernelSize; y++) {
     704        for (psS32 x = 0; x < numCols; x++) {
     705            convolved->data.F32[y][x] = 0.0;
     706        }
     707    }
     708    for (psS32 y = numRows - kernelSize; y < numRows; y++) {
     709        for (psS32 x = 0; x < numCols; x++) {
     710            convolved->data.F32[y][x] = 0.0;
     711        }
     712    }
     713
     714    psTrace("ImageSubtract.ConvolveImage", 5, "Exiting ConvolveImage()\n");
     715    return convolved;
     716}
     717
     718
     719
     720
     721
     722/*******************************************************************************
     723XXX: We should assert that the (footprint, kernelSize, imageSize) stuff
     724ensures that all data is accessed in bounds?
    521725 ******************************************************************************/
    522726bool pmSubtractionCalculateEquation(psArray *stamps,          ///< The stamps for which to calculate the equation,
    … …  
    527731                                   )
    528732{
     733    psTrace("ImageSubtract.pmSubtractionCalculateEquation", 1,
     734            "Calling pmSubtractionCalculateEquation()\n");
    529735    PS_ASSERT_PTR_NON_NULL(stamps, false);
    530736    PS_ASSERT_IMAGE_NON_NULL(reference, false);
    … …  
    623829            psTrace("pmSubtractionCalculateEquation", 5, "subCalcEqn(): stamp %d: generated spatial order terms.\n", s);
    624830
    625             //
    626             // Iterate over all pixels surrounding this stamp.
    627             //
    628             for (psS32 y = stamp->y - footprint; y < stamp->y + footprint; y++) {
    629                 for (psS32 x = stamp->x - footprint; x < stamp->x + footprint; x++) {
    630 
    631                     psF32 invNoise2 = 1.0/reference->data.F32[y][x]; // The inverse of the noise, squared.
    632 
    633                     // XXX: Think about a better decomposition of the computation
    634                     // based on whether the kernels are ISIS or POIS.
    635 
    636                     psTrace("pmSubtractionCalculateEquation", 5, "subCalcEqn(): pixel (%d, %d).\n", y, x);
    637                     if (kernels->type == PM_SUBTRACTION_KERNEL_POIS) {
     831            if (kernels->type == PM_SUBTRACTION_KERNEL_POIS) {
     832                //
     833                // Iterate over all pixels surrounding this stamp.
     834                //
     835                for (psS32 y = stamp->y - footprint; y < stamp->y + footprint; y++) {
     836                    for (psS32 x = stamp->x - footprint; x < stamp->x + footprint; x++) {
     837                        psTrace("pmSubtractionCalculateEquation", 5, "subCalcEqn(): pixel (%d, %d).\n", y, x);
     838
     839                        // The inverse of the noise, squared.
     840                        psF32 invNoise2 = 1.0/reference->data.F32[y][x];
     841
    638842                        //
    639843                        // Iterate over the first convolution */
    … …  
    647851                            //
    648852                            // First convolution.  This will set the value for the stampVector.
    649                             //
    650853                            //
    651854                            // XXX: verify the [y-v2][x-u2] subscript.  This generated errors in
    … …  
    705908                        stampMatrix->data.F64[bgIndex][bgIndex] += invNoise2;
    706909                        stampVector->data.F64[bgIndex] += input->data.F32[y][x] * invNoise2;
    707 
    708                     } else if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) {
     910                    }
     911                }
     912            } else if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) {
     913                //
     914                // Iterate over all pixels surrounding this stamp.
     915                //
     916                // XXX: Why isn't there a polyValues term here?
     917                //
     918
     919                for (psS32 y = stamp->y - footprint; y < stamp->y + footprint; y++) {
     920                    for (psS32 x = stamp->x - footprint; x < stamp->x + footprint; x++) {
     921                        psTrace("pmSubtractionCalculateEquation", 5, "subCalcEqn(): pixel (%d, %d).\n", y, x);
     922                        psF32 invNoise2 = 1.0/reference->data.F32[y][x]; // The inverse of the noise, squared.
     923
    709924                        for (psS32 k1 = 0; k1 < numKernels; k1++) {
    710                             psF32 conv1 = GeneralKernelConvolve(reference, kernels, k1, x, y);
     925                            psF32 conv1 = IsisKernelConvolve(reference, kernels, k1, x, y);
    711926
    712927                            for (psS32 k2 = k1; k2 < numKernels; k2++) {
    713928                                //printf("(k1, k2) is (%d, %d)\n", k1, k2);
    714                                 psF32 conv2 = GeneralKernelConvolve(reference, kernels, k2, x, y);
    715                                 stampMatrix->data.F64[k1][k2] += conv1 * conv2 * invNoise2;
     929                                psF32 conv2 = IsisKernelConvolve(reference, kernels, k2, x, y);
     930                                stampMatrix->data.F64[k1][k2]+= conv1 * conv2 * invNoise2;
    716931                            }
    717                             stampVector->data.F64[k1] += input->data.F32[y][x] * conv1 * invNoise2;
     932                            stampVector->data.F64[k1]+= input->data.F32[y][x] * conv1 * invNoise2;
    718933                            stampMatrix->data.F64[k1][bgIndex] += conv1 * invNoise2;
    719934                        }
    720                     } else {
    721                         psLogMsg(__func__, PS_LOG_WARN, "WARNING: unknown kernel->type.\n");
    722                         return(false);
    723                     } // if-else on kernel type.
    724                 } // y-loop for all pixels.
    725             } // x-loop for all pixels.
     935                        stampMatrix->data.F64[bgIndex][bgIndex] += invNoise2;
     936                        stampVector->data.F64[bgIndex] += input->data.F32[y][x] * invNoise2;
     937                    }
     938                }
     939            } else {
     940                psLogMsg(__func__, PS_LOG_WARN, "WARNING: unknown kernel->type.\n");
     941                return(false);
     942            }
    726943            psFree(polyValues);
     944
     945            // XXX: Generate psTrace()
     946            if (0) {
     947                for (psS32 s = 0; s < stamps->n; s++) {
     948                    pmStamp *stamp = (pmStamp *) stamps->data[s];
     949                    if (stamp->status == PM_STAMP_RECALC) {
     950                        psVector *stampVector = stamp->vector;
     951                        printf("COOL: stamp %d vector:\n", s);
     952                        PS_VECTOR_PRINT_F64(stampVector);
     953                    }
     954                }
     955            }
    727956
    728957            //
    … …  
    756985        }
    757986    }
     987    psTrace("ImageSubtract.pmSubtractionCalculateEquation", 1,
     988            "Exiting pmSubtractionCalculateEquation()\n");
    758989    return(true);
    759990}
    … …  
    763994
    764995/*******************************************************************************
    765 XXX: Assert correct vector matrix sizes.
    766996 ******************************************************************************/
    767997psVector *pmSubtractionSolveEquation(psVector *solution, ///< Solution vector, or NULL
    … …  
    769999                                    )
    7701000{
     1001    psTrace("ImageSubtract.pmSubtractionSolveEquation", 1,
     1002            "Calling pmSubtractionSolveEquation()\n");
    7711003    PS_ASSERT_PTR_NON_NULL(stamps, NULL);
    7721004    psS32 size = -1;
    … …  
    8021034    //
    8031035    // XXX: Test these functions with size=-1.  This caused seg faults during test.
     1036    //      This should be done in the psImage.c and psVector.c test files.  It
     1037    //      should never occur here.
    8041038    //
    8051039    psImage *sumMatrix = psImageAlloc(size, size, PS_TYPE_F64);
    … …  
    8291063        }
    8301064    }
    831 
     1065    psVector *permutation = NULL;
    8321066    // XXX: Check output from these routines.
    833     psVector *permutation = NULL;
     1067
     1068    // XXX: psTrace()
     1069    if (1) {
     1070        PS_IMAGE_PRINT_F64(sumMatrix);
     1071    }
     1072
    8341073    psImage *luMatrix = psMatrixLUD(NULL, &permutation, sumMatrix);
     1074    if (luMatrix == NULL) {
     1075        psError(PS_ERR_UNKNOWN, true, "Failed to LU-Decompose the matrix.\n");
     1076        psFree(sumMatrix);
     1077        psFree(sumVector);
     1078        psFree(luMatrix);
     1079        psFree(permutation);
     1080        return(NULL);
     1081    }
     1082    // XXX: psTrace()
     1083    if (1) {
     1084        PS_IMAGE_PRINT_F64(luMatrix);
     1085    }
     1086
    8351087    solution = psMatrixLUSolve(solution, luMatrix, sumVector, permutation);
     1088    // XXX: psTrace()
     1089    if (1) {
     1090        PS_VECTOR_PRINT_F64(solution);
     1091    }
     1092    if (solution == NULL) {
     1093        psError(PS_ERR_UNKNOWN, true, "Failed to solve the matrix.\n");
     1094        psFree(sumMatrix);
     1095        psFree(sumVector);
     1096        psFree(luMatrix);
     1097        psFree(permutation);
     1098        return(NULL);
     1099    }
    8361100
    8371101    psFree(sumMatrix);
    … …  
    8401104    psFree(permutation);
    8411105
     1106    psTrace("ImageSubtract.pmSubtractionSolveEquation", 1,
     1107            "Exiting pmSubtractionSolveEquation()\n");
    8421108    return(solution);
    843 }
    844 
    845 
    846 /*******************************************************************************
    847 ConvolvePixelPois(input, mask, badStampMaskVal, solution, kernels, col, row):
    848  
    849 This routine takes a single pixel in the psImage input and convolves it with
    850 the set of kernel basis functions and their appropriate weights in solution.
    851 It returns the value of the convolved pixel.
    852  
    853 XXX: Static structure for polyValues?
    854  ******************************************************************************/
    855 static psF32 ConvolvePixelPois(const psImage *input,
    856                                const psImage *mask,
    857                                psU32 badStampMaskVal,
    858                                const psVector *solution,
    859                                const psSubtractionKernels *kernels,
    860                                psS32 col,
    861                                psS32 row)
    862 {
    863     psS32 nBF = kernels->u->n;
    864     psF32 numColsHalf = 0.5 * (psF32) input->numCols;
    865     psF32 numRowsHalf = 0.5 * (psF32) input->numRows;
    866     psF32 background = solution->data.F64[solution->n-1];
    867     psS32 spatialOrder = kernels->p_spatialOrder;
    868     psF32 conv = background; // Initial convolved value
    869 
    870     if ((mask == NULL) || !(mask->data.U8[row][col] & badStampMaskVal)) {
    871         psF32 imageX = (((psF32) col) - numColsHalf) / numColsHalf; // Normalised position
    872         psF32 imageY = (((psF32) row) - numRowsHalf) / numRowsHalf; // Normalised position
    873         psImage *polyValues = GenSpatialOrder(spatialOrder, imageX, imageY);
    874 
    875         // Iterate over the kernel basis functions
    876         for (psS32 k = 0; k < nBF; k++) {
    877             psS32 u = (psS32) kernels->u->data.F32[k];
    878             psS32 v = (psS32) kernels->v->data.F32[k];
    879 
    880             // XXX: What's the story with this?
    881             #if 0
    882 
    883             psS32 xOrder = (psS32) kernels->xOrder->data.F32[k];
    884             psS32 yOrder = (psS32) kernels->yOrder->data.F32[k];
    885             psF32 polyVal = polyValues->data.F32[yOrder][xOrder];
    886             #else
    887 
    888             psF32 polyVal = 1.0;
    889             #endif
    890 
    891             // XXX: Why this?
    892             if (k == 0) {
    893                 conv += solution->data.F64[k] * input->data.F32[row - v][col - u] * polyVal;
    894             } else {
    895                 conv += solution->data.F64[k] *
    896                         (input->data.F32[row - v][col - u] * polyVal - input->data.F32[row][col]);
    897             }
    898         }
    899         psFree(polyValues);
    900     }
    901 
    902     return(conv);
    903 }
    904 
    905 
    906 
    907 /*******************************************************************************
    908 ConvolvePixelIsis(input, mask, badStampMaskVal, solution, kernels, col, row):
    909  
    910 This routine takes a single pixel in the psImage input and convolves it with
    911 the set of kernel basis functions and their appropriate weights in solution.
    912 It returns the value of the convolved pixel.
    913  
    914 XXX: Static structure for polyValues?
    915  ******************************************************************************/
    916 static psF32 ConvolvePixelIsis(const psImage *input,
    917                                const psImage *mask,
    918                                psU32 badStampMaskVal,
    919                                const psVector *solution,
    920                                const psSubtractionKernels *kernels,
    921                                psS32 col,
    922                                psS32 row)
    923 {
    924     psS32 nBF = kernels->u->n;
    925     psF32 numColsHalf = 0.5 * (psF32) input->numCols;
    926     psF32 numRowsHalf = 0.5 * (psF32) input->numRows;
    927     psF32 background = solution->data.F64[solution->n-1];
    928     psS32 spatialOrder = kernels->p_spatialOrder;
    929     psS32 kernelSize = kernels->p_size;
    930     psF32 conv = background; // Initial convolved value
    931 
    932     if ((mask == NULL) || !(mask->data.U8[row][col] & badStampMaskVal)) {
    933         psF32 imageX = (((psF32) col) - numColsHalf) / numColsHalf; // Normalised position
    934         psF32 imageY = (((psF32) row) - numRowsHalf) / numRowsHalf; // Normalised position
    935         psImage *polyValues = GenSpatialOrder(spatialOrder, imageX, imageY);
    936 
    937         // Iterate over the kernel basis functions
    938         for (psS32 k = 0; k < nBF; k++) {
    939             psS32 xOrder = (psS32) kernels->xOrder->data.F32[k];
    940             psS32 yOrder = (psS32) kernels->yOrder->data.F32[k];
    941             psF32 polyVal = polyValues->data.F32[yOrder][xOrder];
    942 
    943             psImage *preCalc = (psImage *) kernels->preCalc->data[k];
    944             for (psS32 yy = -kernelSize ; yy < kernelSize ; yy++) {
    945                 for (psS32 xx = -kernelSize ; xx < kernelSize ; xx++) {
    946                     // XXX: Should I do something special for the k==0 kernel?
    947                     conv += solution->data.F64[k] *
    948                             input->data.F32[yy+row][xx+col] *
    949                             preCalc->data.F32[yy+kernelSize][xx+kernelSize] *
    950                             polyVal;
    951                 }
    952             }
    953         }
    954         psFree(polyValues);
    955     }
    956 
    957     return(conv);
    958 }
    959 
    960 /*******************************************************************************
    961 ConvolveImage(input, mask, badStampMaskVal, solution, kernels): convolves an
    962 arbitrary image with either an ISIS or POIS set of kernel basis functions.
    963  ******************************************************************************/
    964 static psImage *ConvolveImage(const psImage *input,
    965                               const psImage *mask,
    966                               psU32 badStampMaskVal,
    967                               const psVector *solution,
    968                               const psSubtractionKernels *kernels)
    969 {
    970     PS_ASSERT_IMAGE_NON_NULL(input, NULL);
    971     PS_ASSERT_IMAGE_NON_EMPTY(input, NULL);
    972     PS_ASSERT_IMAGE_TYPE(input, PS_TYPE_F32, NULL);
    973     if (mask != NULL) {
    974         PS_ASSERT_IMAGES_SIZE_EQUAL(input, mask, NULL);
    975         PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_U8, NULL);
    976     }
    977     PS_ASSERT_VECTOR_NON_NULL(solution, NULL);
    978     PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, NULL);
    979     PS_ASSERT_PTR_NON_NULL(kernels, NULL);
    980     PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, NULL);
    981     PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, NULL);
    982     PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, NULL);
    983     if (kernels->preCalc != NULL) {
    984         PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, NULL);
    985     } else {
    986         if (kernels->type != PM_SUBTRACTION_KERNEL_POIS) {
    987             psError(PS_ERR_BAD_PARAMETER_NULL, true,
    988                     "Unallowable operation: kernels->preCalc == NULL and kernels->type != PM_SUBTRACTION_KERNEL_POIS.\n");
    989             return(NULL);
    990         }
    991     }
    992     psS32 nBF = kernels->u->n;
    993     PS_ASSERT_VECTOR_SIZE(solution, nBF+1, NULL);
    994 
    995     psS32 numCols = input->numCols;
    996     psS32 numRows = input->numRows;
    997     psS32 kernelSize = kernels->p_size;
    998 
    999     psImage *convolved = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    1000 
    1001     for (psS32 y = kernelSize; y < numRows - kernelSize; y++) {
    1002         for (psS32 x = kernelSize; x < numCols - kernelSize; x++) {
    1003             if (kernels->type == PM_SUBTRACTION_KERNEL_POIS) {
    1004                 convolved->data.F32[y][x] = ConvolvePixelPois(input, mask, badStampMaskVal,
    1005                                             solution, kernels, x, y);
    1006             } else if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) {
    1007                 convolved->data.F32[y][x] = ConvolvePixelIsis(input, mask, badStampMaskVal,
    1008                                             solution, kernels, x, y);
    1009             } else {
    1010                 psLogMsg(__func__, PS_LOG_WARN, "WARNING: unknown kernel type.  Returning NULL\n");
    1011                 return(NULL);
    1012             }
    1013         }
    1014     }
    1015 
    1016     //
    1017     // Pad the rest of the convolved image with 0.0
    1018     //
    1019     for (psS32 y = kernelSize; y < numRows - kernelSize; y++) {
    1020         for (psS32 x = 0; x < kernelSize; x++) {
    1021             convolved->data.F32[y][x] = 0.0;
    1022         }
    1023         for (psS32 x = numCols - kernelSize; x < numCols; x++) {
    1024             convolved->data.F32[y][x] = 0.0;
    1025         }
    1026     }
    1027     for (psS32 y = 0; y < kernelSize; y++) {
    1028         for (psS32 x = 0; x < numCols; x++) {
    1029             convolved->data.F32[y][x] = 0.0;
    1030         }
    1031     }
    1032     for (psS32 y = numRows - kernelSize; y < numRows; y++) {
    1033         for (psS32 x = 0; x < numCols; x++) {
    1034             convolved->data.F32[y][x] = 0.0;
    1035         }
    1036     }
    1037 
    1038     return convolved;
    10391109}
    10401110
    … …  
    10521122                                     const psVector *solution)
    10531123{
     1124    psTrace("ImageSubtract.CalculateDeviations", 5,
     1125            "Calling CalculateDeviations()\n");
    10541126    PS_ASSERT_PTR_NON_NULL(stamps, NULL);
    10551127    if (deviations != NULL) {
    … …  
    11151187            psImage *subStampTrim = psImageSubset((psImage *) subStamp, myReg);
    11161188            psImage *maskStampTrim = psImageSubset((psImage *) maskStamp, myReg);
    1117             // Copy image to workaround bug 305
    1118             // XXX: Modify.  This bug is now fixed.
    1119             //            psImage *tempImage = psImageCopy(NULL, subStampTrim, PS_TYPE_F32);
    1120             //            psImage *tempMask = psImageCopy(NULL, maskStampTrim, PS_TYPE_U8);
    1121             //            (void)psImageStats(stats, tempImage, tempMask, badStampMaskVal);
    1122             (void)psImageStats(stats, subStampTrim, maskStampTrim, badStampMaskVal);
     1189            psImageStats(stats, subStampTrim, maskStampTrim, badStampMaskVal);
    11231190
    11241191            deviations->data.F32[s] = stats->sampleMean * (psF32)footprint * (psF32)footprint * 4.0;
    … …  
    11361203    psFree(subStamp);
    11371204
     1205    psTrace("ImageSubtract.CalculateDeviations", 5,
     1206            "Exiting CalculateDeviations()\n");
    11381207    return deviations;
    11391208}
    … …  
    11521221                              )
    11531222{
     1223    psTrace("ImageSubtract.pmSubtractionRejectStamps", 5,
     1224            "Calling pmSubtractionRejectStamps()\n");
    11541225    PS_ASSERT_PTR_NON_NULL(stamps, false);
    11551226    PS_ASSERT_IMAGE_NON_NULL(refImage, false);
    … …  
    12251296
    12261297    psFree(deviations);
     1298    psTrace("ImageSubtract.pmSubtractionRejectStamps", 5,
     1299            "Exiting pmSubtractionRejectStamps()\n");
    12271300    return(true);
    12281301}
    … …  
    12371310                                 )
    12381311{
     1312    psTrace("ImageSubtract.pmSubtractionKernelImage", 5,
     1313            "Calling pmSubtractionKernelImage()\n");
    12391314    PS_ASSERT_VECTOR_NON_NULL(solution, NULL);
    12401315    PS_ASSERT_PTR_NON_NULL(kernels, NULL);
    … …  
    12801355            psS32 xOrder = (psS32) kernels->xOrder->data.F32[k];
    12811356            psS32 yOrder = (psS32) kernels->yOrder->data.F32[k];
    1282             psF32 polyVal = polyValues->data.F32[yOrder][xOrder];
     1357            psF64 polyVal = polyValues->data.F64[yOrder][xOrder];
    12831358
    12841359            // XXX: Verify that this is correct.
    … …  
    13041379
    13051380            out->data.F32[kernelSize - v][kernelSize - u]+=
    1306                 solution->data.F64[k] * polyValues->data.F32[yOrder][xOrder];
     1381                solution->data.F64[k] * polyValues->data.F64[yOrder][xOrder];
    13071382        }
    13081383    }
    13091384    psFree(polyValues);
    13101385
     1386    psTrace("ImageSubtract.pmSubtractionKernelImage", 5,
     1387            "Exiting pmSubtractionKernelImage()\n");
    13111388    return(out);
    13121389}

  • trunk/psModules/src/pmImageSubtract.h

    r4322 r4417  
    11/** @file  ImageSubtract.h
    22 *
    3  *  This file will ...
     3 *  This file will contain code which creates a set of kernel basis
     4 *  functions, solves for their solution, and applies them to an image.
    45 *
     6 *  @author Paul Price, IfA (original prototype)
    57 *  @author GLG, MHPCC
    68 *
    7  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2005-06-20 23:21:05 $
     9 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2005-06-28 23:23:55 $
    911 *
    1012 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2628typedef enum {
    2729    PM_SUBTRACTION_KERNEL_POIS,    ///< POIS kernel --- delta functions
    28     PM_SUBTRACTION_KERNEL_ISIS    ///< ISIS kernel --- gaussians modified by polynomials
     30    PM_SUBTRACTION_KERNEL_ISIS     ///< ISIS kernel --- gaussians modified by polynomials
    2931} pmSubtractionKernelsType;
    3032
    … …  
    3234{
    3335    pmSubtractionKernelsType type;   ///< Type ofKernels --- allowing the use of multiple kernels
    34     psVector *u, *v;     ///< Offset (for POIS) or polynomial order (for ISIS)
    35     psVector *sigma;     ///< Width of Gaussian (for ISIS)
    36     psVector *xOrder, *yOrder;    ///< Spatial Polynomial order (for all)
    37     int subIndex;     ///< Index of kernel to be subtracted to maintain flux conservation
    38     psArray *preCalc;     ///< Array of images containing pre-calculated kernel (to
     36    psVector *u, *v;                 ///< Offset (for POIS) or polynomial order (for ISIS)
     37    psVector *sigma;                 ///< Width of Gaussian (for ISIS)
     38    psVector *xOrder, *yOrder;       ///< Spatial Polynomial order (for all)
     39    int subIndex;                    ///< Index of kernel to be subtracted to maintain flux conservation
     40    psArray *preCalc;                ///< Array of images containing pre-calculated kernel (to
    3941    ///< accelerate ISIS; don't use for POIS)
    40     psS32 p_size;                                       ///< The halfsize of the kernel
    41     psS32 p_spatialOrder;                               ///< The spatial order of the kernels
     42    psS32 p_size;                    ///< The halfsize of the kernel
     43    psS32 p_spatialOrder;            ///< The spatial order of the kernels
    4244}
    4345psSubtractionKernels;
    … …  
    5254
    5355typedef enum {
    54     PM_STAMP_INIT,     ///< Initial state
    55     PM_STAMP_USED,     ///< Use this stamp
    56     PM_STAMP_REJECTED,     ///< This stamp has been rejected
    57     PM_STAMP_RECALC,     ///< Having been reset, this stamp is to be recalculated
    58     PM_STAMP_NONE     ///< No stamp in this region
     56    PM_STAMP_INIT,                   ///< Initial state
     57    PM_STAMP_USED,                   ///< Use this stamp
     58    PM_STAMP_REJECTED,               ///< This stamp has been rejected
     59    PM_STAMP_RECALC,                 ///< Having been reset, this stamp is to be recalculated
     60    PM_STAMP_NONE                    ///< No stamp in this region
    5961} pmStampStatus;
    6062
Note: See TracChangeset for help on using the changeset viewer.