IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 16, 2007, 10:32:38 AM (19 years ago)
Author:
Paul Price
Message:

Using FWHM instead of Gaussian sigmas. I figure FWHM is more intuitive and common for astronomers.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r14525 r14532  
    9999
    100100pmSubtractionKernels *pmSubtractionKernelsISIS(int size, int spatialOrder,
    101                                                const psVector *sigmas, const psVector *orders)
    102 {
    103     PS_ASSERT_VECTOR_NON_NULL(sigmas, NULL);
    104     PS_ASSERT_VECTOR_TYPE(sigmas, PS_TYPE_F32, NULL);
     101                                               const psVector *fwhms, const psVector *orders)
     102{
     103    PS_ASSERT_VECTOR_NON_NULL(fwhms, NULL);
     104    PS_ASSERT_VECTOR_TYPE(fwhms, PS_TYPE_F32, NULL);
    105105    PS_ASSERT_VECTOR_NON_NULL(orders, NULL);
    106106    PS_ASSERT_VECTOR_TYPE(orders, PS_TYPE_S32, NULL);
    107     PS_ASSERT_VECTORS_SIZE_EQUAL(sigmas, orders, NULL);
     107    PS_ASSERT_VECTORS_SIZE_EQUAL(fwhms, orders, NULL);
    108108    PS_ASSERT_INT_POSITIVE(size, NULL);
    109109    PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL);
    110110
    111     int numGaussians = sigmas->n;       // Number of Gaussians
     111    int numGaussians = fwhms->n;       // Number of Gaussians
    112112
    113113    int num = 0;                        // Number of basis functions
    … …  
    115115    for (int i = 0; i < numGaussians; i++) {
    116116        int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian
    117         psStringAppend(&params, "(%.2f,%d)", sigmas->data.F32[i], orders->data.S32[i]);
     117        psStringAppend(&params, "(%.2f,%d)", fwhms->data.F32[i], orders->data.S32[i]);
    118118        num += (gaussOrder + 1) * (gaussOrder + 2) / 2;
    119119    }
    … …  
    133133    // Set the kernel parameters
    134134    for (int i = 0, index = 0; i < numGaussians; i++) {
    135         float norm = 1.0 / (M_2_PI * sqrtf(sigmas->data.F32[i])); // Normalisation for Gaussian
     135        float sigma = fwhms->data.F32[i] / (2.0 * sqrtf(2.0 * logf(2.0))); // Gaussian sigma
     136        float norm = 1.0 / (M_2_PI * sqrtf(sigma)); // Normalisation for Gaussian
     137
    136138        // Iterate over (u,v) order
    137139        for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
    … …  
    144146                    for (int u = -size; u <= size; u++) {
    145147                        sum += preCalc->kernel[v][u] = norm * power(u, uOrder) * power(v, vOrder) *
    146                             expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigmas->data.F32[i]));
     148                            expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigma));
    147149                    }
    148150                }
    … …  
    163165                }
    164166
    165                 kernels->widths->data.F32[index] = sigmas->data.F32[i];
     167                kernels->widths->data.F32[index] = fwhms->data.F32[i];
    166168                kernels->u->data.S32[index] = uOrder;
    167169                kernels->v->data.S32[index] = vOrder;
    … …  
    169171
    170172                psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index,
    171                         sigmas->data.F32[i], uOrder, vOrder);
     173                        fwhms->data.F32[i], uOrder, vOrder);
    172174            }
    173175        }
    … …  
    378380
    379381// Grid United with Normal Kernel
    380 pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, int spatialOrder, const psVector *sigmas,
     382pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, int spatialOrder, const psVector *fwhms,
    381383                                               const psVector *orders, int inner)
    382384{
    383385    PS_ASSERT_INT_POSITIVE(size, NULL);
    384386    PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL);
    385     PS_ASSERT_VECTOR_NON_NULL(sigmas, NULL);
    386     PS_ASSERT_VECTOR_TYPE(sigmas, PS_TYPE_F32, NULL);
     387    PS_ASSERT_VECTOR_NON_NULL(fwhms, NULL);
     388    PS_ASSERT_VECTOR_TYPE(fwhms, PS_TYPE_F32, NULL);
    387389    PS_ASSERT_VECTOR_NON_NULL(orders, NULL);
    388390    PS_ASSERT_VECTOR_TYPE(orders, PS_TYPE_S32, NULL);
    389     PS_ASSERT_VECTORS_SIZE_EQUAL(sigmas, orders, NULL);
     391    PS_ASSERT_VECTORS_SIZE_EQUAL(fwhms, orders, NULL);
    390392    PS_ASSERT_INT_NONNEGATIVE(inner, NULL);
    391393    PS_ASSERT_INT_LESS_THAN(inner, size, NULL);
    392394
    393     int numGaussians = sigmas->n;       // Number of Gaussians
     395    int numGaussians = fwhms->n;       // Number of Gaussians
    394396    int numGaussianVars = 0;            // Number of Gaussian variant functions in the kernel
    395397    psString params = NULL;             // List of params
    … …  
    397399        int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian
    398400        numGaussianVars += (gaussOrder + 1) * (gaussOrder + 2) / 2;
    399         psStringAppend(&params, "(%.2f,%d)", sigmas->data.F32[i], orders->data.S32[i]);
     401        psStringAppend(&params, "(%.2f,%d)", fwhms->data.F32[i], orders->data.S32[i]);
    400402    }
    401403
    … …  
    419421    // Set the Gaussian kernel parameters
    420422    for (int i = 0, index = 0; i < numGaussians; i++) {
    421         float norm = 1.0 / (M_2_PI * sqrtf(sigmas->data.F32[i])); // Normalisation for Gaussian
     423        float sigma = fwhms->data.F32[i] / (2.0 * sqrtf(2.0 * logf(2.0))); // Gaussian sigma
     424        float norm = 1.0 / (M_2_PI * sqrtf(sigma)); // Normalisation for Gaussian
    422425        // Iterate over (u,v) order
    423426        for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
    … …  
    429432                    for (int u = -size; u <= size; u++) {
    430433                        sum += preCalc->kernel[v][u] = norm * power(u, uOrder) * power(v, vOrder) *
    431                             expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigmas->data.F32[i]));
     434                            expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(fwhms->data.F32[i]));
    432435                    }
    433436                }
    … …  
    448451                }
    449452
    450                 kernels->widths->data.F32[index] = sigmas->data.F32[i];
     453                kernels->widths->data.F32[index] = fwhms->data.F32[i];
    451454                kernels->u->data.S32[index] = uOrder;
    452455                kernels->v->data.S32[index] = vOrder;
    … …  
    454457
    455458                psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index,
    456                         sigmas->data.F32[i], uOrder, vOrder);
     459                        fwhms->data.F32[i], uOrder, vOrder);
    457460            }
    458461        }
    … …  
    616619
    617620pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, int size, int spatialOrder,
    618                                                    const psVector *sigmas, const psVector *orders, int inner,
     621                                                   const psVector *fwhms, const psVector *orders, int inner,
    619622                                                   int binning, int ringsOrder)
    620623{
    … …  
    623626        return pmSubtractionKernelsPOIS(size, spatialOrder);
    624627      case PM_SUBTRACTION_KERNEL_ISIS:
    625         return pmSubtractionKernelsISIS(size, spatialOrder, sigmas, orders);
     628        return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders);
    626629      case PM_SUBTRACTION_KERNEL_SPAM:
    627630        return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning);
    … …  
    629632        return pmSubtractionKernelsFRIES(size, spatialOrder, inner);
    630633      case PM_SUBTRACTION_KERNEL_GUNK:
    631         return pmSubtractionKernelsGUNK(size, spatialOrder, sigmas, orders, inner);
     634        return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner);
    632635      case PM_SUBTRACTION_KERNEL_RINGS:
    633636        return pmSubtractionKernelsRINGS(size, ringsOrder, inner, spatialOrder);
Note: See TracChangeset for help on using the changeset viewer.