IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 29597


Ignore:
Timestamp:
Oct 28, 2010, 5:49:12 PM (16 years ago)
Author:
eugene
Message:

compile-time options to (a) use a central delta function for kernel zeroing; (b) use standard AL kernel zeroing; (c) use a smaller gaussian for kernel zeroing; add function pmSubtractionKernelsImageMosaic

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r29543 r29597  
    183183}
    184184
    185 # define CENTRAL_DELTA 0
    186 
    187185# if (CENTRAL_DELTA)
    188 
    189 // XXX *** this code used the central pixel to force zero net flux,
    190 // Alard actually uses kernel(0) for this purpose (for even order, kernel[i] = kernel'[i] - kernel[0])
     186// This version of the code uses the central pixel to force zero net flux, Alard actually uses
     187// kernel(0) for this purpose (for even order, kernel[i] = kernel'[i] - kernel[0])
    191188static bool pmSubtractionKernelPreCalcNormalize(pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc,
    192189                                                int index, int uOrder, int vOrder, float fwhm,
    … …  
    301298static bool pmSubtractionKernelPreCalcNormalize(pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc,
    302299                                                int index, int uOrder, int vOrder, float fwhm,
    303                                                 bool AlardLuptonStyle, bool forceZeroNull)
     300                                                pmSubtractionKernelPreCalc *zeroKernel)
    304301{
    305302    // 1) for odd functions: no renormalization
    … …  
    339336            }
    340337        }
    341         if (index == 0) {
     338# if (ZERO_KERNEL_ZERO_FLUX)
     339        int firstZeroIndex = 0;
     340# else
     341        int firstZeroIndex = 1;
     342# endif
     343        if (index < firstZeroIndex) {
    342344            zeroIsNormal = true;
    343345        } else {
    344346            psAssert(zeroIsNormal, "failed to normalize zero kernel first");
    345             pmSubtractionKernelPreCalc *zeroKernel = kernels->preCalc->data[0];
    346347            psAssert(zeroKernel, "failed to supply zero kernel");
    347348            for (int v = preCalc->kernel->yMin; v <= preCalc->kernel->yMax; v++) {
    … …  
    379380#endif
    380381
    381     kernels->widths->data.F32[index] = fwhm;
    382     kernels->u->data.S32[index] = uOrder;
    383     kernels->v->data.S32[index] = vOrder;
    384     if (kernels->preCalc->data[index]) {
    385         psFree(kernels->preCalc->data[index]);
    386     }
    387     kernels->preCalc->data[index] = preCalc;
     382    if (kernels) {
     383        kernels->widths->data.F32[index] = fwhm;
     384        kernels->u->data.S32[index] = uOrder;
     385        kernels->v->data.S32[index] = vOrder;
     386        if (kernels->preCalc->data[index]) {
     387            psFree(kernels->preCalc->data[index]);
     388        }
     389        kernels->preCalc->data[index] = preCalc;
     390    }
    388391    psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index, fwhm, uOrder, vOrder);
    389392
    … …  
    431434    psFree(params);
    432435
     436# if (!CENTRAL_DELTA && !ZERO_KERNEL_ZERO_FLUX)
     437    // in this case, subtract the 0th kernel from everyone else
     438    float zeroFWHM = fwhms->data.F32[0];
     439    float zeroSigma = zeroFWHM / (2.0 * sqrtf(2.0 * logf(2.0))); // Gaussian sigma
     440    pmSubtractionKernelPreCalc *zeroKernel = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS, 0, 0, size, zeroSigma); // structure to hold precalculated values
     441    pmSubtractionKernelPreCalcNormalize (NULL, zeroKernel, -1, 0, 0, zeroFWHM, NULL);
     442# endif
     443
    433444    // Set the kernel parameters
    434445    for (int i = 0, index = 0; i < numGaussians; i++) {
    435446        float sigma = fwhms->data.F32[i] / (2.0 * sqrtf(2.0 * logf(2.0))); // Gaussian sigma
     447
     448# if (!CENTRAL_DELTA && ZERO_KERNEL_ZERO_FLUX)
     449        // in this case, subtract the a 1/2 width Gaussian from each series
     450        float zeroFWHM = 0.50 * fwhms->data.F32[0];
     451        float zeroSigma = zeroFWHM / (2.0 * sqrtf(2.0 * logf(2.0))); // Gaussian sigma
     452        pmSubtractionKernelPreCalc *zeroKernel = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS, 0, 0, size, zeroSigma); // structure to hold precalculated values
     453        pmSubtractionKernelPreCalcNormalize (NULL, zeroKernel, -1, 0, 0, zeroFWHM, NULL);
     454# endif
    436455        // Iterate over (u,v) order
    437456        for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
    … …  
    439458
    440459                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS, uOrder, vOrder, size, sigma); // structure to hold precalculated values
    441                 pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false);
    442                 // pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], false, false);
     460# if (CENTRAL_DELTA)
     461                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], false, false);
     462# else
     463                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], zeroKernel);
     464# endif
    443465            }
    444466        }
    445     }
    446 
     467# if (!CENTRAL_DELTA && ZERO_KERNEL_ZERO_FLUX)
     468        psFree(zeroKernel);
     469# endif
     470    }
     471
     472# if (!CENTRAL_DELTA && !ZERO_KERNEL_ZERO_FLUX)
     473    psFree(zeroKernel);
     474# endif
    447475    psFree(orders);
    448476    psFree(fwhms);
    … …  
    497525            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
    498526                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS, uOrder, vOrder, size, sigma); // structure to hold precalculated values
     527# if (CENTRAL_DELTA)
    499528                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false);
     529# else
     530                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], NULL);
     531# endif
    500532            }
    501533        }
    … …  
    503535            // XXX modify size for hermitians to account for sqrt(2) in Hermitian definition (relative to ISIS Gaussian)
    504536            pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS_RADIAL, order, order, size, sigma / sqrt(2.0)); // structure to hold precalculated values
     537# if (CENTRAL_DELTA)
    505538            pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, order, order, fwhms->data.F32[i], true, true);
     539# else
     540            pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, order, order, fwhms->data.F32[i], NULL);
     541# endif
    506542        }
    507543    }
    … …  
    555591            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
    556592                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_HERM, uOrder, vOrder, size, sigma); // structure to hold precalculated values
     593# if (CENTRAL_DELTA)
    557594                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false);
     595# else
     596                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], NULL);
     597# endif
    558598            }
    559599        }
    … …  
    624664
    625665                // XXX do we use Alard-Lupton normalization (last param true) or not?
     666# if (CENTRAL_DELTA)
    626667                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false);
    627 
     668# else         
     669                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], NULL);
     670# endif
    628671                // XXXX test demo that deconvolved kernel is valid
    629672# if 1
    … …  
    13961439    return out;
    13971440}
     1441
     1442#define KERNEL_MOSAIC 2                 // Half-number of kernel instances in the mosaic image
     1443psImage *pmSubtractionKernelsImageMosaic(pmSubtractionKernels *kernels) {
     1444
     1445    psTrace("psModules.imcombine", 2, "Generating diagnostic image..\n");
     1446
     1447    // Generate image with convolution kernels
     1448    int size = kernels->size;       // Half-size of kernel
     1449    int fullSize = 2 * size + 1 + 1; // Full size of kernel
     1450    int imageSize = (2 * KERNEL_MOSAIC + 1) * fullSize;
     1451    psImage *convKernels = psImageAlloc((kernels->mode == PM_SUBTRACTION_MODE_DUAL ? 2 : 1) *
     1452                                        imageSize - 1 +
     1453                                        (kernels->mode == PM_SUBTRACTION_MODE_DUAL ? 4 : 0),
     1454                                        imageSize - 1, PS_TYPE_F32);
     1455    psImageInit(convKernels, NAN);
     1456    for (int j = -KERNEL_MOSAIC; j <= KERNEL_MOSAIC; j++) {
     1457        for (int i = -KERNEL_MOSAIC; i <= KERNEL_MOSAIC; i++) {
     1458            psImage *kernel = pmSubtractionKernelImage(kernels, (float)i / (float)KERNEL_MOSAIC,
     1459                                                       (float)j / (float)KERNEL_MOSAIC,
     1460                                                       false); // Image of the kernel
     1461            if (!kernel) {
     1462                psError(psErrorCodeLast(), false, "Unable to generate kernel image.");
     1463                psFree(convKernels);
     1464                return NULL;
     1465            }
     1466
     1467            if (psImageOverlaySection(convKernels, kernel, (i + KERNEL_MOSAIC) * fullSize,
     1468                                      (j + KERNEL_MOSAIC) * fullSize, "=") == 0) {
     1469                psError(psErrorCodeLast(), false, "Unable to overlay kernel image.");
     1470                psFree(kernel);
     1471                psFree(convKernels);
     1472                return NULL;
     1473            }
     1474            psFree(kernel);
     1475
     1476            if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
     1477                kernel = pmSubtractionKernelImage(kernels, (float)i / (float)KERNEL_MOSAIC,
     1478                                                  (float)j / (float)KERNEL_MOSAIC,
     1479                                                  true); // Image of the kernel
     1480                if (!kernel) {
     1481                    psError(psErrorCodeLast(), false, "Unable to generate kernel image.");
     1482                    psFree(convKernels);
     1483                    return NULL;
     1484                }
     1485
     1486                if (psImageOverlaySection(convKernels, kernel,
     1487                                          (2 * KERNEL_MOSAIC + 1 + i + KERNEL_MOSAIC) * fullSize + 4,
     1488                                          (j + KERNEL_MOSAIC) * fullSize, "=") == 0) {
     1489                    psError(psErrorCodeLast(), false, "Unable to overlay kernel image.");
     1490                    psFree(kernel);
     1491                    psFree(convKernels);
     1492                    return NULL;
     1493                }
     1494                psFree(kernel);
     1495            }
     1496        }
     1497    }
     1498    return convKernels;
     1499}
Note: See TracChangeset for help on using the changeset viewer.