IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 7, 2008, 12:47:04 PM (18 years ago)
Author:
Paul Price
Message:

Statistics are only required when fitting a trend. If we load a trend (e.g., pmPSFmodelRead) then we don't need the statistics. Similar to psImageMap.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/objects/pmTrend2D.c

    r16065 r19961  
    33 *  @author EAM, IfA
    44 *
    5  *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
    6  *  @date $Date: 2008-01-15 02:47:51 $
     5 *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
     6 *  @date $Date: 2008-10-07 22:47:04 $
    77 *  Copyright 2004 Institute for Astronomy, University of Hawaii
    88 *
    … …  
    1313#endif
    1414
    15 # include <strings.h>
    16 # include <pslib.h>
    17 # include "pmTrend2D.h"
    18 
    19 static void pmTrend2DFree (pmTrend2D *trend) {
    20 
    21     if (trend == NULL)
    22         return;
    23 
    24     psFree (trend->stats);
    25     psFree (trend->poly);
    26     psFree (trend->map);
     15#include <strings.h>
     16#include <pslib.h>
     17#include "pmTrend2D.h"
     18
     19static void pmTrend2DFree(pmTrend2D *trend)
     20{
     21    psFree(trend->stats);
     22    psFree(trend->poly);
     23    psFree(trend->map);
    2724    return;
    2825}
    2926
    30 pmTrend2D *pmTrend2DAlloc (pmTrend2DMode mode, psImage *image, int nXtrend, int nYtrend, psStats *stats)
    31 {
    32     PS_ASSERT_PTR_NON_NULL(stats, NULL);
     27pmTrend2D *pmTrend2DAlloc(pmTrend2DMode mode, psImage *image, int nXtrend, int nYtrend, psStats *stats)
     28{
    3329    if (mode == PM_TREND_MAP) {
    34         PS_ASSERT_PTR_NON_NULL(image, NULL);
    35     }
    36 
    37     pmTrend2D *trend = (pmTrend2D *) psAlloc(sizeof(pmTrend2D));
    38     psMemSetDeallocator(trend, (psFreeFunc) pmTrend2DFree);
     30        psAssert(image, "Need an image for MAP trend mode");
     31    }
     32
     33    pmTrend2D *trend = psAlloc(sizeof(pmTrend2D));
     34    psMemSetDeallocator(trend, (psFreeFunc)pmTrend2DFree);
    3935
    4036    trend->map = NULL;
    4137    trend->poly = NULL;
    42     trend->stats = psMemIncrRefCounter (stats);
     38    trend->stats = psMemIncrRefCounter(stats);
    4339    trend->mode = mode;
    44        
    45     switch (mode) {
    46       case PM_TREND_POLY_ORD:
    47         trend->poly = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXtrend, nYtrend);
    48         // set masking somehow
    49         for (int nx = 0; nx < trend->poly->nX + 1; nx++) {
    50             for (int ny = 0; ny < trend->poly->nY + 1; ny++) {
    51                 if (nx + ny >= PS_MAX (trend->poly->nX, trend->poly->nY) + 1) {
    52                     trend->poly->coeffMask[nx][ny] = PS_POLY_MASK_SET;
    53                 } else {
    54                     trend->poly->coeffMask[nx][ny] = PS_POLY_MASK_NONE;
    55                 }
    56             }
    57         }
    58         break;
    59 
    60       case PM_TREND_POLY_CHEB:
    61         trend->poly = psPolynomial2DAlloc (PS_POLYNOMIAL_CHEB, nXtrend, nYtrend);
    62         break;
     40
     41    switch (mode) {
     42      case PM_TREND_POLY_ORD:
     43        trend->poly = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, nXtrend, nYtrend);
     44        // set masking somehow
     45        for (int nx = 0; nx < trend->poly->nX + 1; nx++) {
     46            for (int ny = 0; ny < trend->poly->nY + 1; ny++) {
     47                if (nx + ny >= PS_MAX (trend->poly->nX, trend->poly->nY) + 1) {
     48                    trend->poly->coeffMask[nx][ny] = PS_POLY_MASK_SET;
     49                } else {
     50                    trend->poly->coeffMask[nx][ny] = PS_POLY_MASK_NONE;
     51                }
     52            }
     53        }
     54        break;
     55
     56      case PM_TREND_POLY_CHEB:
     57        trend->poly = psPolynomial2DAlloc(PS_POLYNOMIAL_CHEB, nXtrend, nYtrend);
     58        break;
    6359
    6460      case PM_TREND_MAP: {
    65           // binning defines the map scale relationship
    66           psImageBinning *binning = psImageBinningAlloc();
    67           binning->nXruff = nXtrend;
    68           binning->nYruff = nYtrend;
    69           binning->nXfine = image->numCols;
    70           binning->nYfine = image->numRows;
    71 
    72           trend->map = psImageMapAlloc (image, binning, stats);
    73           psFree (binning);
    74           break;
     61          // binning defines the map scale relationship
     62          psImageBinning *binning = psImageBinningAlloc();
     63          binning->nXruff = nXtrend;
     64          binning->nYruff = nYtrend;
     65          binning->nXfine = image->numCols;
     66          binning->nYfine = image->numRows;
     67
     68          trend->map = psImageMapAlloc(image, binning, stats);
     69          psFree(binning);
     70          break;
    7571      }
    7672      // XXX: Put a more graceful error here.
    7773      default:
    78         psAbort ("error");
    79     }
    80     return (trend);
     74        psAbort("error");
     75    }
     76    return trend;
    8177}
    8278
    … …  
    8783}
    8884
    89 pmTrend2D *pmTrend2DNoImageAlloc (pmTrend2DMode mode, psImageBinning *binning, psStats *stats)
     85pmTrend2D *pmTrend2DNoImageAlloc(pmTrend2DMode mode, psImageBinning *binning, psStats *stats)
    9086{
    9187    if (mode == PM_TREND_MAP) {
    92         PS_ASSERT_PTR_NON_NULL(binning, NULL);
    93         PS_ASSERT_PTR_NON_NULL(stats, NULL);
    94     }
    95     pmTrend2D *trend = (pmTrend2D *) psAlloc(sizeof(pmTrend2D));
    96     psMemSetDeallocator(trend, (psFreeFunc) pmTrend2DFree);
     88        psAssert(binning, "Need binning for MAP mode");
     89    }
     90    pmTrend2D *trend = psAlloc(sizeof(pmTrend2D));
     91    psMemSetDeallocator(trend, (psFreeFunc)pmTrend2DFree);
    9792
    9893    trend->map = NULL;
    9994    trend->poly = NULL;
    100     trend->stats = psMemIncrRefCounter (stats);
     95    trend->stats = psMemIncrRefCounter(stats);
    10196    trend->mode = mode;
    102        
    103     switch (mode) {
    104       case PM_TREND_POLY_ORD:
    105         trend->poly = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, binning->nXruff, binning->nYruff);
    106         // set masking somehow
    107         for (int nx = 0; nx < trend->poly->nX + 1; nx++) {
    108             for (int ny = 0; ny < trend->poly->nY + 1; ny++) {
    109                 if (nx + ny >= PS_MAX (trend->poly->nX, trend->poly->nY) + 1) {
    110                     trend->poly->coeffMask[nx][ny] = PS_POLY_MASK_SET;
    111                 } else {
    112                     trend->poly->coeffMask[nx][ny] = PS_POLY_MASK_NONE;
    113                 }
    114             }
    115         }
    116         break;
    117 
    118       case PM_TREND_POLY_CHEB:
    119         trend->poly = psPolynomial2DAlloc (PS_POLYNOMIAL_CHEB, binning->nXruff, binning->nYruff);
    120         break;
     97
     98    switch (mode) {
     99      case PM_TREND_POLY_ORD:
     100        trend->poly = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, binning->nXruff, binning->nYruff);
     101        // set masking somehow
     102        for (int nx = 0; nx < trend->poly->nX + 1; nx++) {
     103            for (int ny = 0; ny < trend->poly->nY + 1; ny++) {
     104                if (nx + ny >= PS_MAX (trend->poly->nX, trend->poly->nY) + 1) {
     105                    trend->poly->coeffMask[nx][ny] = PS_POLY_MASK_SET;
     106                } else {
     107                    trend->poly->coeffMask[nx][ny] = PS_POLY_MASK_NONE;
     108                }
     109            }
     110        }
     111        break;
     112
     113      case PM_TREND_POLY_CHEB:
     114        trend->poly = psPolynomial2DAlloc(PS_POLYNOMIAL_CHEB, binning->nXruff, binning->nYruff);
     115        break;
    121116
    122117      case PM_TREND_MAP: {
    123           // binning defines the map scale relationship
    124           trend->map = psImageMapNoImageAlloc (binning, stats);
    125           break;
     118          // binning defines the map scale relationship
     119          trend->map = psImageMapNoImageAlloc(binning, stats);
     120          break;
    126121      }
    127122
    128123      default:
    129         psAbort ("error");
    130     }
    131     return (trend);
    132 }
    133 
    134 pmTrend2D *pmTrend2DFieldAlloc (pmTrend2DMode mode, int nXfield, int nYfield, int nXtrend, int nYtrend, psStats *stats)
    135 {
    136     PS_ASSERT_PTR_NON_NULL(stats, NULL);
    137     pmTrend2D *trend = (pmTrend2D *) psAlloc(sizeof(pmTrend2D));
    138     psMemSetDeallocator(trend, (psFreeFunc) pmTrend2DFree);
     124        psAbort("error");
     125    }
     126    return trend;
     127}
     128
     129pmTrend2D *pmTrend2DFieldAlloc(pmTrend2DMode mode, int nXfield, int nYfield,
     130                               int nXtrend, int nYtrend, psStats *stats)
     131{
     132    psAssert(stats, "Require statistics");
     133
     134    pmTrend2D *trend = psAlloc(sizeof(pmTrend2D));
     135    psMemSetDeallocator(trend, (psFreeFunc)pmTrend2DFree);
    139136
    140137    trend->map = NULL;
    141138    trend->poly = NULL;
    142     trend->stats = psMemIncrRefCounter (stats);
     139    trend->stats = psMemIncrRefCounter(stats);
    143140    trend->mode = mode;
    144        
    145     switch (mode) {
    146       case PM_TREND_POLY_ORD:
    147         trend->poly = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXtrend, nYtrend);
    148         // set masking somehow
    149         for (int nx = 0; nx < trend->poly->nX + 1; nx++) {
    150             for (int ny = 0; ny < trend->poly->nY + 1; ny++) {
    151                 if (nx + ny >= PS_MAX (trend->poly->nX, trend->poly->nY) + 1) {
    152                     trend->poly->coeffMask[nx][ny] = PS_POLY_MASK_SET;
    153                 } else {
    154                     trend->poly->coeffMask[nx][ny] = PS_POLY_MASK_NONE;
    155                 }
    156             }
    157         }
    158         break;
    159 
    160       case PM_TREND_POLY_CHEB:
    161         trend->poly = psPolynomial2DAlloc (PS_POLYNOMIAL_CHEB, nXtrend, nYtrend);
    162         break;
     141
     142    switch (mode) {
     143      case PM_TREND_POLY_ORD:
     144        trend->poly = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, nXtrend, nYtrend);
     145        // set masking somehow
     146        for (int nx = 0; nx < trend->poly->nX + 1; nx++) {
     147            for (int ny = 0; ny < trend->poly->nY + 1; ny++) {
     148                if (nx + ny >= PS_MAX (trend->poly->nX, trend->poly->nY) + 1) {
     149                    trend->poly->coeffMask[nx][ny] = PS_POLY_MASK_SET;
     150                } else {
     151                    trend->poly->coeffMask[nx][ny] = PS_POLY_MASK_NONE;
     152                }
     153            }
     154        }
     155        break;
     156
     157      case PM_TREND_POLY_CHEB:
     158        trend->poly = psPolynomial2DAlloc(PS_POLYNOMIAL_CHEB, nXtrend, nYtrend);
     159        break;
    163160
    164161      case PM_TREND_MAP: {
    165           // binning defines the map scale relationship
    166           psImageBinning *binning = psImageBinningAlloc();
    167           binning->nXfine = nXfield;
    168           binning->nYfine = nYfield;
    169           binning->nXruff = nXtrend;
    170           binning->nYruff = nYtrend;
    171 
    172           trend->map = psImageMapAlloc (NULL, binning, stats);
    173           psFree (binning);
    174           break;
     162          // binning defines the map scale relationship
     163          psImageBinning *binning = psImageBinningAlloc();
     164          binning->nXfine = nXfield;
     165          binning->nYfine = nYfield;
     166          binning->nXruff = nXtrend;
     167          binning->nYruff = nYtrend;
     168
     169          trend->map = psImageMapAlloc(NULL, binning, stats);
     170          psFree (binning);
     171          break;
    175172      }
    176173
    177174      default:
    178175        // XXX: Put a more graceful error here.
    179         psAbort ("error");
    180     }
    181     return (trend);
    182 }
    183 
    184 bool pmTrend2DFit (pmTrend2D *trend, psVector *mask, psMaskType maskVal, psVector *x,
    185                    psVector *y, psVector *f, psVector *df)
    186 {
    187     PS_ASSERT_PTR_NON_NULL(trend, false);
     176        psAbort("error");
     177    }
     178    return trend;
     179}
     180
     181bool pmTrend2DFit(pmTrend2D *trend, psVector *mask, psMaskType maskVal, const psVector *x,
     182                  const psVector *y, const psVector *f, const psVector *df)
     183{
     184    PM_ASSERT_TREND2D_NON_NULL(trend, false);
     185    PM_ASSERT_TREND2D_STATS(trend, false);
    188186    PS_ASSERT_VECTOR_NON_NULL(x, false);
    189187    PS_ASSERT_VECTOR_NON_NULL(y, false);
    … …  
    194192      case PM_TREND_POLY_ORD:
    195193      case PM_TREND_POLY_CHEB:
    196         status = psVectorClipFitPolynomial2D (trend->poly, trend->stats, mask, maskVal, f, df, x, y);
    197         // we can use the API here which adjusts the polynomial order based on the number
    198         // of points in the image, and potentially based on the fractional range of the
    199         // data?
    200         break;
    201 
    202       case PM_TREND_MAP:
    203         // XXX supply fraction from trend elements
    204         // XXX need to add the API which adjusts the scale
    205         status = psImageMapClipFit (trend->map, trend->stats, mask, maskVal, x, y, f, df);
    206         break;
    207 
    208       default:
    209         psAbort ("error");
     194        status = psVectorClipFitPolynomial2D(trend->poly, trend->stats, mask, maskVal, f, df, x, y);
     195        // we can use the API here which adjusts the polynomial order based on the number
     196        // of points in the image, and potentially based on the fractional range of the
     197        // data?
     198        break;
     199
     200      case PM_TREND_MAP:
     201        // XXX supply fraction from trend elements
     202        // XXX need to add the API which adjusts the scale
     203        status = psImageMapClipFit(trend->map, trend->stats, mask, maskVal, x, y, f, df);
     204        break;
     205
     206      default:
     207        psAbort ("error");
    210208    }
    211209    return status;
    212210}
    213211
    214 double pmTrend2DEval (pmTrend2D *trend, float x, float y)
    215 {
    216     if (!trend) return 0.0;
     212double pmTrend2DEval(const pmTrend2D *trend, float x, float y)
     213{
     214    // This might be in a tight loop, so no complicated assertions
     215    if (!trend) {
     216        return 0.0;
     217    }
    217218
    218219    double result;
    … …  
    220221      case PM_TREND_POLY_ORD:
    221222      case PM_TREND_POLY_CHEB:
    222         result = psPolynomial2DEval (trend->poly, x, y);
    223         break;
    224 
    225       case PM_TREND_MAP:
    226         result = psImageMapEval (trend->map, x, y);
    227         break;
    228 
    229       default:
    230         psAbort ("error");
     223        result = psPolynomial2DEval(trend->poly, x, y);
     224        break;
     225
     226      case PM_TREND_MAP:
     227        result = psImageMapEval(trend->map, x, y);
     228        break;
     229
     230      default:
     231        psAbort ("error");
    231232    }
    232233    return result;
    233234}
    234235
    235 psVector *pmTrend2DEvalVector (pmTrend2D *trend, psVector *x, psVector *y)
    236 {
    237     PS_ASSERT_PTR_NON_NULL(trend, NULL);
     236psVector *pmTrend2DEvalVector(const pmTrend2D *trend, const psVector *x, const psVector *y)
     237{
     238    PM_ASSERT_TREND2D_NON_NULL(trend, NULL);
    238239    PS_ASSERT_VECTOR_NON_NULL(x, false);
    239240    PS_ASSERT_VECTOR_NON_NULL(y, false);
    … …  
    243244      case PM_TREND_POLY_ORD:
    244245      case PM_TREND_POLY_CHEB:
    245         result = psPolynomial2DEvalVector (trend->poly, x, y);
    246         break;
    247 
    248       case PM_TREND_MAP:
    249         result = psImageMapEvalVector (trend->map, x, y);
    250         break;
    251 
    252       default:
    253         psAbort ("error");
     246        result = psPolynomial2DEvalVector (trend->poly, x, y);
     247        break;
     248
     249      case PM_TREND_MAP:
     250        result = psImageMapEvalVector (trend->map, x, y);
     251        break;
     252
     253      default:
     254        psAbort ("error");
    254255    }
    255256    return result;
    256257}
    257258
    258 psString pmTrend2DModeToString (pmTrend2DMode mode) {
    259    
    260     psString name;
    261 
     259psString pmTrend2DModeToString(pmTrend2DMode mode)
     260{
    262261    switch (mode) {
    263262      case PM_TREND_NONE:
    264         name = psStringCopy ("NONE");
    265         break;
    266       case PM_TREND_POLY_ORD:
    267         name = psStringCopy ("POLY_ORD");
    268         break;
    269       case PM_TREND_POLY_CHEB:
    270         name = psStringCopy ("POLY_CHEB");
    271         break;
    272       case PM_TREND_MAP:
    273         name = psStringCopy ("MAP");
    274         break;
    275       default:
    276         psError (PS_ERR_UNKNOWN, true, "Unknown pmTrend2D mode\n");
    277         psAbort ("invalid mode %d", mode);
    278     }
    279     return name;
    280 }
    281 
    282 pmTrend2DMode pmTrend2DModeFromString (psString name) {
    283 
    284     if (!name) return PM_TREND_NONE;
    285 
    286     if (!strcasecmp (name, "NONE")) {
    287         return PM_TREND_NONE;
    288     }
    289     if (!strcasecmp (name, "POLY_ORD")) {
    290         return PM_TREND_POLY_ORD;
    291     }
    292     if (!strcasecmp (name, "POLY_CHEB")) {
    293         return PM_TREND_POLY_CHEB;
    294     }
    295     if (!strcasecmp (name, "MAP")) {
    296         return PM_TREND_MAP;
    297     }
    298     psError (PS_ERR_UNKNOWN, true, "Unknown pmTrend2D mode %s\n", name);
     263        return psStringCopy("NONE");
     264      case PM_TREND_POLY_ORD:
     265        return psStringCopy("POLY_ORD");
     266        break;
     267      case PM_TREND_POLY_CHEB:
     268        return psStringCopy("POLY_CHEB");
     269      case PM_TREND_MAP:
     270        return psStringCopy("MAP");
     271        break;
     272      default:
     273        psError(PS_ERR_UNKNOWN, true, "Unknown pmTrend2D mode");
     274    }
     275    psAbort("invalid mode %d", mode);
     276}
     277
     278pmTrend2DMode pmTrend2DModeFromString(psString name)
     279 {
     280    if (!name) {
     281        return PM_TREND_NONE;
     282    }
     283
     284    if (!strcasecmp(name, "NONE")) {
     285        return PM_TREND_NONE;
     286    }
     287    if (!strcasecmp(name, "POLY_ORD")) {
     288        return PM_TREND_POLY_ORD;
     289    }
     290    if (!strcasecmp(name, "POLY_CHEB")) {
     291        return PM_TREND_POLY_CHEB;
     292    }
     293    if (!strcasecmp(name, "MAP")) {
     294        return PM_TREND_MAP;
     295    }
     296    psError(PS_ERR_UNKNOWN, true, "Unknown pmTrend2D mode %s", name);
    299297    return PM_TREND_NONE;
    300298}
Note: See TracChangeset for help on using the changeset viewer.