IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 13, 2011, 12:19:53 PM (15 years ago)
Author:
eugene
Message:

some code reorganzation (move all types into pmSubtractionTypes.h); add stage to check on different mode and order options, choosing best based on chisq of subtraction; careful with the window and the kernel sizes: if the measured kron radius is too large, allow the window to grow as needed; use the 1st radial moment measured on the stamps to set the kernel scaling; remove some cruft from the code; calling program now passes in kernel scaling options to be used when 1st radial moment is measured; new function to update the kernel description string used for kernel I/O; update the kernel description after the spatial order and direction is chosen; need to carry the kernel fwhms and spatial order to allow for update; calculate the psf-match solution errors; psf solution now uses the weights to get sensible chisq values, window is used to calculate the moments; penalty scale is now adjusted to be consistent with sensible reduced chisq; improved visualizations; use range-limiter in SVD of 1e-10; calculate chisq and moments for the solution and use in the evaluation; split out the stamp selection / convolution steps; reject sources after fitting a model to the chisq assuming non-zero systematic error

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r29597 r30622  
    88#include <pslib.h>
    99
     10#include "pmFPA.h"
     11#include "pmSubtractionTypes.h"
    1012#include "pmSubtraction.h"
    1113#include "pmSubtractionKernels.h"
    … …  
    2022{
    2123    psFree(kernels->description);
     24    psFree(kernels->fwhms);
     25    psFree(kernels->orders);
    2226    psFree(kernels->u);
    2327    psFree(kernels->v);
    … …  
    3034    psFree(kernels->solution1);
    3135    psFree(kernels->solution2);
     36    psFree(kernels->solution1err);
     37    psFree(kernels->solution2err);
    3238    psFree(kernels->sampleStamps);
    3339}
    … …  
    419425
    420426    int num = 0;                        // Number of basis functions
    421     psString params = NULL;             // List of parameters
    422427    for (int i = 0; i < numGaussians; i++) {
    423428        int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian
    424         psStringAppend(&params, "(%.1f,%d)", fwhms->data.F32[i], orders->data.S32[i]);
    425429        num += (gaussOrder + 1) * (gaussOrder + 2) / 2;
    426430    }
    427431
    428     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, size,
    429                                                               spatialOrder, penalty, bounds, mode); // Kernels
    430     psStringAppend(&kernels->description, "ISIS(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
    431 
    432     psLogMsg("psModules.imcombine", PS_LOG_INFO, "ISIS kernel: %s,%d --> %d elements",
    433              params, spatialOrder, num);
    434     psFree(params);
     432    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, size, fwhms, orders, spatialOrder, penalty, bounds, mode); // Kernels
     433    pmSubtractionKernelsMakeDescription(kernels);
     434    psLogMsg("psModules.imcombine", PS_LOG_INFO, "kernel: %s --> %d elements", kernels->description, num);
    435435
    436436# if (!CENTRAL_DELTA && !ZERO_KERNEL_ZERO_FLUX)
    … …  
    503503
    504504    int num = 0;                        // Number of basis functions
    505     psString params = NULL;             // List of parameters
    506505    for (int i = 0; i < numGaussians; i++) {
    507506        int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian
    508         psStringAppend(&params, "(%.1f,%d)", fwhms->data.F32[i], orders->data.S32[i]);
    509507        num += (gaussOrder + 1) * (gaussOrder + 2) / 2;
    510508        num += (11 - gaussOrder - 1);   // include all higher order radial terms
    511509    }
    512510
    513     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS_RADIAL, size,
    514                                                               spatialOrder, penalty, bounds, mode); // Kernels
    515     psStringAppend(&kernels->description, "ISIS_RADIAL(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
    516 
    517     psLogMsg("psModules.imcombine", PS_LOG_INFO, "ISIS_RADIAL kernel: %s,%d --> %d elements", params, spatialOrder, num);
    518     psFree(params);
     511    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS_RADIAL, size, fwhms, orders, spatialOrder, penalty, bounds, mode); // Kernels
     512    pmSubtractionKernelsMakeDescription(kernels);
     513    psLogMsg("psModules.imcombine", PS_LOG_INFO, "kernel: %s --> %d elements", kernels->description, num);
    519514
    520515    // Set the kernel parameters
    … …  
    569564
    570565    int num = 0;                        // Number of basis functions
    571     psString params = NULL;             // List of parameters
    572566    for (int i = 0; i < numGaussians; i++) {
    573567        int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian
    574         psStringAppend(&params, "(%.1f,%d)", fwhms->data.F32[i], orders->data.S32[i]);
    575568        num += (gaussOrder + 1) * (gaussOrder + 2) / 2;
    576569    }
    577570
    578     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_HERM, size,
    579                                                               spatialOrder, penalty, bounds, mode); // Kernels
    580     psStringAppend(&kernels->description, "HERM(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
    581 
    582     psLogMsg("psModules.imcombine", PS_LOG_INFO, "HERM kernel: %s,%d --> %d elements",
    583              params, spatialOrder, num);
    584     psFree(params);
     571    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_HERM, size, fwhms, orders, spatialOrder, penalty, bounds, mode); // Kernels
     572    pmSubtractionKernelsMakeDescription(kernels);
     573    psLogMsg("psModules.imcombine", PS_LOG_INFO, "kernel: %s --> %d elements", kernels->description, num);
    585574
    586575    // Set the kernel parameters
    … …  
    627616
    628617    int num = 0;                        // Number of basis functions
    629     psString params = NULL;             // List of parameters
    630618    for (int i = 0; i < numGaussians; i++) {
    631619        int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian
    632         psStringAppend(&params, "(%.1f,%d)", fwhms->data.F32[i], orders->data.S32[i]);
    633620        num += PS_SQR(gaussOrder + 1);
    634621    }
    635622
    636     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_DECONV_HERM, size,
    637                                                               spatialOrder, penalty, bounds, mode); // Kernels
    638     psStringAppend(&kernels->description, "DECONV_HERM(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
    639 
    640     psLogMsg("psModules.imcombine", PS_LOG_INFO, "DECONVOLVED HERM kernel: %s,%d --> %d elements", params, spatialOrder, num);
    641     psFree(params);
     623    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_DECONV_HERM, size, fwhms, orders, spatialOrder, penalty, bounds, mode); // Kernels
     624    pmSubtractionKernelsMakeDescription(kernels);
     625    psLogMsg("psModules.imcombine", PS_LOG_INFO, "kernel: %s --> %d elements", kernels->description, num);
    642626
    643627    // XXXXX hard-wired reference sigma for now of 1.7 pix (== 4.0 pix fwhm == 1.0 arcsec in simtest)
    … …  
    713697
    714698pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, pmSubtractionKernelsType type,
    715                                                 int size, int spatialOrder, float penalty, psRegion bounds,
     699                                                int size, psVector *fwhms, psVector *orders, int spatialOrder, float penalty, psRegion bounds,
    716700                                                pmSubtractionMode mode)
    717701{
    … …  
    722706    kernels->description = NULL;
    723707    kernels->num = numBasisFunctions;
     708    kernels->fwhms = psMemIncrRefCounter(fwhms);
     709    kernels->orders = psMemIncrRefCounter(orders);
    724710    kernels->u = psVectorAlloc(numBasisFunctions, PS_TYPE_S32);
    725711    kernels->v = psVectorAlloc(numBasisFunctions, PS_TYPE_S32);
    … …  
    740726    kernels->size = size;
    741727    kernels->inner = 0;
     728    kernels->binning = 0;
     729    kernels->ringsOrder = 0;
    742730    kernels->spatialOrder = spatialOrder;
    743731    kernels->bgOrder = 0;
    … …  
    745733    kernels->solution1 = NULL;
    746734    kernels->solution2 = NULL;
     735    kernels->solution1err = NULL;
     736    kernels->solution2err = NULL;
    747737    kernels->mean = NAN;
    748738    kernels->rms = NAN;
    … …  
    823813    int num = PS_SQR(2 * size + 1) - 1; // Number of basis functions
    824814
    825     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(0, PM_SUBTRACTION_KERNEL_POIS, size,
    826                                                               spatialOrder, penalty, bounds, mode); // Kernels
    827     psStringAppend(&kernels->description, "POIS(%d,%d,%.2e)", size, spatialOrder, penalty);
    828     psLogMsg("psModules.imcombine", PS_LOG_INFO, "POIS kernel: %d,%d --> %d elements",
    829              size, spatialOrder, num);
     815    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(0, PM_SUBTRACTION_KERNEL_POIS, size, NULL, NULL, spatialOrder, penalty, bounds, mode); // Kernels
     816    pmSubtractionKernelsMakeDescription(kernels);
     817    psLogMsg("psModules.imcombine", PS_LOG_INFO, "kernel: %s --> %d elements", kernels->description, num);
    830818
    831819    if (!p_pmSubtractionKernelsAddGrid(kernels, 0, size)) {
    … …  
    871859    psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num);
    872860
    873     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, size,
    874                                                               spatialOrder, penalty, bounds, mode); // Kernels
     861    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, size, NULL, NULL, spatialOrder, penalty, bounds, mode); // Kernels
    875862    kernels->inner = inner;
    876     psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d,%.2e)", size, inner, binning, spatialOrder,
    877                    penalty);
    878 
    879     psLogMsg("psModules.imcombine", PS_LOG_INFO, "SPAM kernel: %d,%d,%d,%d --> %d elements",
    880              size, inner, binning, spatialOrder, num);
     863    kernels->binning = binning;
     864    pmSubtractionKernelsMakeDescription(kernels);
     865    psLogMsg("psModules.imcombine", PS_LOG_INFO, "kernel: %s --> %d elements", kernels->description, num);
    881866
    882867    kernels->uStop = psVectorAlloc(num, PS_TYPE_S32);
    … …  
    970955    psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num);
    971956
    972     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES, size,
    973                                                               spatialOrder, penalty, bounds, mode); // Kernels
     957    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES, size, NULL, NULL, spatialOrder, penalty, bounds, mode); // Kernels
    974958    kernels->inner = inner;
    975     psStringAppend(&kernels->description, "FRIES(%d,%d,%d,%.2e)", size, inner, spatialOrder, penalty);
    976 
    977     psLogMsg("psModules.imcombine", PS_LOG_INFO, "FRIES kernel: %d,%d,%d --> %d elements",
    978              size, inner, spatialOrder, num);
     959    pmSubtractionKernelsMakeDescription(kernels);
     960    psLogMsg("psModules.imcombine", PS_LOG_INFO, "kernel: %s --> %d elements", kernels->description, num);
    979961
    980962    kernels->uStop = psVectorAlloc(num, PS_TYPE_S32);
    … …  
    10531035    PS_ASSERT_INT_LESS_THAN(inner, size, NULL);
    10541036
    1055     pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders,
    1056                                                                   penalty, bounds, mode); // Kernels
     1037    pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders, penalty, bounds, mode); // Kernels
    10571038    kernels->type = PM_SUBTRACTION_KERNEL_GUNK;
    1058     psStringPrepend(&kernels->description, "GUNK=");
    1059     psStringAppend(&kernels->description, "+POIS(%d,%d)", inner, spatialOrder);
     1039    kernels->inner = inner;
     1040    pmSubtractionKernelsMakeDescription(kernels);
     1041    psLogMsg("psModules.imcombine", PS_LOG_INFO, "kernel: %s --> %d elements", kernels->description, (int) kernels->num);
    10601042
    10611043    int numISIS = kernels->num;         // Number of ISIS kernels
    … …  
    11001082    int num = numRings * numPoly; // Total number of basis functions
    11011083
    1102     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS, size,
    1103                                                               spatialOrder, penalty, bounds, mode); // Kernels
     1084    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS, size, NULL, NULL, spatialOrder, penalty, bounds, mode); // Kernels
    11041085    kernels->inner = inner;
    1105     psStringAppend(&kernels->description, "RINGS(%d,%d,%d,%d,%.2e)", size, inner, ringsOrder, spatialOrder,
    1106                    penalty);
    1107 
    1108     psLogMsg("psModules.imcombine", PS_LOG_INFO, "RINGS kernel: %d,%d,%d,%d --> %d elements",
    1109              size, inner, ringsOrder, spatialOrder, num);
     1086    kernels->ringsOrder = ringsOrder;
     1087    pmSubtractionKernelsMakeDescription(kernels);
     1088    psLogMsg("psModules.imcombine", PS_LOG_INFO, "kernel: %s --> %d elements", kernels->description, num);
    11101089
    11111090    // Set the Gaussian kernel parameters
    … …  
    14051384}
    14061385
     1386bool pmSubtractionKernelsMakeDescription(pmSubtractionKernels *kernels) {
     1387
     1388    // free if it exists
     1389    psFree (kernels->description);
     1390
     1391    // generate the description parameter string
     1392    psString params = NULL;
     1393    if (kernels->fwhms) {
     1394        for (int i = 0; i < kernels->fwhms->n; i++) {
     1395            psStringAppend(&params, "(%.1f,%d)", kernels->fwhms->data.F32[i], kernels->orders->data.S32[i]);
     1396        }
     1397    }
     1398
     1399    switch (kernels->type) {
     1400      case PM_SUBTRACTION_KERNEL_ISIS:
     1401        psStringAppend (&kernels->description, "ISIS(%d,%s,%d,%.2e)", kernels->size, params, kernels->spatialOrder, kernels->penalty);
     1402        break;
     1403
     1404      case PM_SUBTRACTION_KERNEL_ISIS_RADIAL:
     1405        psStringAppend(&kernels->description, "ISIS_RADIAL(%d,%s,%d,%.2e)", kernels->size, params, kernels->spatialOrder, kernels->penalty);
     1406        break;
     1407
     1408      case PM_SUBTRACTION_KERNEL_HERM:
     1409        psStringAppend(&kernels->description, "HERM(%d,%s,%d,%.2e)", kernels->size, params, kernels->spatialOrder, kernels->penalty);
     1410        break;
     1411
     1412      case PM_SUBTRACTION_KERNEL_DECONV_HERM:
     1413        psStringAppend(&kernels->description, "DECONV_HERM(%d,%s,%d,%.2e)", kernels->size, params, kernels->spatialOrder, kernels->penalty);
     1414        break;
     1415
     1416      case PM_SUBTRACTION_KERNEL_POIS:
     1417        psStringAppend(&kernels->description, "POIS(%d,%d,%.2e)", kernels->size, kernels->spatialOrder, kernels->penalty);
     1418        break;
     1419
     1420      case PM_SUBTRACTION_KERNEL_SPAM:
     1421        psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d,%.2e)", kernels->size, kernels->inner, kernels->binning, kernels->spatialOrder, kernels->penalty);
     1422        break;
     1423
     1424      case PM_SUBTRACTION_KERNEL_FRIES:
     1425        psStringAppend(&kernels->description, "FRIES(%d,%d,%d,%.2e)", kernels->size, kernels->inner, kernels->spatialOrder, kernels->penalty);
     1426        break;
     1427
     1428        // Grid United with Normal Kernel [description: GUNK=ISIS(...)+POIS(...)]
     1429      case PM_SUBTRACTION_KERNEL_GUNK:
     1430        psStringAppend(&kernels->description, "GUNK=ISIS(%d,%s,%d,%.2e)", kernels->size, params, kernels->spatialOrder, kernels->penalty);
     1431        psStringAppend(&kernels->description, "+POIS(%d,%d)", kernels->inner, kernels->spatialOrder);
     1432        break;
     1433       
     1434      case PM_SUBTRACTION_KERNEL_RINGS:
     1435        psStringAppend(&kernels->description, "RINGS(%d,%d,%d,%d,%.2e)", kernels->size, kernels->inner, kernels->ringsOrder, kernels->spatialOrder, kernels->penalty);
     1436        break;
     1437
     1438      default:
     1439        psAbort("unknown kernel");
     1440    }
     1441    psFree (params);
     1442    return true;
     1443}
     1444
    14071445pmSubtractionKernels *pmSubtractionKernelsCopy(const pmSubtractionKernels *in)
    14081446{
    … …  
    14351473    out->solution1 = in->solution1 ? psVectorCopy(NULL, in->solution1, PS_TYPE_F64) : NULL;
    14361474    out->solution2 = in->solution2 ? psVectorCopy(NULL, in->solution2, PS_TYPE_F64) : NULL;
     1475    out->solution1err = in->solution1err ? psVectorCopy(NULL, in->solution1err, PS_TYPE_F64) : NULL;
     1476    out->solution2err = in->solution2err ? psVectorCopy(NULL, in->solution2err, PS_TYPE_F64) : NULL;
    14371477    out->sampleStamps = psMemIncrRefCounter(in->sampleStamps);
    14381478
Note: See TracChangeset for help on using the changeset viewer.