IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 29, 2009, 8:04:00 AM (17 years ago)
Author:
eugene
Message:

modify preCalc to carry a specific structure to avoid confusion of array elements; update to use deconvolved hermitians

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r26486 r26491  
    4949// deconvolve kernelTarget by kernelConv to get the kernel which, when convolved
    5050// by kernelConv results in kernelTarget...
     51// XXX using complex to complex, explicitly setting the imaginary part to zero
    5152psKernel *pmSubtractionDeconvolveKernel (psKernel *kernelTarg, psKernel *kernelConv) {
    5253
    … …  
    7778    // Create data array containing the image and kernel
    7879    FFTW_LOCK;
    79     psF32 *dataTarg = fftwf_malloc(numPixels * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); // Data for FFTW
    80     psF32 *dataConv = fftwf_malloc(numPixels * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); // Data for FFTW
    81     FFTW_UNLOCK;
    82 
    83     size_t numBytes = numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes per image row
     80    // psF32 *dataTarg = fftwf_malloc(numPixels * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); // Data for FFTW
     81    // psF32 *dataConv = fftwf_malloc(numPixels * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); // Data for FFTW
     82    fftwf_complex *dataTarg = fftwf_malloc(numPixels * sizeof(fftwf_complex)); // Data for FFTW
     83    fftwf_complex *dataConv = fftwf_malloc(numPixels * sizeof(fftwf_complex)); // Data for FFTW
     84    FFTW_UNLOCK;
     85
     86    // size_t numBytes = numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes per image row
    8487
    8588    // copy data from kernelTarg image to dataTarg array
    8689    for (int y = 0; y < numRows; y++) {
    87         memcpy(&dataTarg[y*numCols], kernelTarg->image->data.F32[y], numBytes);
    88     }
     90        for (int x = 0; x < numCols; x++) {
     91            dataTarg[x + y*numCols][0] = kernelTarg->image->data.F32[y][x];
     92            dataTarg[x + y*numCols][1] = 0.0;
     93        }
     94    }
     95   
     96    // kernel must be copied to corners of image (0,0 pixel is center of kernel)
    8997    // copy data from kernelConv image to dataConv array
    90     for (int y = 0; y < numRows; y++) {
    91         memcpy(&dataConv[y*numCols], kernelConv->image->data.F32[y], numBytes);
     98    int oy = 0;
     99    for (int iy = 0; iy <= yMax; iy++, oy++) {
     100        int ox = 0;
     101        for (int ix = 0; ix <= xMax; ix++, ox++) {
     102            dataConv[ox + oy*numCols][0] = kernelConv->kernel[iy][ix];
     103            dataConv[ox + oy*numCols][1] = 0.0;
     104        }
     105        for (int ix = xMin; ix <= -1; ix++, ox++) {
     106            dataConv[ox + oy*numCols][0] = kernelConv->kernel[iy][ix];
     107            dataConv[ox + oy*numCols][1] = 0.0;
     108        }
     109    }
     110    for (int iy = yMin; iy <= -1; iy++, oy++) {
     111        int ox = 0;
     112        for (int ix = 0; ix <= xMax; ix++, ox++) {
     113            dataConv[ox + oy*numCols][0] = kernelConv->kernel[iy][ix];
     114            dataConv[ox + oy*numCols][1] = 0.0;
     115        }
     116        for (int ix = xMin; ix <= -1; ix++, ox++) {
     117            dataConv[ox + oy*numCols][0] = kernelConv->kernel[iy][ix];
     118            dataConv[ox + oy*numCols][1] = 0.0;
     119        }
    92120    }
    93121
    … …  
    95123    // Note that the FFT images have different size from the input
    96124    FFTW_LOCK;
    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);
     125    fftwf_complex *fftTarg = fftwf_malloc(numCols * numRows * sizeof(fftwf_complex)); // FFT
     126    fftwf_complex *fftConv = fftwf_malloc(numCols * numRows * sizeof(fftwf_complex)); // FFT
     127    FFTW_UNLOCK;
     128
     129    FFTW_LOCK;
     130    fftwf_plan forwardTarg = fftwf_plan_dft_2d(numRows, numCols, dataTarg, fftTarg, FFTW_FORWARD, FFTW_PLAN_RIGOR);
     131    fftwf_plan forwardConv = fftwf_plan_dft_2d(numRows, numCols, dataConv, fftConv, FFTW_FORWARD, FFTW_PLAN_RIGOR);
    104132    FFTW_UNLOCK;
    105133
    … …  
    121149
    122150    // the X dimension is halved by FFTW
    123     int numColsOut = numCols / 2 + 1;
     151    // int numColsOut = numCols / 2 + 1;
    124152
    125153    // generate Det = Cr^2 - Ci^2
    126154    float maxValue = 0.0;
    127     psImage *det = psImageAlloc(numColsOut, numRows, PS_TYPE_F32);
    128     for (int iy = 0; iy < numRows; iy++) {
    129         for (int ix = 0; ix < numColsOut; ix++) {
    130             float convReal = fftConv[ix + iy*numColsOut][0];
    131             float convImag = fftConv[ix + iy*numColsOut][1];
     155    psImage *det = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     156    psImage *tR  = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     157    psImage *tI  = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     158    psImage *cR  = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     159    psImage *cI  = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     160    for (int iy = 0; iy < numRows; iy++) {
     161        for (int ix = 0; ix < numCols; ix++) {
     162            float convReal = fftConv[ix + iy*numCols][0];
     163            float convImag = fftConv[ix + iy*numCols][1];
    132164            det->data.F32[iy][ix] = convReal*convReal - convImag*convImag;
    133165            maxValue = PS_MAX(maxValue, fabs(det->data.F32[iy][ix]));
    134         }
    135     }
     166
     167            tR->data.F32[iy][ix] = fftTarg[ix + iy*numCols][0];
     168            tI->data.F32[iy][ix] = fftTarg[ix + iy*numCols][1];
     169            cR->data.F32[iy][ix] = fftConv[ix + iy*numCols][0];
     170            cI->data.F32[iy][ix] = fftConv[ix + iy*numCols][1];
     171        }
     172    }
     173
     174    // pmSubtractionVisualShowSubtraction (det, tR, tI);
     175    // pmSubtractionVisualShowSubtraction (det, cR, cI);
     176
     177# if 1
    136178# define TOL 1e-7
    137179    float limit = TOL*maxValue;
    138 
    139180    // generate Deco = targ * conv^* / (Cr^2 - Ci^2)
    140181    for (int iy = 0; iy < numRows; iy++) {
    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];
     182        for (int ix = 0; ix < numCols; ix++) {
     183            float targReal = fftTarg[ix + iy*numCols][0];
     184            float targImag = fftTarg[ix + iy*numCols][1];
     185            float convReal = fftConv[ix + iy*numCols][0];
     186            float convImag = fftConv[ix + iy*numCols][1];
    146187            if (fabs(det->data.F32[iy][ix]) < limit) {
    147                 fftTarg[ix + iy*numColsOut][0] = 0.0;
    148                 fftTarg[ix + iy*numColsOut][1] = 0.0;
     188                fftTarg[ix + iy*numCols][0] = 0.0;
     189                fftTarg[ix + iy*numCols][1] = 0.0;
    149190            } else {
    150                 fftTarg[ix + iy*numColsOut][0] = targReal*convReal + targImag*convImag;
    151                 fftTarg[ix + iy*numColsOut][1] = targImag*convReal - targReal*convImag;
     191                fftTarg[ix + iy*numCols][0] = (targReal*convReal + targImag*convImag) / det->data.F32[iy][ix];
     192                fftTarg[ix + iy*numCols][1] = (targImag*convReal - targReal*convImag) / det->data.F32[iy][ix];
     193                // fftTarg[ix + iy*numCols][0] = (targReal*convReal + targImag*convImag);
     194                // fftTarg[ix + iy*numCols][1] = (targImag*convReal - targReal*convImag);
    152195            }
    153196        }
    154197    }
     198# else
     199    for (int iy = 0; iy < numRows; iy++) {
     200        for (int ix = 0; ix < numCols; ix++) {
     201            float targReal = fftTarg[ix + iy*numCols][0];
     202            float targImag = fftTarg[ix + iy*numCols][1];
     203            float convReal = fftConv[ix + iy*numCols][0];
     204            float convImag = fftConv[ix + iy*numCols][1];
     205            fftTarg[ix + iy*numCols][0] = targReal*convReal - targImag*convImag;
     206            fftTarg[ix + iy*numCols][1] = targImag*convReal + targReal*convImag;
     207        }
     208    }
     209# endif
     210
     211    for (int iy = 0; iy < numRows; iy++) {
     212        for (int ix = 0; ix < numCols; ix++) {
     213            tR->data.F32[iy][ix] = fftTarg[ix + iy*numCols][0];
     214            tI->data.F32[iy][ix] = fftTarg[ix + iy*numCols][1];
     215        }
     216    }
     217    // pmSubtractionVisualShowSubtraction (det, tR, tI);
    155218
    156219    // Do the backward FFT
    157220    FFTW_LOCK;
    158     fftwf_plan backward = fftwf_plan_dft_c2r_2d(numRows, numCols, fftTarg, dataTarg, FFTW_PLAN_RIGOR);
     221    fftwf_plan backward = fftwf_plan_dft_2d(numRows, numCols, fftTarg, dataTarg, FFTW_BACKWARD, FFTW_PLAN_RIGOR);
    159222    FFTW_UNLOCK;
    160223
    … …  
    169232    psKernel *output = psKernelAlloc (kernelTarg->xMin, kernelTarg->xMax, kernelTarg->yMin, kernelTarg->yMax);
    170233    for (int y = 0; y < numRows; y++) {
    171         memcpy(output->image->data.F32[y], &dataTarg[y*numCols], numBytes);
     234        for (int x = 0; x < numCols; x++) {
     235            output->image->data.F32[y][x] = dataTarg[x + y*numCols][0];
     236        }
    172237    }
    173238
    … …  
    185250}
    186251
    187 bool pmSubtractionDeconvolutionTest () {
    188 
    189     float sigma = 2.0;
    190     int size = 15;
     252bool pmSubtractionDeconvolutionTest (int order) {
     253
     254    float sigma = 1.0;
     255    int size = 31;
    191256
    192257    // generate a Hermite polynomial
    193     psVector *xKernel = pmSubtractionKernelHERM(sigma, 2, size); // x Kernel
    194     psVector *yKernel = pmSubtractionKernelHERM(sigma, 2, size); // y Kernel
     258    psVector *xKernel = pmSubtractionKernelHERM(sigma, order, size); // x Kernel
     259    psVector *yKernel = pmSubtractionKernelHERM(sigma, order, size); // y Kernel
    195260    psKernel *kernelTarget = psKernelAlloc(-size, size, -size, size);   // Kernel
    196261
Note: See TracChangeset for help on using the changeset viewer.