IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26486


Ignore:
Timestamp:
Dec 27, 2009, 8:21:02 AM (16 years ago)
Author:
eugene
Message:

set up deconvolved hermite basis functions (not quite ready)

Location:
branches/eam_branches/20091201/psModules/src
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionDeconvolve.c

    r26479 r26486  
    66#include <string.h>
    77#include <strings.h>
     8#include <fftw3.h>
    89#include <pslib.h>
    910
     11#include "pmFPA.h"
     12#include "pmSubtractionKernels.h"
    1013#include "pmSubtractionDeconvolve.h"
     14#include "pmSubtractionStamps.h"
     15#include "pmSubtractionVisual.h"
     16
     17// Lock FFTW access
     18#define FFTW_LOCK \
     19if (threaded) { \
     20    psFFTLock(); \
     21}
     22// Unlock FFTW access
     23#define FFTW_UNLOCK \
     24if (threaded) { \
     25    psFFTUnlock(); \
     26}
     27
     28#define FFTW_PLAN_RIGOR FFTW_ESTIMATE   // How rigorous the FFTW planning is
    1129
    1230psKernel *pmSubtractionDeconvolveGauss (int size, float sigma) {
    … …  
    1533
    1634    // build the gaussian from 2 1-D Gaussians
    17     psVector *vector = subtractionKernelISIS(sigma, 0, size);
     35    psVector *vector = pmSubtractionKernelISIS(sigma, 0, size);
    1836
    1937    // generate 2D kernel, calculate moments
    … …  
    6381    FFTW_UNLOCK;
    6482
     83    size_t numBytes = numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes per image row
     84
    6585    // copy data from kernelTarg image to dataTarg array
    66     size_t numBytes = numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes per image row
    6786    for (int y = 0; y < numRows; y++) {
    6887        memcpy(&dataTarg[y*numCols], kernelTarg->image->data.F32[y], numBytes);
    6988    }
    70 
    7189    // copy data from kernelConv image to dataConv array
    72     size_t numBytes = numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes per image row
    7390    for (int y = 0; y < numRows; y++) {
    7491        memcpy(&dataConv[y*numCols], kernelConv->image->data.F32[y], numBytes);
    … …  
    7895    // Note that the FFT images have different size from the input
    7996    FFTW_LOCK;
    80     fftwf_complex *fftTarg = fftwf_malloc((numCols/2 + 1) * paddedRows * sizeof(fftwf_complex)); // FFT
    81     fftwf_complex *fftConv = fftwf_malloc((numCols/2 + 1) * paddedRows * sizeof(fftwf_complex)); // FFT
    82     FFTW_UNLOCK;
    83 
    84     FFTW_LOCK;
    85     fftwf_plan forwardTarg = fftwf_plan_dft_r2c_2d(numRows, numCols, dataTarg, fftTarg, FFT_PLAN_RIGOR);
    86     fftwf_plan forwardConv = fftwf_plan_dft_r2c_2d(numRows, numCols, dataConv, fftConv, FFT_PLAN_RIGOR);
     97    fftwf_complex *fftTarg = fftwf_malloc((numCols/2 + 1) * numRows * sizeof(fftwf_complex)); // FFT
     98    fftwf_complex *fftConv = fftwf_malloc((numCols/2 + 1) * numRows * sizeof(fftwf_complex)); // FFT
     99    FFTW_UNLOCK;
     100
     101    FFTW_LOCK;
     102    fftwf_plan forwardTarg = fftwf_plan_dft_r2c_2d(numRows, numCols, dataTarg, fftTarg, FFTW_PLAN_RIGOR);
     103    fftwf_plan forwardConv = fftwf_plan_dft_r2c_2d(numRows, numCols, dataConv, fftConv, FFTW_PLAN_RIGOR);
    87104    FFTW_UNLOCK;
    88105
    … …  
    103120    // but anywhere Cr^2 - Ci^2 < 1e-7 of the max, mask it
    104121
     122    // the X dimension is halved by FFTW
     123    int numColsOut = numCols / 2 + 1;
     124
    105125    // generate Det = Cr^2 - Ci^2
    106126    float maxValue = 0.0;
    107     psImage *det = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     127    psImage *det = psImageAlloc(numColsOut, numRows, PS_TYPE_F32);
    108128    for (int iy = 0; iy < numRows; iy++) {
    109         for (int ix = 0; ix < numCols; ix++) {
    110             float convReal = fftConv[ix + iy*numCols][0];
    111             float convImag = fftConv[ix + iy*numCols][1];
     129        for (int ix = 0; ix < numColsOut; ix++) {
     130            float convReal = fftConv[ix + iy*numColsOut][0];
     131            float convImag = fftConv[ix + iy*numColsOut][1];
    112132            det->data.F32[iy][ix] = convReal*convReal - convImag*convImag;
    113             maxValue = PS_MAX(fabs(det->data.F32[iy][ix]));
    114         }
    115     }
     133            maxValue = PS_MAX(maxValue, fabs(det->data.F32[iy][ix]));
     134        }
     135    }
     136# define TOL 1e-7
    116137    float limit = TOL*maxValue;
    117138
    118139    // generate Deco = targ * conv^* / (Cr^2 - Ci^2)
    119140    for (int iy = 0; iy < numRows; iy++) {
    120         for (int ix = 0; ix < numCols; ix++) {
    121             float targReal = fftTarg[ix + iy*numCols][0];
    122             float targImag = fftTarg[ix + iy*numCols][1];
    123             float convReal = fftConv[ix + iy*numCols][0];
    124             float convImag = fftConv[ix + iy*numCols][1];
     141        for (int ix = 0; ix < numColsOut; ix++) {
     142            float targReal = fftTarg[ix + iy*numColsOut][0];
     143            float targImag = fftTarg[ix + iy*numColsOut][1];
     144            float convReal = fftConv[ix + iy*numColsOut][0];
     145            float convImag = fftConv[ix + iy*numColsOut][1];
    125146            if (fabs(det->data.F32[iy][ix]) < limit) {
    126                 fftTarg[ix + iy*numCols][0] = 0.0;
    127                 fftTarg[ix + iy*numCols][1] = 0.0;
     147                fftTarg[ix + iy*numColsOut][0] = 0.0;
     148                fftTarg[ix + iy*numColsOut][1] = 0.0;
    128149            } else {
    129                 fftTarg[ix + iy*numCols][0] = targReal*convReal + targImag*convImag;
    130                 fftTarg[ix + iy*numCols][1] = targImag*convReal - targReal*convImag;
     150                fftTarg[ix + iy*numColsOut][0] = targReal*convReal + targImag*convImag;
     151                fftTarg[ix + iy*numColsOut][1] = targImag*convReal - targReal*convImag;
    131152            }
    132153        }
    … …  
    147168
    148169    psKernel *output = psKernelAlloc (kernelTarg->xMin, kernelTarg->xMax, kernelTarg->yMin, kernelTarg->yMax);
    149     for (int y = 0; y < numRows; y++, outData++, dataPtr += paddedCols) {
     170    for (int y = 0; y < numRows; y++) {
    150171        memcpy(output->image->data.F32[y], &dataTarg[y*numCols], numBytes);
    151172    }
    … …  
    157178
    158179    return output;
     180
     181 escape:
     182    psError(PS_ERR_BAD_PARAMETER_VALUE, true, "mismatch between kernel and image");
     183    return NULL;
     184
    159185}
    160186
    161187bool pmSubtractionDeconvolutionTest () {
    162188
     189    float sigma = 2.0;
    163190    int size = 15;
    164191
    165192    // generate a Hermite polynomial
    166     psVector *xKernel = subtractionKernelHERM(sigma, 2, size); // x Kernel
    167     psVector *yKernel = subtractionKernelHERM(sigma, 2, size); // y Kernel
     193    psVector *xKernel = pmSubtractionKernelHERM(sigma, 2, size); // x Kernel
     194    psVector *yKernel = pmSubtractionKernelHERM(sigma, 2, size); // y Kernel
    168195    psKernel *kernelTarget = psKernelAlloc(-size, size, -size, size);   // Kernel
    169196
    … …  
    183210
    184211    // re-convolve the kernel
    185     psImage *kernelConv = psImageConvolveFFT(NULL, kernel, NULL, 0, kernelGauss);
    186     pmSubtractionVisualShowSubtraction (kernelTarget->image, kernel->image, kernelConv->image);
     212    psImage *kernelConv = psImageConvolveFFT(NULL, kernel->image, NULL, 0, kernelGauss);
     213    pmSubtractionVisualShowSubtraction (kernelTarget->image, kernel->image, kernelConv);
    187214
    188215    return true;

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c

    r26479 r26486  
    4545
    4646// Generate 1D convolution kernel for ISIS
    47 static psVector *subtractionKernelISIS(float sigma, // Gaussian width
     47psVector *pmSubtractionKernelISIS(float sigma, // Gaussian width
    4848                                       int order, // Polynomial order
    4949                                       int size // Kernel half-size
    … …  
    6363
    6464// Generate 1D convolution kernel for HERM (normalized for 2D)
    65 static psVector *subtractionKernelHERM(float sigma, // Gaussian width
     65psVector *pmSubtractionKernelHERM(float sigma, // Gaussian width
    6666                                       int order, // Polynomial order
    6767                                       int size // Kernel half-size
    … …  
    179179            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
    180180                psArray *preCalc = psArrayAlloc(3); // Array to hold precalculated values
    181                 psVector *xKernel = preCalc->data[0] = subtractionKernelISIS(sigma, uOrder, size); // x Kernel
    182                 psVector *yKernel = preCalc->data[1] = subtractionKernelISIS(sigma, vOrder, size); // y Kernel
     181                psVector *xKernel = preCalc->data[0] = pmSubtractionKernelISIS(sigma, uOrder, size); // x Kernel
     182                psVector *yKernel = preCalc->data[1] = pmSubtractionKernelISIS(sigma, vOrder, size); // y Kernel
    183183                psKernel *kernel = preCalc->data[2] = psKernelAlloc(-size, size, -size, size);      // Kernel
    184184
    … …  
    283283            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
    284284                psArray *preCalc = psArrayAlloc(3); // Array to hold precalculated values
    285                 psVector *xKernel = preCalc->data[0] = subtractionKernelHERM(sigma, uOrder, size); // x Kernel
    286                 psVector *yKernel = preCalc->data[1] = subtractionKernelHERM(sigma, vOrder, size); // y Kernel
     285                psVector *xKernel = preCalc->data[0] = pmSubtractionKernelHERM(sigma, uOrder, size); // x Kernel
     286                psVector *yKernel = preCalc->data[1] = pmSubtractionKernelHERM(sigma, vOrder, size); // y Kernel
    287287                psKernel *kernel = preCalc->data[2] = psKernelAlloc(-size, size, -size, size);      // Kernel
    288288
    … …  
    398398            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
    399399                psArray *preCalc  = psArrayAlloc(3); // Array to hold precalculated values
    400                 psVector *xKernel = preCalc->data[0] = subtractionKernelHERM(sigma, uOrder, size); // x Kernel
    401                 psVector *yKernel = preCalc->data[1] = subtractionKernelHERM(sigma, vOrder, size); // y Kernel
     400                psVector *xKernel = preCalc->data[0] = pmSubtractionKernelHERM(sigma, uOrder, size); // x Kernel
     401                psVector *yKernel = preCalc->data[1] = pmSubtractionKernelHERM(sigma, vOrder, size); // y Kernel
    402402                psKernel *kernelTarget = psKernelAlloc(-size, size, -size, size);       // Kernel
    403403

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.h

    r26479 r26486  
    111111}
    112112
     113psVector *pmSubtractionKernelISIS(float sigma, // Gaussian width
     114                                       int order, // Polynomial order
     115                                       int size // Kernel half-size
     116    );
     117
     118psVector *pmSubtractionKernelHERM(float sigma, // Gaussian width
     119                                       int order, // Polynomial order
     120                                       int size // Kernel half-size
     121    );
     122
    113123/// Generate a delta-function grid for subtraction kernels (like the POIS kernel)
    114124bool p_pmSubtractionKernelsAddGrid(pmSubtractionKernels *kernels, ///< The subtraction kernels to append to

  • branches/eam_branches/20091201/psModules/src/psmodules.h

    r25383 r26486  
    9696#include <pmSubtractionStamps.h>
    9797#include <pmSubtractionKernels.h>
     98#include <pmSubtractionDeconvolve.h>
    9899#include <pmSubtractionAnalysis.h>
    99100#include <pmSubtractionMatch.h>
Note: See TracChangeset for help on using the changeset viewer.