IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jan 29, 2010, 6:02:38 PM (16 years ago)
Author:
Paul Price
Message:

Fairly extensive change to make subtractions more stable on images with limited lit area in the presence of spatial variation. Took the polynomials which were scaled over the entire image and scale them only over the lit area. Removed some fprintf statements. Fixed the problem causing the background term to always come out 1.0.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r26734 r26739  
    197197    }
    198198
    199 #if 1
     199#if 0
    200200    fprintf(stderr, "%d raw: %lf, null: %f, min: %lf, max: %lf, moment: %lf\n", index, sum, preCalc->kernel->kernel[0][0], min, max, penalty);
    201201#endif
    … …  
    245245    }
    246246
    247 #if 1
     247#if 0
    248248    {
    249249        double sum = 0.0;   // Sum of kernel component
    … …  
    276276pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, int spatialOrder,
    277277                                                    const psVector *fwhmsIN, const psVector *ordersIN,
    278                                                     float penalty, pmSubtractionMode mode)
     278                                                    float penalty, psRegion bounds, pmSubtractionMode mode)
    279279{
    280280    PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL);
    … …  
    305305    }
    306306
    307     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, size, spatialOrder, penalty, mode); // The kernels
     307    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, size,
     308                                                              spatialOrder, penalty, bounds, mode); // Kernels
    308309    psStringAppend(&kernels->description, "ISIS(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
    309310
    … …  
    331332pmSubtractionKernels *pmSubtractionKernelsISIS_RADIAL(int size, int spatialOrder,
    332333                                                      const psVector *fwhmsIN, const psVector *ordersIN,
    333                                                       float penalty, pmSubtractionMode mode)
     334                                                      float penalty, psRegion bounds, pmSubtractionMode mode)
    334335{
    335336    PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL);
    … …  
    361362    }
    362363
    363     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS_RADIAL, size, spatialOrder, penalty, mode); // The kernels
     364    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS_RADIAL, size,
     365                                                              spatialOrder, penalty, bounds, mode); // Kernels
    364366    psStringAppend(&kernels->description, "ISIS_RADIAL(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
    365367
    … …  
    388390pmSubtractionKernels *pmSubtractionKernelsHERM(int size, int spatialOrder,
    389391                                               const psVector *fwhmsIN, const psVector *ordersIN,
    390                                                float penalty, pmSubtractionMode mode)
     392                                               float penalty, psRegion bounds, pmSubtractionMode mode)
    391393{
    392394    PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL);
    … …  
    417419    }
    418420
    419     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_HERM, size, spatialOrder, penalty, mode); // The kernels
     421    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_HERM, size,
     422                                                              spatialOrder, penalty, bounds, mode); // Kernels
    420423    psStringAppend(&kernels->description, "HERM(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
    421424
    422     psLogMsg("psModules.imcombine", PS_LOG_INFO, "HERM kernel: %s,%d --> %d elements", params, spatialOrder, num);
     425    psLogMsg("psModules.imcombine", PS_LOG_INFO, "HERM kernel: %s,%d --> %d elements",
     426             params, spatialOrder, num);
    423427    psFree(params);
    424428
    … …  
    439443
    440444pmSubtractionKernels *pmSubtractionKernelsDECONV_HERM(int size, int spatialOrder,
    441                                                      const psVector *fwhmsIN, const psVector *ordersIN,
    442                                                      float penalty, pmSubtractionMode mode)
     445                                                      const psVector *fwhmsIN, const psVector *ordersIN,
     446                                                      float penalty, psRegion bounds, pmSubtractionMode mode)
    443447{
    444448    PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL);
    … …  
    469473    }
    470474
    471     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_DECONV_HERM, size, spatialOrder, penalty, mode); // The kernels
     475    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_DECONV_HERM, size,
     476                                                              spatialOrder, penalty, bounds, mode); // Kernels
    472477    psStringAppend(&kernels->description, "DECONV_HERM(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
    473478
    … …  
    544549
    545550pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, pmSubtractionKernelsType type,
    546                                                 int size, int spatialOrder, float penalty,
     551                                                int size, int spatialOrder, float penalty, psRegion bounds,
    547552                                                pmSubtractionMode mode)
    548553{
    … …  
    558563    kernels->uStop = NULL;
    559564    kernels->vStop = NULL;
     565    kernels->xMin = bounds.x0;
     566    kernels->xMax = bounds.x1;
     567    kernels->yMin = bounds.y0;
     568    kernels->yMax = bounds.y1;
    560569    kernels->preCalc = psArrayAlloc(numBasisFunctions);
    561570    kernels->penalty = penalty;
    … …  
    566575    kernels->bgOrder = 0;
    567576    kernels->mode = mode;
    568     kernels->numCols = 0;
    569     kernels->numRows = 0;
    570577    kernels->solution1 = NULL;
    571578    kernels->solution2 = NULL;
    … …  
    640647}
    641648
    642 pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, float penalty,
     649pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, float penalty, psRegion bounds,
    643650                                               pmSubtractionMode mode)
    644651{
    … …  
    649656
    650657    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(0, PM_SUBTRACTION_KERNEL_POIS, size,
    651                                                               spatialOrder, penalty, mode); // The kernels
     658                                                              spatialOrder, penalty, bounds, mode); // Kernels
    652659    psStringAppend(&kernels->description, "POIS(%d,%d,%.2e)", size, spatialOrder, penalty);
    653660    psLogMsg("psModules.imcombine", PS_LOG_INFO, "POIS kernel: %d,%d --> %d elements",
    … …  
    664671pmSubtractionKernels *pmSubtractionKernelsISIS(int size, int spatialOrder,
    665672                                               const psVector *fwhms, const psVector *orders,
    666                                                float penalty, pmSubtractionMode mode)
     673                                               float penalty, psRegion bounds, pmSubtractionMode mode)
    667674{
    668675    pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders,
    669                                                                   penalty, mode); // Kernels
     676                                                                  penalty, bounds, mode); // Kernels
    670677    if (!kernels) {
    671678        return NULL;
    … …  
    676683
    677684pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, int spatialOrder, int inner, int binning,
    678                                                float penalty, pmSubtractionMode mode)
     685                                               float penalty, psRegion bounds, pmSubtractionMode mode)
    679686{
    680687    PS_ASSERT_INT_POSITIVE(size, NULL);
    … …  
    697704
    698705    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, size,
    699                                                               spatialOrder, penalty, mode); // The kernels
     706                                                              spatialOrder, penalty, bounds, mode); // Kernels
    700707    kernels->inner = inner;
    701708    psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d,%.2e)", size, inner, binning, spatialOrder,
    … …  
    768775
    769776pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, int spatialOrder, int inner, float penalty,
    770                                                 pmSubtractionMode mode)
     777                                                psRegion bounds, pmSubtractionMode mode)
    771778{
    772779    PS_ASSERT_INT_POSITIVE(size, NULL);
    … …  
    795802
    796803    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES, size,
    797                                                               spatialOrder, penalty, mode); // The kernels
     804                                                              spatialOrder, penalty, bounds, mode); // Kernels
    798805    kernels->inner = inner;
    799806    psStringAppend(&kernels->description, "FRIES(%d,%d,%d,%.2e)", size, inner, spatialOrder, penalty);
    … …  
    864871pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, int spatialOrder, const psVector *fwhms,
    865872                                               const psVector *orders, int inner, float penalty,
    866                                                pmSubtractionMode mode)
     873                                               psRegion bounds, pmSubtractionMode mode)
    867874{
    868875    PS_ASSERT_INT_POSITIVE(size, NULL);
    … …  
    877884
    878885    pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders,
    879                                                                   penalty, mode); // Kernels
     886                                                                  penalty, bounds, mode); // Kernels
    880887    kernels->type = PM_SUBTRACTION_KERNEL_GUNK;
    881888    psStringPrepend(&kernels->description, "GUNK=");
    … …  
    893900// RINGS --- just what it says
    894901pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, int spatialOrder, int inner, int ringsOrder,
    895                                                 float penalty, pmSubtractionMode mode)
     902                                                float penalty, psRegion bounds, pmSubtractionMode mode)
    896903{
    897904    PS_ASSERT_INT_POSITIVE(size, NULL);
    … …  
    924931
    925932    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS, size,
    926                                                               spatialOrder, penalty, mode); // The kernels
     933                                                              spatialOrder, penalty, bounds, mode); // Kernels
    927934    kernels->inner = inner;
    928935    psStringAppend(&kernels->description, "RINGS(%d,%d,%d,%d,%.2e)", size, inner, ringsOrder, spatialOrder,
    … …  
    10461053pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, int size, int spatialOrder,
    10471054                                                   const psVector *fwhms, const psVector *orders, int inner,
    1048                                                    int binning, int ringsOrder, float penalty,
     1055                                                   int binning, int ringsOrder, float penalty, psRegion bounds,
    10491056                                                   pmSubtractionMode mode)
    10501057{
    10511058    switch (type) {
    10521059      case PM_SUBTRACTION_KERNEL_POIS:
    1053         return pmSubtractionKernelsPOIS(size, spatialOrder, penalty, mode);
     1060        return pmSubtractionKernelsPOIS(size, spatialOrder, penalty, bounds, mode);
    10541061      case PM_SUBTRACTION_KERNEL_ISIS:
    1055         return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, penalty, mode);
     1062        return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, penalty, bounds, mode);
    10561063      case PM_SUBTRACTION_KERNEL_ISIS_RADIAL:
    1057         return pmSubtractionKernelsISIS_RADIAL(size, spatialOrder, fwhms, orders, penalty, mode);
     1064        return pmSubtractionKernelsISIS_RADIAL(size, spatialOrder, fwhms, orders, penalty, bounds, mode);
    10581065      case PM_SUBTRACTION_KERNEL_HERM:
    1059         return pmSubtractionKernelsHERM(size, spatialOrder, fwhms, orders, penalty, mode);
     1066        return pmSubtractionKernelsHERM(size, spatialOrder, fwhms, orders, penalty, bounds, mode);
    10601067      case PM_SUBTRACTION_KERNEL_DECONV_HERM:
    1061         return pmSubtractionKernelsDECONV_HERM(size, spatialOrder, fwhms, orders, penalty, mode);
     1068        return pmSubtractionKernelsDECONV_HERM(size, spatialOrder, fwhms, orders, penalty, bounds, mode);
    10621069      case PM_SUBTRACTION_KERNEL_SPAM:
    1063         return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, penalty, mode);
     1070        return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, penalty, bounds, mode);
    10641071      case PM_SUBTRACTION_KERNEL_FRIES:
    1065         return pmSubtractionKernelsFRIES(size, spatialOrder, inner, penalty, mode);
     1072        return pmSubtractionKernelsFRIES(size, spatialOrder, inner, penalty, bounds, mode);
    10661073      case PM_SUBTRACTION_KERNEL_GUNK:
    1067         return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner, penalty, mode);
     1074        return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner, penalty, bounds, mode);
    10681075      case PM_SUBTRACTION_KERNEL_RINGS:
    1069         return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, penalty, mode);
     1076        return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, penalty, bounds, mode);
    10701077      default:
    10711078        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown kernel type: %x", type);
    … …  
    11031110
    11041111pmSubtractionKernels *pmSubtractionKernelsFromDescription(const char *description, int bgOrder,
    1105                                                           pmSubtractionMode mode)
     1112                                                          psRegion bounds, pmSubtractionMode mode)
    11061113{
    11071114    PS_ASSERT_STRING_NON_EMPTY(description, NULL);
    … …  
    11561163        PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt);
    11571164        penalty = parseStringFloat(ptr);
    1158 
    1159         return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode);
    1160 
     1165        break;
    11611166      case PM_SUBTRACTION_KERNEL_RINGS:
    11621167        PARSE_STRING_NUMBER(size, ptr, ',', parseStringInt);
    … …  
    11651170        PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt);
    11661171        PARSE_STRING_NUMBER(penalty, ptr, ')', parseStringInt);
    1167         return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode);
     1172        break;
    11681173      default:
    11691174        psAbort("Deciphering kernels other than ISIS, HERM, DECONV_HERM or RINGS is not currently supported.");
    11701175    }
    1171     return NULL;
     1176
     1177    return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning,
     1178                                        ringsOrder, penalty, bounds, mode);
    11721179}
    11731180
    … …  
    12451252    out->bgOrder = in->bgOrder;
    12461253    out->mode = in->mode;
    1247     out->numCols = in->numCols;
    1248     out->numRows = in->numRows;
     1254    out->xMin = in->xMin;
     1255    out->xMax = in->xMax;
     1256    out->yMin = in->yMin;
     1257    out->yMax = in->yMax;
    12491258    out->solution1 = in->solution1 ? psVectorCopy(NULL, in->solution1, PS_TYPE_F64) : NULL;
    12501259    out->solution2 = in->solution2 ? psVectorCopy(NULL, in->solution2, PS_TYPE_F64) : NULL;
Note: See TracChangeset for help on using the changeset viewer.