IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26491


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

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

Location:
branches/eam_branches/20091201/psModules/src/imcombine
Files:
7 edited

Legend:

Unmodified
Added
Removed

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

    r26474 r26491  
    7272{
    7373    int size = kernels->size;           // Kernel half-size
    74     psArray *preCalc = kernels->preCalc->data[index]; // Precalculated values
     74    pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[index]; // Precalculated values
    7575#if 0
    76     psVector *xKernel = preCalc->data[0]; // Kernel in x
    77     psVector *yKernel = preCalc->data[1]; // Kernel in y
    7876    // Iterating over the kernel
    7977    for (int y = 0, v = -size; v <= size; y++, v++) {
    80         float yValue = value * yKernel->data.F32[y];
     78        float yValue = value * preCalc->yKernel->data.F32[y];
    8179        for (int x = 0, u = -size; u <= size; x++, u++) {
    82             kernel->kernel[v][u] +=  yValue * xKernel->data.F32[x];
     80            kernel->kernel[v][u] +=  yValue * preCalc->xKernel->data.F32[x];
    8381        }
    8482    }
    … …  
    8886    }
    8987#else
    90     psKernel *k = preCalc->data[2]; // Kernel image
    9188    for (int v = -size; v <= size; v++) {
    9289        for (int u = -size; u <= size; u++) {
    93             kernel->kernel[v][u] +=  value * k->kernel[v][u];
     90            kernel->kernel[v][u] +=  value * preCalc->kernel->kernel[v][u];
    9491        }
    9592    }
    … …  
    168165              break;
    169166          }
     167          case PM_SUBTRACTION_KERNEL_DECONV_HERM: {
     168              solvedKernelISIS(kernel, kernels, value, i);
     169              break;
     170          }
    170171          case PM_SUBTRACTION_KERNEL_RINGS: {
    171               psArray *preCalc = kernels->preCalc->data[i]; // Precalculated data
    172               psVector *uCoords = preCalc->data[0]; // u coordinates
    173               psVector *vCoords = preCalc->data[1]; // v coordinates
    174               psVector *poly = preCalc->data[2]; // Polynomial values
    175               int num = uCoords->n;     // Number of pixels
     172              pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[i]; // Precalculated kernels
     173              int num = preCalc->uCoords->n;     // Number of pixels
    176174
    177175              for (int j = 0; j < num; j++) {
    178                   int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; // Kernel coordinates
    179                   kernel->kernel[v][u] += poly->data.F32[j] * value;
     176                  int u = preCalc->uCoords->data.S32[j];
     177                  int v = preCalc->vCoords->data.S32[j]; // Kernel coordinates
     178                  kernel->kernel[v][u] += preCalc->poly->data.F32[j] * value;
    180179              }
    181180              // Photometric scaling is built into the kernel --- no subtraction!
    … …  
    466465    )
    467466{
    468     psArray *preCalc = kernels->preCalc->data[index]; // Precalculated data
     467    pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[index]; // Precalculated data
    469468#if 0
    470469    // Convolving using separable convolution
    471     psVector *xKernel = preCalc->data[0]; // Kernel in x
    472     psVector *yKernel = preCalc->data[1]; // Kernel in y
    473470    int size = kernels->size;     // Size of kernel
    474471
    … …  
    482479            float value = 0.0;    // Value of convolved pixel
    483480            int uMin = x - size, uMax = x + size; // Range for u
    484             psF32 *xKernelData = &xKernel->data.F32[xKernel->n - 1]; // Kernel values
     481            psF32 *xKernelData = &preCalc->xKernel->data.F32[xKernel->n - 1]; // Kernel values
    485482            psF32 *imageData = &image->kernel[y][uMin]; // Image values
    486483            for (int u = uMin; u <= uMax; u++, xKernelData--, imageData++) {
    … …  
    497494            float value = 0.0;    // Value of convolved pixel
    498495            int vMin = y - size, vMax = y + size; // Range for v
    499             psF32 *yKernelData = &yKernel->data.F32[yKernel->n - 1]; // Kernel values
     496            psF32 *yKernelData = &preCalc->yKernel->data.F32[yKernel->n - 1]; // Kernel values
    500497            psF32 *imageData = &temp->kernel[x][vMin]; // Image values; NOTE: wrong way!
    501498            for (int v = vMin; v <= vMax; v++, yKernelData--, imageData++) {
    … …  
    514511#else
    515512    // Convolving using precalculated kernel
    516     return p_pmSubtractionConvolveStampPrecalc(image, preCalc->data[2]);
     513    return p_pmSubtractionConvolveStampPrecalc(image, preCalc->kernel);
    517514#endif
    518515}
    … …  
    687684          return convolveStampISIS(image, kernels, index, footprint);
    688685      }
     686      case PM_SUBTRACTION_KERNEL_DECONV_HERM: {
     687          return convolveStampISIS(image, kernels, index, footprint);
     688      }
    689689      case PM_SUBTRACTION_KERNEL_RINGS: {
    690           psKernel *convolved = psKernelAlloc(-footprint, footprint,
    691                                               -footprint, footprint); // Convolved image
    692           psArray *preCalc = kernels->preCalc->data[index]; // Precalculated data
    693           psVector *uCoords = preCalc->data[0]; // u coordinates
    694           psVector *vCoords = preCalc->data[1]; // v coordinates
    695           psVector *poly = preCalc->data[2]; // Polynomial values
    696           int num = uCoords->n;         // Number of pixels
    697           psS32 *uData = uCoords->data.S32, *vData = vCoords->data.S32; // Dereference u,v coordinates
    698           psF32 *polyData = poly->data.F32; // Dereference polynomial values
     690          psKernel *convolved = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Convolved image
     691          pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[index]; // Precalculated data
     692
     693          int num = preCalc->uCoords->n;         // Number of pixels
     694          psS32 *uData = preCalc->uCoords->data.S32; // Dereference v coordinate
     695          psS32 *vData = preCalc->vCoords->data.S32; // Dereference u coordinate
     696          psF32 *polyData = preCalc->poly->data.F32; // Dereference polynomial values
    699697          psF32 **imageData = image->kernel;  // Dereference image
    700698          psF32 **convData = convolved->kernel; // Dereference convolved image

  • 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

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

    r26479 r26491  
    1 #ifndef PM_SUBTRACTION_HERMITIAN_H
    2 #define PM_SUBTRACTION_HERMITIAN_H
     1#ifndef PM_SUBTRACTION_DECONVOLVE_H
     2#define PM_SUBTRACTION_DECONVOLVE_H
    33
    44/* these function support deconvolution operations used to generate deconvolved kernels.  These
    … …  
    1111psKernel *pmSubtractionDeconvolveKernel (psKernel *kernelTarg, psKernel *kernelConv);
    1212
    13 bool pmSubtractionDeconvolutionTest ();
     13bool pmSubtractionDeconvolutionTest (int order);
    1414
    1515# endif

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

    r26486 r26491  
    1111#include "pmSubtractionKernels.h"
    1212#include "pmSubtractionHermitian.h"
     13#include "pmSubtractionDeconvolve.h"
     14#include "pmSubtractionVisual.h"
    1315
    1416#define RINGS_BUFFER 10                 // Buffer size for RINGS data
    … …  
    2729    psFree(kernels->solution1);
    2830    psFree(kernels->solution2);
     31}
     32
     33// Free function for pmSubtractionPreCalcKernel
     34static void pmSubtractionKernelPreCalcFree(pmSubtractionKernelPreCalc *kernel)
     35{
     36    psFree(kernel->xKernel);
     37    psFree(kernel->yKernel);
     38    psFree(kernel->kernel);
     39
     40    psFree(kernel->uCoords);
     41    psFree(kernel->vCoords);
     42    psFree(kernel->poly);
    2943}
    3044
    … …  
    132146}
    133147
     148bool pmSubtractionKernelPreCalcNormalize (pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc, int index, int size, int uOrder, int vOrder, float fwhm) {
     149
     150    // Calculate moments
     151    double moment = 0.0;    // Moment, for penalty
     152    for (int v = -size; v <= size; v++) {
     153        for (int u = -size; u <= size; u++) {
     154            double value = preCalc->kernel->kernel[v][u];
     155            moment += value * PS_SQR((PS_SQR(u) + PS_SQR(v)));
     156        }
     157    }
     158
     159    // Normalize sum of kernel component to unity for even functions
     160    if (uOrder % 2 == 0 && vOrder % 2 == 0) {
     161        double sum = 0.0;   // Sum of kernel component
     162        for (int v = -size; v <= size; v++) {
     163            for (int u = -size; u <= size; u++) {
     164                sum += preCalc->kernel->kernel[v][u];
     165            }
     166        }
     167        sum = 1.0 / sqrt(sum);
     168        psBinaryOp(preCalc->xKernel, preCalc->xKernel, "*", psScalarAlloc(sum, PS_TYPE_F32));
     169        psBinaryOp(preCalc->yKernel, preCalc->yKernel, "*", psScalarAlloc(sum, PS_TYPE_F32));
     170        psBinaryOp(preCalc->kernel->image, preCalc->kernel->image, "*", psScalarAlloc(PS_SQR(sum), PS_TYPE_F32));
     171
     172#if 1
     173        fprintf(stderr, "%d norm: %lf, null: %f\n", index, sum, preCalc->kernel->kernel[0][0]);
     174#endif
     175                   
     176        preCalc->kernel->kernel[0][0] -= 1.0;
     177        moment *= PS_SQR(sum);
     178    }
     179
     180#if 1
     181    double sum = 0.0;   // Sum of kernel component
     182    for (int v = -size; v <= size; v++) {
     183        for (int u = -size; u <= size; u++) {
     184            sum += preCalc->kernel->kernel[v][u];
     185        }
     186    }
     187    fprintf(stderr, "%d sum: %lf\n", index, sum);
     188#endif
     189
     190    kernels->widths->data.F32[index] = fwhm;
     191    kernels->u->data.S32[index] = uOrder;
     192    kernels->v->data.S32[index] = vOrder;
     193    if (kernels->preCalc->data[index]) {
     194        psFree(kernels->preCalc->data[index]);
     195    }
     196    kernels->preCalc->data[index] = preCalc;
     197    kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment);
     198
     199    psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index,
     200            fwhm, uOrder, vOrder, fabsf(moment));
     201
     202    return true;
     203}
     204
    134205pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, int spatialOrder,
    135206                                                    const psVector *fwhmsIN, const psVector *ordersIN,
    … …  
    172243
    173244    // Set the kernel parameters
    174     int fullSize = 2 * size + 1;        // Full size of kernels
    175245    for (int i = 0, index = 0; i < numGaussians; i++) {
    176246        float sigma = fwhms->data.F32[i] / (2.0 * sqrtf(2.0 * logf(2.0))); // Gaussian sigma
    … …  
    178248        for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
    179249            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
    180                 psArray *preCalc = psArrayAlloc(3); // Array to hold precalculated values
    181                 psVector *xKernel = preCalc->data[0] = pmSubtractionKernelISIS(sigma, uOrder, size); // x Kernel
    182                 psVector *yKernel = preCalc->data[1] = pmSubtractionKernelISIS(sigma, vOrder, size); // y Kernel
    183                 psKernel *kernel = preCalc->data[2] = psKernelAlloc(-size, size, -size, size);      // Kernel
    184 
    185                 // Calculate moments
    186                 double moment = 0.0;    // Moment, for penalty
    187                 for (int v = -size, y = 0; v <= size; v++, y++) {
    188                     for (int u = -size, x = 0; u <= size; u++, x++) {
    189                         double value = xKernel->data.F32[x] * yKernel->data.F32[y]; // Value of kernel
    190                         kernel->kernel[v][u] = value;
    191                         moment += value * PS_SQR((PS_SQR(u) + PS_SQR(v)));
    192                     }
    193                 }
    194 
    195                 // Normalise sum of kernel component to unity for even functions
    196                 if (uOrder % 2 == 0 && vOrder % 2 == 0) {
    197                     double sum = 0.0;   // Sum of kernel component
    198                     for (int v = 0; v < fullSize; v++) {
    199                         for (int u = 0; u < fullSize; u++) {
    200                             sum += xKernel->data.F32[u] * yKernel->data.F32[v];
    201                         }
    202                     }
    203                     sum = 1.0 / sqrt(sum);
    204                     psBinaryOp(xKernel, xKernel, "*", psScalarAlloc(sum, PS_TYPE_F32));
    205                     psBinaryOp(yKernel, yKernel, "*", psScalarAlloc(sum, PS_TYPE_F32));
    206                     psBinaryOp(kernel->image, kernel->image, "*", psScalarAlloc(PS_SQR(sum), PS_TYPE_F32));
    207                     kernel->kernel[0][0] -= 1.0;
    208                     moment *= PS_SQR(sum);
    209                 }
    210 
    211 
    212 #if 0
    213                 double sum = 0.0;   // Sum of kernel component
    214                 for (int v = -size; v <= size; v++) {
    215                     for (int u = -size; u <= size; u++) {
    216                         sum += kernel->kernel[v][u];
    217                     }
    218                 }
    219                 fprintf(stderr, "%d sum: %lf\n", index, sum);
    220 #endif
    221 
    222                 kernels->widths->data.F32[index] = fwhms->data.F32[i];
    223                 kernels->u->data.S32[index] = uOrder;
    224                 kernels->v->data.S32[index] = vOrder;
    225                 if (kernels->preCalc->data[index]) {
    226                     psFree(kernels->preCalc->data[index]);
    227                 }
    228                 kernels->preCalc->data[index] = preCalc;
    229                 kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment);
    230 
    231                 psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index,
    232                         fwhms->data.F32[i], uOrder, vOrder, fabsf(moment));
     250
     251                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS, uOrder, vOrder, size, sigma); // structure to hold precalculated values
     252                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i]);
    233253            }
    234254        }
    … …  
    276296
    277297    // Set the kernel parameters
    278     int fullSize = 2 * size + 1;        // Full size of kernels
    279298    for (int i = 0, index = 0; i < numGaussians; i++) {
    280299        float sigma = fwhms->data.F32[i] / (2.0 * sqrtf(2.0 * logf(2.0))); // Gaussian sigma
    … …  
    282301        for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
    283302            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
    284                 psArray *preCalc = psArrayAlloc(3); // Array to hold precalculated values
    285                 psVector *xKernel = preCalc->data[0] = pmSubtractionKernelHERM(sigma, uOrder, size); // x Kernel
    286                 psVector *yKernel = preCalc->data[1] = pmSubtractionKernelHERM(sigma, vOrder, size); // y Kernel
    287                 psKernel *kernel = preCalc->data[2] = psKernelAlloc(-size, size, -size, size);      // Kernel
    288 
    289                 // Calculate moments
    290                 double moment = 0.0;    // Moment, for penalty
    291                 for (int v = -size, y = 0; v <= size; v++, y++) {
    292                     for (int u = -size, x = 0; u <= size; u++, x++) {
    293                         double value = xKernel->data.F32[x] * yKernel->data.F32[y]; // Value of kernel
    294                         kernel->kernel[v][u] = value;
    295                         moment += value * PS_SQR((PS_SQR(u) + PS_SQR(v)));
    296                     }
    297                 }
    298 
    299                 // Normalise sum of kernel component to unity for even functions
    300                 if (uOrder % 2 == 0 && vOrder % 2 == 0) {
    301                     double sum = 0.0;   // Sum of kernel component
    302                     for (int v = 0; v < fullSize; v++) {
    303                         for (int u = 0; u < fullSize; u++) {
    304                             sum += xKernel->data.F32[u] * yKernel->data.F32[v];
    305                         }
    306                     }
    307                     sum = 1.0 / sqrt(sum);
    308                     psBinaryOp(xKernel, xKernel, "*", psScalarAlloc(sum, PS_TYPE_F32));
    309                     psBinaryOp(yKernel, yKernel, "*", psScalarAlloc(sum, PS_TYPE_F32));
    310                     psBinaryOp(kernel->image, kernel->image, "*", psScalarAlloc(PS_SQR(sum), PS_TYPE_F32));
    311 
    312 #if 1
    313                     fprintf(stderr, "%d norm: %lf, null: %f\n", index, sum,kernel->kernel[0][0]);
    314 #endif
    315                    
    316                     kernel->kernel[0][0] -= 1.0;
    317                     moment *= PS_SQR(sum);
    318                 }
    319 
    320 
    321 #if 1
    322                 double sum = 0.0;   // Sum of kernel component
    323                 for (int v = -size; v <= size; v++) {
    324                     for (int u = -size; u <= size; u++) {
    325                         sum += kernel->kernel[v][u];
    326                     }
    327                 }
    328                 fprintf(stderr, "%d sum: %lf\n", index, sum);
    329 #endif
    330 
    331                 kernels->widths->data.F32[index] = fwhms->data.F32[i];
    332                 kernels->u->data.S32[index] = uOrder;
    333                 kernels->v->data.S32[index] = vOrder;
    334                 if (kernels->preCalc->data[index]) {
    335                     psFree(kernels->preCalc->data[index]);
    336                 }
    337                 kernels->preCalc->data[index] = preCalc;
    338                 kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment);
    339 
    340                 psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index,
    341                         fwhms->data.F32[i], uOrder, vOrder, fabsf(moment));
     303                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_HERM, uOrder, vOrder, size, sigma); // structure to hold precalculated values
     304                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i]);
    342305            }
    343306        }
    … …  
    347310}
    348311
    349 # if (0)
    350 pmSubtractionKernels *pmSubtractionKernelsDECON_HERM(int size, int spatialOrder,
     312pmSubtractionKernels *pmSubtractionKernelsDECONV_HERM(int size, int spatialOrder,
    351313                                                     const psVector *fwhmsIN, const psVector *ordersIN,
    352314                                                     float penalty, pmSubtractionMode mode)
    … …  
    379341    }
    380342
    381     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_DECON_HERM, size, spatialOrder, penalty, mode); // The kernels
    382     psStringAppend(&kernels->description, "DECON_HERM(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
     343    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_DECONV_HERM, size, spatialOrder, penalty, mode); // The kernels
     344    psStringAppend(&kernels->description, "DECONV_HERM(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
    383345
    384346    psLogMsg("psModules.imcombine", PS_LOG_INFO, "DECONVOLVED HERM kernel: %s,%d --> %d elements", params, spatialOrder, num);
    385347    psFree(params);
    386348
    387     // XXXXX hard-wired reference sigma for now of 1.7 pix
     349    // XXXXX hard-wired reference sigma for now of 1.7 pix (== 4.0 pix fwhm == 1.0 arcsec in simtest)
    388350    // generate the Gaussian deconvolution kernel
    389     # define DECON_SIGMA 1.7
    390     psKernel *kernelGauss = pmSubtractionDeconvolveGauss (size, DECON_SIGMA);
     351    # define DECONV_SIGMA 1.7
     352    psKernel *kernelGauss = pmSubtractionDeconvolveGauss (size, DECONV_SIGMA);
    391353
    392354    // Set the kernel parameters
    393     int fullSize = 2 * size + 1;        // Full size of kernels
    394355    for (int i = 0, index = 0; i < numGaussians; i++) {
    395356        float sigma = fwhms->data.F32[i] / (2.0 * sqrtf(2.0 * logf(2.0))); // Gaussian sigma
    … …  
    397358        for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
    398359            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
    399                 psArray *preCalc  = psArrayAlloc(3); // Array to hold precalculated values
    400                 psVector *xKernel = preCalc->data[0] = pmSubtractionKernelHERM(sigma, uOrder, size); // x Kernel
    401                 psVector *yKernel = preCalc->data[1] = pmSubtractionKernelHERM(sigma, vOrder, size); // y Kernel
    402                 psKernel *kernelTarget = psKernelAlloc(-size, size, -size, size);       // Kernel
    403 
    404                 // generate 2D kernel, calculate moments
    405                 for (int v = -size, y = 0; v <= size; v++, y++) {
    406                     for (int u = -size, x = 0; u <= size; u++, x++) {
    407                         double value = xKernel->data.F32[x] * yKernel->data.F32[y]; // Value of kernel
    408                         kernelTarget->kernel[v][u] = value;
    409                     }
    410                 }
    411 
    412                 // deconvolve the target by the gaussian:
    413                 psKernel *kernel = pmSubtractionDeconvolveKernel(kernelTarget, kernelGauss); // Kernel
    414                 preCalc->data[2] = kernel;
    415 
    416                 psImage *kernelConv = psImageConvolveFFT(NULL, kernel, NULL, 0, kernelGauss);
    417                 pmSubtractionVisualShowSubtraction (kernelTarget->image, kernel->image, kernelConv->image);
    418 
    419                 // Normalise sum of kernel component to unity for even functions
    420                 if (uOrder % 2 == 0 && vOrder % 2 == 0) {
    421                     double sum = 0.0;   // Sum of kernel component
    422                     for (int v = 0; v < fullSize; v++) {
    423                         for (int u = 0; u < fullSize; u++) {
    424                             sum += xKernel->data.F32[u] * yKernel->data.F32[v];
    425                         }
    426                     }
    427                     sum = 1.0 / sqrt(sum);
    428                     psBinaryOp(xKernel, xKernel, "*", psScalarAlloc(sum, PS_TYPE_F32));
    429                     psBinaryOp(yKernel, yKernel, "*", psScalarAlloc(sum, PS_TYPE_F32));
    430                     psBinaryOp(kernel->image, kernel->image, "*", psScalarAlloc(PS_SQR(sum), PS_TYPE_F32));
    431 
    432 #if 1
    433                     fprintf(stderr, "%d norm: %e, null: %e\n", index, sum, kernel->kernel[0][0]);
    434 #endif
    435                    
    436                     kernel->kernel[0][0] -= 1.0;
    437                     moment *= PS_SQR(sum);
    438                 }
    439 
    440 
    441 #if 1
    442                 double sum = 0.0;   // Sum of kernel component
    443                 for (int v = -size; v <= size; v++) {
    444                     for (int u = -size; u <= size; u++) {
    445                         sum += kernel->kernel[v][u];
    446                     }
    447                 }
    448                 fprintf(stderr, "%d sum: %e\n", index, sum);
    449 #endif
    450 
    451                 kernels->widths->data.F32[index] = fwhms->data.F32[i];
    452                 kernels->u->data.S32[index] = uOrder;
    453                 kernels->v->data.S32[index] = vOrder;
    454                 if (kernels->preCalc->data[index]) {
    455                     psFree(kernels->preCalc->data[index]);
    456                 }
    457                 kernels->preCalc->data[index] = preCalc;
    458                 kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment);
    459 
    460                 psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index,
    461                         fwhms->data.F32[i], uOrder, vOrder, fabsf(moment));
     360
     361                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_HERM, uOrder, vOrder, size, sigma); // structure to hold precalculated values
     362
     363                // save the generated 2D kernel as the target, deconvolve it by Gaussian, replacing the generated 2D kernel
     364                psKernel *kernelTarget = preCalc->kernel;
     365                preCalc->kernel = pmSubtractionDeconvolveKernel(kernelTarget, kernelGauss); // Kernel
     366
     367                // XXXX test demo that deconvolved kernel is valid
     368                // psImage *kernelConv = psImageConvolveFFT(NULL, preCalc->kernel->image, NULL, 0, kernelGauss);
     369                // pmSubtractionVisualShowSubtraction (kernelTarget->image, preCalc->kernel->image, kernelConv);
     370
     371                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i]);
    462372            }
    463373        }
    … …  
    466376    return kernels;
    467377}
    468 # endif
    469378
    470379//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    … …  
    503412}
    504413
     414pmSubtractionKernelPreCalc *pmSubtractionKernelPreCalcAlloc(pmSubtractionKernelsType type, int uOrder, int vOrder, int size, float sigma) {
     415
     416    pmSubtractionKernelPreCalc *preCalc = psAlloc(sizeof(pmSubtractionKernelPreCalc)); // Kernels, to return
     417    psMemSetDeallocator(preCalc, (psFreeFunc)pmSubtractionKernelPreCalcFree);
     418
     419    // 1D kernel realizations:
     420    switch (type) {
     421      case PM_SUBTRACTION_KERNEL_ISIS:
     422        preCalc->xKernel = pmSubtractionKernelISIS(sigma, uOrder, size);
     423        preCalc->yKernel = pmSubtractionKernelISIS(sigma, vOrder, size);
     424        preCalc->uCoords = NULL;
     425        preCalc->vCoords = NULL;
     426        preCalc->poly    = NULL;
     427        break;
     428      case PM_SUBTRACTION_KERNEL_HERM:
     429        preCalc->xKernel = pmSubtractionKernelHERM(sigma, uOrder, size);
     430        preCalc->yKernel = pmSubtractionKernelHERM(sigma, vOrder, size);
     431        preCalc->uCoords = NULL;
     432        preCalc->vCoords = NULL;
     433        preCalc->poly    = NULL;
     434        break;
     435      case PM_SUBTRACTION_KERNEL_DECONV_HERM:
     436        preCalc->xKernel = pmSubtractionKernelHERM(sigma, uOrder, size);
     437        preCalc->yKernel = pmSubtractionKernelHERM(sigma, vOrder, size);
     438        preCalc->uCoords = NULL;
     439        preCalc->vCoords = NULL;
     440        preCalc->poly    = NULL;
     441        break;
     442      case PM_SUBTRACTION_KERNEL_RINGS:
     443        // the RINGS kernel uses the uCoords, vCoords, and poly elements of the structure
     444        // we allocate these vectors here, but leave the kernel generation to the main function
     445        preCalc->xKernel = NULL;
     446        preCalc->yKernel = NULL;
     447        preCalc->kernel  = NULL;
     448        preCalc->uCoords = psVectorAllocEmpty(size, PS_TYPE_S32); // u coords
     449        preCalc->vCoords = psVectorAllocEmpty(size, PS_TYPE_S32); // v coords
     450        preCalc->poly    = psVectorAllocEmpty(size, PS_TYPE_F32); // Polynomial
     451        return preCalc;
     452      default:
     453        psAbort("programming error: invalid type for PreCalc kernel");
     454    }
     455
     456    preCalc->kernel = psKernelAlloc(-size, size, -size, size); // 2D Kernel
     457
     458    // generate 2D kernel from 1D realizations
     459    for (int v = -size, y = 0; v <= size; v++, y++) {
     460        for (int u = -size, x = 0; u <= size; u++, x++) {
     461            preCalc->kernel->kernel[v][u] = preCalc->xKernel->data.F32[x] * preCalc->yKernel->data.F32[y]; // Value of kernel
     462        }
     463    }
     464   
     465    return preCalc;
     466}
     467
    505468pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, float penalty,
    506469                                               pmSubtractionMode mode)
    … …  
    827790            for (int vOrder = 0; vOrder <= (i == 0 ? 0 : ringsOrder - uOrder); vOrder++, index++) {
    828791
    829                 psArray *data = psArrayAlloc(3); // Container for data
    830                 psVector *uCoords = data->data[0] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_S32); // u coords
    831                 psVector *vCoords = data->data[1] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_S32); // v coords
    832                 psVector *poly = data->data[2] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_F32); // Polynomial
     792                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc (PM_SUBTRACTION_KERNEL_RINGS, 0, 0, RINGS_BUFFER, 0.0);
    833793                double moment = 0.0;    // Moment, for penalty
    834794
    835795                if (i == 0) {
    836796                    // Central pixel is easy
    837                     uCoords->data.S32[0] = vCoords->data.S32[0] = 0;
    838                     poly->data.F32[0] = 1.0;
    839                     uCoords->n = vCoords->n = poly->n = 1;
     797                    preCalc->uCoords->data.S32[0] = 0;
     798                    preCalc->vCoords->data.S32[0] = 0;
     799                    preCalc->poly->data.F32[0] = 1.0;
     800                    preCalc->uCoords->n = 1;
     801                    preCalc->vCoords->n = 1;
     802                    preCalc->poly->n = 1;
    840803                    radiusLast = 0;
    841804                    moment = 0.0;
    … …  
    855818                                float polyVal = uPoly * vPoly; // Value of polynomial
    856819                                if (polyVal != 0) { // No point adding it otherwise
    857                                     uCoords->data.S32[j] = u;
    858                                     vCoords->data.S32[j] = v;
    859                                     poly->data.F32[j] = polyVal;
     820                                    preCalc->uCoords->data.S32[j] = u;
     821                                    preCalc->vCoords->data.S32[j] = v;
     822                                    preCalc->poly->data.F32[j] = polyVal;
    860823                                    norm += polyVal;
    861824                                    moment += polyVal * PS_SQR(PS_SQR(u) + PS_SQR(v));
    862825
    863                                     psVectorExtend(uCoords, RINGS_BUFFER, 1);
    864                                     psVectorExtend(vCoords, RINGS_BUFFER, 1);
    865                                     psVectorExtend(poly, RINGS_BUFFER, 1);
     826                                    psVectorExtend(preCalc->uCoords, RINGS_BUFFER, 1);
     827                                    psVectorExtend(preCalc->vCoords, RINGS_BUFFER, 1);
     828                                    psVectorExtend(preCalc->poly, RINGS_BUFFER, 1);
    866829                                    psTrace("psModules.imcombine", 9, "u = %d, v = %d, poly = %f\n",
    867                                             u, v, poly->data.F32[j]);
     830                                            u, v, preCalc->poly->data.F32[j]);
    868831                                    j++;
    869832                                }
    … …  
    873836                    // Normalise kernel component to unit sum
    874837                    if (uOrder % 2 == 0 && vOrder % 2 == 0) {
    875                         psBinaryOp(poly, poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32));
     838                        psBinaryOp(preCalc->poly, preCalc->poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32));
    876839                        // Add subtraction of 0,0 component to preserve photometric scaling
    877                         uCoords->data.S32[j] = 0;
    878                         vCoords->data.S32[j] = 0;
    879                         poly->data.F32[j] = -1.0;
    880                         psVectorExtend(uCoords, RINGS_BUFFER, 1);
    881                         psVectorExtend(vCoords, RINGS_BUFFER, 1);
    882                         psVectorExtend(poly, RINGS_BUFFER, 1);
     840                        preCalc->uCoords->data.S32[j] = 0;
     841                        preCalc->vCoords->data.S32[j] = 0;
     842                        preCalc->poly->data.F32[j] = -1.0;
     843                        psVectorExtend(preCalc->uCoords, RINGS_BUFFER, 1);
     844                        psVectorExtend(preCalc->vCoords, RINGS_BUFFER, 1);
     845                        psVectorExtend(preCalc->poly, RINGS_BUFFER, 1);
    883846                    } else {
    884847                        norm = powf(size, uOrder) * powf(size, vOrder);
    885                         psBinaryOp(poly, poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32));
     848                        psBinaryOp(preCalc->poly, preCalc->poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32));
    886849                    }
    887850                    moment /= norm;
    888851                }
    889852
    890                 psTrace("psModules.imcombine", 8, "%ld pixels in kernel\n", uCoords->n);
    891 
    892                 kernels->preCalc->data[index] = data;
     853                psTrace("psModules.imcombine", 8, "%ld pixels in kernel\n", preCalc->uCoords->n);
     854
     855                kernels->preCalc->data[index] = preCalc;
    893856                kernels->u->data.S32[index] = uOrder;
    894857                kernels->v->data.S32[index] = vOrder;
    … …  
    916879      case PM_SUBTRACTION_KERNEL_HERM:
    917880        return pmSubtractionKernelsHERM(size, spatialOrder, fwhms, orders, penalty, mode);
     881      case PM_SUBTRACTION_KERNEL_DECONV_HERM:
     882        return pmSubtractionKernelsDECONV_HERM(size, spatialOrder, fwhms, orders, penalty, mode);
    918883      case PM_SUBTRACTION_KERNEL_SPAM:
    919884        return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, penalty, mode);
    … …  
    978943    float penalty = 0.0;                // Penalty for wideness
    979944
    980     // ISIS and HERM have the same description layout
    981     if (!strncmp(description, "ISIS", 4) || !strcmp (description, "HERM")) {
     945    // ISIS, HERM, and DECONV_HERM have the same description layout
     946    if (!strncmp(description, "ISIS", 4) || !strcmp (description, "HERM") || !strcmp (description, "DECONV_HERM")) {
    982947        // XXX Support for GUNK (not yet supported)
    983948        if (strstr(description, "+POIS")) {
    … …  
    987952
    988953        type = pmSubtractionKernelsTypeFromString (description);
    989         psAssert (type != PM_SUBTRACTION_KERNEL_NONE, "must  be ISIS or HERM");
    990 
    991         char *ptr = (char*)description + 5;    // Eat "ISIS(" or "HERM("
     954        psAssert (type != PM_SUBTRACTION_KERNEL_NONE, "must  be ISIS, HERM or DECONV_HERM");
     955
     956        char *ptr = NULL;
     957        switch (type) {
     958          case PM_SUBTRACTION_KERNEL_ISIS:
     959          case PM_SUBTRACTION_KERNEL_HERM:
     960            ptr = (char*) description + 5;    // Eat "ISIS(" or "HERM("
     961            break;
     962          case PM_SUBTRACTION_KERNEL_DECONV_HERM:
     963            ptr = (char*) description + 12;    // Eat "DECONV_HERM("
     964            break;
     965          default:
     966            psAbort("programming error: invalid kernel type");
     967        }
    992968        PARSE_STRING_NUMBER(size, ptr, ',', parseStringInt);
    993969
    … …  
    10251001    }
    10261002
    1027     psAbort("Deciphering kernels other than ISIS, HERM and RINGS is not currently supported.");
     1003    psAbort("Deciphering kernels other than ISIS, HERM, DECONV_HERM or RINGS is not currently supported.");
    10281004
    10291005    return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders,
    … …  
    10421018    if (strcasecmp(type, "HERM") == 0) {
    10431019        return PM_SUBTRACTION_KERNEL_HERM;
     1020    }
     1021    if (strcasecmp(type, "DECONV_HERM") == 0) {
     1022        return PM_SUBTRACTION_KERNEL_DECONV_HERM;
    10441023    }
    10451024    if (strcasecmp(type, "SPAM") == 0) {

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

    r26486 r26491  
    1111    PM_SUBTRACTION_KERNEL_ISIS,         ///< Traditional kernel --- gaussians modified by polynomials
    1212    PM_SUBTRACTION_KERNEL_HERM,         ///< Hermitian polynomial kernels
    13     PM_SUBTRACTION_KERNEL_DECON_HERM,   ///< Deconvolved Hermitian polynomial kernels
     13    PM_SUBTRACTION_KERNEL_DECONV_HERM,  ///< Deconvolved Hermitian polynomial kernels
    1414    PM_SUBTRACTION_KERNEL_SPAM,         ///< Summed Pixels for Advanced Matching --- summed delta functions
    1515    PM_SUBTRACTION_KERNEL_FRIES,        ///< Fibonacci Radius Increases Excellence of Subtraction
    … …  
    3232    psString description;               ///< Description of the kernel parameters
    3333    long num;                           ///< Number of kernel components (not including the spatial ones)
    34     psVector *u, *v;                    ///< Offset (for POIS) or polynomial order (for ISIS, HERM or DECON_HERM)
    35     psVector *widths;                   ///< Gaussian FWHMs (ISIS, HERM or DECON_HERM)
     34    psVector *u, *v;                    ///< Offset (for POIS) or polynomial order (for ISIS, HERM or DECONV_HERM)
     35    psVector *widths;                   ///< Gaussian FWHMs (ISIS, HERM or DECONV_HERM)
    3636    psVector *uStop, *vStop;            ///< Width of kernel element (SPAM,FRIES only)
    37     psArray *preCalc;                   ///< Array of images containing pre-calculated kernel (for ISIS, HERM or DECON_HERM)
     37    psArray *preCalc;                   ///< Array of images containing pre-calculated kernel (for ISIS, HERM or DECONV_HERM)
    3838    float penalty;                      ///< Penalty for wideness
    3939    psVector *penalties;                ///< Penalty for each kernel component
    … …  
    4949    int numStamps;                      ///< Number of good stamps
    5050} pmSubtractionKernels;
     51
     52// pmSubtractionKernels->preCalc is an array of pmSubtractionKernelPreCalc structures
     53typedef struct {
     54    psVector *uCoords;                  // used by RINGS
     55    psVector *vCoords;                  // used by RINGS
     56    psVector *poly;                     // used by RINGS
     57
     58    psVector *xKernel;                  // used by ISIS, HERM, DECONV_HERM
     59    psVector *yKernel;                  // used by ISIS, HERM, DECONV_HERM
     60    psKernel *kernel;                   // used by ISIS, HERM, DECONV_HERM
     61} pmSubtractionKernelPreCalc;
    5162
    5263// Assertion to check pmSubtractionKernels
    … …  
    7182        PS_ASSERT_VECTOR_SIZE((KERNELS)->widths, (KERNELS)->num, RETURNVALUE); \
    7283    } \
    73     if ((KERNELS)->type == PM_SUBTRACTION_KERNEL_DECON_HERM) { \
     84    if ((KERNELS)->type == PM_SUBTRACTION_KERNEL_DECONV_HERM) { \
    7485        PS_ASSERT_VECTOR_NON_NULL((KERNELS)->widths, RETURNVALUE); \
    7586        PS_ASSERT_VECTOR_TYPE((KERNELS)->widths, PS_TYPE_F32, RETURNVALUE); \
    … …  
    139150    );
    140151
     152/// Allocator for pre-calculated kernel data structure
     153pmSubtractionKernelPreCalc *pmSubtractionKernelPreCalcAlloc(
     154    pmSubtractionKernelsType type, ///< type of kernel to allocate (not all can be pre-calculated)
     155    int uOrder,                    ///< order in x-direction
     156    int vOrder,                    ///< order in x-direction
     157    int size,                      ///< Half-size of the kernel
     158    float sigma                    ///< sigma of gaussian kernel
     159    );
     160
     161
    141162/// Generate POIS kernels
    142163pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, ///< Half-size of the kernel (in both dims)
    … …  
    173194                                               );
    174195
    175 /// Generate DECON_HERM kernels
    176 pmSubtractionKernels *pmSubtractionKernelsDECON_HERM(int size, ///< Half-size of the kernel
     196/// Generate DECONV_HERM kernels
     197pmSubtractionKernels *pmSubtractionKernelsDECONV_HERM(int size, ///< Half-size of the kernel
    177198                                                     int spatialOrder, ///< Order of spatial variations
    178199                                                     const psVector *fwhms, ///< Gaussian FWHMs

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

    r26035 r26491  
    482482    // Maintain photometric scaling
    483483    if (type == PM_SUBTRACTION_KERNEL_ISIS) {
    484         psKernel *subtract = kernels->preCalc->data[0]; // Kernel to subtract from the rest
     484
     485        // XXX in r26035, this code was just wrong.  we had:
     486
     487        // psKernel *subtract = kernels->preCalc->data[0]
     488
     489        // but, kernels->preCalc was an array of psArray, not an array of kernels.  It is now
     490        // an array of pmSubtractionKernelPreCalc.
     491
     492        pmSubtractionKernelPreCalc *subtract = kernels->preCalc->data[0]; // Kernel to subtract from the rest
     493
    485494        for (int i = 1; i < newSize; i++) {
    486495            if (kernels->u->data.S32[i] % 2 == 0 && kernels->v->data.S32[i] % 2 == 0) {
    487                 psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest
    488                 psBinaryOp(kernel->image, kernel->image, "-", subtract->image);
     496                pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[i]; // Kernel of interest
     497                psBinaryOp(preCalc->kernel->image, preCalc->kernel->image, "-", subtract->kernel->image);
    489498            }
    490499        }
    … …  
    495504        for (int i = 0; i < newSize; i++) {
    496505            if (kernels->u->data.S32[i] % 2 == 0 && kernels->v->data.S32[i] % 2 == 0) {
    497                 psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest
    498                 kernel->kernel[0][0] -= 1.0;
     506                pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[i]; // Kernel of interest
     507                preCalc->kernel->kernel[0][0] -= 1.0;
    499508            }
    500509        }

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

    r26472 r26491  
    242242
    243243    for (int i = 0; i < kernels->num; i++) {
    244         psArray *preCalc = kernels->preCalc->data[i];
    245         psKernel *kernel = preCalc->data[2];
     244        pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[i];
     245        psKernel *kernel = preCalc->kernel;
    246246
    247247        int xSub = i % NXsub;
Note: See TracChangeset for help on using the changeset viewer.