IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 7, 2005, 6:32:40 PM (21 years ago)
Author:
eugene
Message:

cleanup, organization

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psPolynomials.c

    r4954 r4977  
    1 # include "psphot.h"
     1# include <pslib.h>
     2# include "psLibUtils.h"
    23
    34// XXX EAM : this file defines alternate versions of the polynomial fitting
    … …  
    910
    1011// write out the terms of the given 1D polynomial
    11 void psPolynomial1DDump (psPolynomial1D *poly) {
     12void Polynomial1DDump_EAM (psPolynomial1D *poly) {
    1213
    1314    for (int i = 0; i < poly->n + 1; i++) {
    … …  
    1617}   
    1718
    18 psF64 Polynomial1DEval_EAM(psF64 x, const psPolynomial1D* myPoly)
    19 {
    20     psS32 loop_x = 0;
    21     psF64 polySum = 0.0;
    22     psF64 xSum = 1.0;
    23 
    24     for (loop_x = 0; loop_x < myPoly->n + 1; loop_x++) {
    25         if (myPoly->mask[loop_x] == 0) {
    26             polySum += xSum * myPoly->coeff[loop_x];
    27         }
    28         xSum *= x;
    29     }
    30 
    31     return(polySum);
    32 }
    33 
    34 psVector *Polynomial1DEvalVector_EAM(const psPolynomial1D *myPoly,
    35                                        const psVector *x)
    36 {
    37     // PS_POLY_CHECK_NULL(myPoly, NULL);
    38     // PS_VECTOR_CHECK_NULL(x, NULL);
    39     // PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
    40 
    41     psVector *tmp;
    42 
    43     tmp = psVectorAlloc(x->n, PS_TYPE_F64);
    44     for (psS32 i=0;i<x->n;i++) {
    45         tmp->data.F64[i] = Polynomial1DEval_EAM(x->data.F64[i], myPoly);
    46     }
    47 
    48     return(tmp);
    49 }
    50 
    51 static void psPolynomial1DFree (psPolynomial1D *poly) {
    52 
    53   if (poly == NULL) return;
    54   psFree (poly->coeff);
    55   psFree (poly->coeffErr);
    56   psFree (poly->mask);
    57   return;
     19static void Polynomial1DFree_EAM (psPolynomial1D *poly) {
     20
     21    if (poly == NULL) return;
     22    psFree (poly->coeff);
     23    psFree (poly->coeffErr);
     24    psFree (poly->mask);
     25    return;
    5826}
    5927
    6028// XXX EAM : use Nterm = Norder + 1 definition 
    6129// XXX EAM : should we provide both order and nterms in struct?
    62 psPolynomial1D* Polynomial1DAlloc(psS32 nOrder,
    63                                   psPolynomialType type)
     30psPolynomial1D* Polynomial1DAlloc_EAM(psPolynomialType type,
     31                                      psS32 nOrder)
    6432{
    6533    // PS_INT_CHECK_NON_NEGATIVE(nOrder, NULL);
    … …  
    8048        newPoly->mask[i] = 0;
    8149    }
    82     psMemSetDeallocator(newPoly, (psFreeFunc) psPolynomial1DFree);
     50    psMemSetDeallocator(newPoly, (psFreeFunc) Polynomial1DFree_EAM);
    8351    return(newPoly);
    8452}
    8553
     54psF64 Polynomial1DEval_EAM(const psPolynomial1D* myPoly, psF64 x)
     55{
     56    psS32 loop_x = 0;
     57    psF64 polySum = 0.0;
     58    psF64 xSum = 1.0;
     59
     60    for (loop_x = 0; loop_x < myPoly->n + 1; loop_x++) {
     61        if (myPoly->mask[loop_x] == 0) {
     62            polySum += xSum * myPoly->coeff[loop_x];
     63        }
     64        xSum *= x;
     65    }
     66
     67    return(polySum);
     68}
     69
     70psVector *Polynomial1DEvalVector_EAM(const psPolynomial1D *myPoly,
     71                                     const psVector *x)
     72{
     73    // PS_POLY_CHECK_NULL(myPoly, NULL);
     74    // PS_VECTOR_CHECK_NULL(x, NULL);
     75    // PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
     76
     77    psVector *tmp;
     78
     79    tmp = psVectorAlloc(x->n, PS_TYPE_F64);
     80    for (psS32 i=0;i<x->n;i++) {
     81        tmp->data.F64[i] = Polynomial1DEval_EAM(myPoly, x->data.F64[i]);
     82    }
     83
     84    return(tmp);
     85}
     86
    8687// XXX EAM : my alternate BuildSums1D
    87 static psVector *BuildSums1D(psVector* sums,
    88                             psF64 x,
    89                             psS32 nTerm)
     88static psVector *BuildSums1D_EAM(psVector* sums,
     89                                psF64 x,
     90                                psS32 nTerm)
    9091{
    9192    psS32 nSum = 0;
    … …  
    109110
    110111// XXX EAM : test version of 1d fitting
    111 psPolynomial1D* VectorFitPolynomial1DOrd_EAM(psPolynomial1D* myPoly,
    112                                              psVector *mask,
    113                                              const psVector *x,
    114                                              const psVector *y,
    115                                              const psVector *yErr)
     112psPolynomial1D* VectorFitPolynomial1D_EAM(psPolynomial1D* myPoly,
     113                                          psVector *mask,
     114                                          psMaskType maskValue,
     115                                          const psVector *y,
     116                                          const psVector *yErr,
     117                                          const psVector *x)
    116118{
    117119    // I think this is 1 dimension down
    … …  
    130132    // XXX EAM : change from FILE to fd breaks this code:
    131133    if (psTraceGetLevel (".psLib.dataManip.VectorFitPolynomial1DOrd") >= 5) {
    132       FILE *f = psTraceGetDestination ();
    133       fprintf (f, "VectorFitPolynomial1D()\n");
    134       for (int i = 0; i < x->n; i++) {
    135         fprintf (f, "(x, y, yErr) is (%f, %f, %f)\n", x->data.F64[i], y->data.F64[i], yErr->data.F64[i]);
    136       }
     134        FILE *f = psTraceGetDestination ();
     135        fprintf (f, "VectorFitPolynomial1D()\n");
     136        for (int i = 0; i < x->n; i++) {
     137            fprintf (f, "(x, y, yErr) is (%f, %f, %f)\n", x->data.F64[i], y->data.F64[i], yErr->data.F64[i]);
     138        }
    137139    }
    138140# endif
    … …  
    155157    // Build the B and A data structs.
    156158    for (int k = 0; k < x->n; k++) {
    157         if ((mask != NULL) && mask->data.U8[k]) continue;
    158         xSums = BuildSums1D(xSums, x->data.F64[k], nTerm);
     159        if ((mask != NULL) && (mask->data.U8[k] & maskValue)) continue;
     160        xSums = BuildSums1D_EAM(xSums, x->data.F64[k], nTerm);
    159161     
    160162        if (yErr == NULL) {
    … …  
    190192    }
    191193    else
    192     // LUD version of the fit
     194        // LUD version of the fit
    193195    {
    194196        psImage *ALUD = NULL;
    … …  
    218220// ********************** 2D polynomial functions ******************
    219221
    220 // XXX EAM : this version uses myPoly->nX as Norder, not Nterms
    221 psVector *Polynomial2DEvalVector(const psPolynomial2D *myPoly,
    222                                  const psVector *x,
    223                                  const psVector *y)
    224 
    225 {
    226     // PS_POLY_CHECK_NULL(myPoly, NULL);
    227     // PS_VECTOR_CHECK_NULL(x, NULL);
    228     // PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
    229     // PS_VECTOR_CHECK_NULL(y, NULL);
    230     // PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
    231 
    232     psVector *tmp;
    233     psS32 vecLen=x->n;
    234 
    235     // Determine the length of the output vector to by the minimum of the x,y vectors
    236     if (y->n < vecLen) {
    237         vecLen = y->n;
    238     }
    239 
    240     // Create output vector to return
    241     tmp = psVectorAlloc(vecLen, PS_TYPE_F32);
    242 
    243     // Evaluate the polynomial at the specified points
    244     for (psS32 i=0; i<vecLen; i++) {
    245         tmp->data.F32[i] = Polynomial2DEval(myPoly,x->data.F32[i],y->data.F32[i]);
    246     }
    247 
    248     // Return output vector
    249     return(tmp);
    250 }
    251 
    252 // XXX EAM : this version uses the F64 vectors
    253 psVector *Polynomial2DEvalVectorD(const psPolynomial2D *myPoly,
    254                                   const psVector *x,
    255                                   const psVector *y)
    256 
    257 {
    258     // PS_POLY_CHECK_NULL(myPoly, NULL);
    259     // PS_VECTOR_CHECK_NULL(x, NULL);
    260     // PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
    261     // PS_VECTOR_CHECK_NULL(y, NULL);
    262     // PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
    263 
    264     psVector *tmp;
    265     psS32 vecLen=x->n;
    266 
    267     // Determine the length of the output vector to by the minimum of the x,y vectors
    268     if (y->n < vecLen) {
    269         vecLen = y->n;
    270     }
    271 
    272     // Create output vector to return
    273     tmp = psVectorAlloc(vecLen, PS_TYPE_F64);
    274 
    275     // Evaluate the polynomial at the specified points
    276     for (psS32 i=0; i<vecLen; i++) {
    277         tmp->data.F64[i] = Polynomial2DEval(myPoly,x->data.F64[i],y->data.F64[i]);
    278     }
    279 
    280     // Return output vector
    281     return(tmp);
    282 }
    283 
    284 static void psPolynomial2DFree(psPolynomial2D *poly) {
    285 
    286   if (poly == NULL) return;
    287   for (int i = 0; i < poly->nX + 1; i++) {
    288     psFree (poly->coeff[i]);
    289     psFree (poly->coeffErr[i]);
    290     psFree (poly->mask[i]);
    291   }
    292   psFree (poly->coeff);
    293   psFree (poly->coeffErr);
    294   psFree (poly->mask);
    295   return;
     222static void Polynomial2DFree_EAM(psPolynomial2D *poly) {
     223
     224    if (poly == NULL) return;
     225    for (int i = 0; i < poly->nX + 1; i++) {
     226        psFree (poly->coeff[i]);
     227        psFree (poly->coeffErr[i]);
     228        psFree (poly->mask[i]);
     229    }
     230    psFree (poly->coeff);
     231    psFree (poly->coeffErr);
     232    psFree (poly->mask);
     233    return;
    296234}
    297235
    298236// XXX EAM : use Nterm = Norder + 1 definition 
    299237// the user requests a polynomial of order Norder
    300 psPolynomial2D* Polynomial2DAlloc(psS32 nXorder, psS32 nYorder,
    301                                     psPolynomialType type)
     238psPolynomial2D* Polynomial2DAlloc_EAM(psPolynomialType type,
     239                                      psS32 nXorder,
     240                                      psS32 nYorder)
    302241{
    303242    // PS_INT_CHECK_NON_NEGATIVE(nXorder, NULL);
    … …  
    332271        }
    333272    }
    334     psMemSetDeallocator(newPoly, (psFreeFunc) psPolynomial2DFree);
     273    psMemSetDeallocator(newPoly, (psFreeFunc) Polynomial2DFree_EAM);
    335274    return(newPoly);
    336275}
    337276
    338 // XXX EAM : BuildSums2D in analogy with BuildSums1D
    339 static psImage *BuildSums2D(psImage* sums,
    340                             psF64 x,      psF64 y,
    341                             psS32 nXterm, psS32 nYterm)
    342 {
    343     psS32 nXsum = 0;
    344     psS32 nYsum = 0;
    345     psF64 xSum = 1.0;
    346     psF64 ySum = 1.0;
    347 
    348     nXsum = 2*nXterm;
    349     nYsum = 2*nYterm;
    350     if (sums == NULL) {
    351         sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
    352     }
    353     if ((nXsum != sums->numCols) || (nYsum != sums->numRows)) {
    354         psFree (sums);
    355         sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
    356     }
    357 
    358     ySum = 1.0;
    359     for (int j = 0; j < nYsum; j++) {
    360         xSum = ySum;
    361         for (int i = 0; i < nXsum; i++) {
    362             sums->data.F64[i][j] = xSum;
    363             xSum *= x;
    364         }
    365         ySum *= y;
    366     }
    367     return (sums);
    368 }
    369 
    370 // XXX EAM : test version of 2d fitting
    371 psPolynomial2D* VectorFitPolynomial2DOrd_EAM(psPolynomial2D* myPoly,
    372                                              psVector* mask,
    373                                              const psVector* x,
    374                                              const psVector* y,
    375                                              const psVector* z,
    376                                              const psVector* zErr)
    377 {
    378     // I think this is 1 dimension down
    379     psImage*      A = NULL;
    380     psVector*     B = NULL;
    381     psImage*   Sums = NULL;
    382     psF64 wt;
    383     psS32 nTerm, nXterm, nYterm;
    384 
    385     nXterm = myPoly->nX + 1;
    386     nYterm = myPoly->nY + 1;
    387     nTerm = nXterm * nYterm;
    388 
    389     A = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
    390     B = psVectorAlloc(nTerm, PS_TYPE_F64);
    391 
    392     // Initialize data structures (why is this not a function!)
    393     for (int i = 0; i < nTerm; i++) {
    394         B->data.F64[i] = 0.0;
    395         for (int j = 0; j < nTerm; j++) {
    396             A->data.F64[i][j] = 0.0;
    397         }
    398     }
    399 
    400     // Sums look like: 1, x, x^2, ... x^(2n+1), y, xy, x^2y, ... x^(2n+1)
    401 
    402     // Build the B and A data structs.
    403     for (int k  = 0; k < x->n; k++) {
    404       if ((mask != NULL) && mask->data.U8[k]) continue;
    405         Sums = BuildSums2D(Sums, x->data.F64[k], y->data.F64[k], nXterm, nYterm);
    406      
    407         if (zErr == NULL) {
    408             wt = 1.0;
    409         } else {
    410             // this should probably by zErr^2 !!
    411             // this should filter zErr == 0 values
    412             wt = 1.0 / zErr->data.F64[k];
    413         }
    414 
    415         // we could skip half of the array and assign at the end
    416         // we must handle masked orders
    417         for (int n = 0; n < nXterm; n++) {
    418             for (int m = 0; m < nYterm; m++) {
    419                 B->data.F64[n+m*nXterm] += z->data.F64[k] * Sums->data.F64[n][m] * wt;
    420             }
    421         }
    422 
    423         for (int i = 0; i < nXterm; i++) {
    424             for (int j = 0; j < nYterm; j++) {
    425                 for (int n = 0; n < nXterm; n++) {
    426                     for (int m = 0; m < nYterm; m++) {
    427                         A->data.F64[i+j*nXterm][n+m*nXterm] += Sums->data.F64[i+n][j+m] * wt;
    428                     }
    429                 }
    430             }
    431         }
    432     }
    433  
    434     // does the solution in place
    435     psGaussJordan (A, B);
    436    
    437     // XXX: How do we know if these routines were successful?
    438     // ALUD = psMatrixLUD(ALUD, &outPerm, A);
    439     // coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
    440 
    441     for (int n = 0; n < nXterm; n++) {
    442         for (int m = 0; m < nYterm; m++) {
    443             myPoly->coeff[n][m] = B->data.F64[n+m*nXterm];
    444         }
    445     }
    446 
    447     psFree(A);
    448     psFree(B);
    449     psFree(Sums);
    450 
    451     psTrace(".psLib.dataManip.VectorFitPolynomial2DOrd", 4,
    452             "---- VectorFitPolynomial2DOrd() begin ----\n");
    453     return (myPoly);
    454 }
    455 
    456277// write out the terms of the given 2D polynomial
    457 void psPolynomial2DDump (psPolynomial2D *poly) {
     278void Polynomial2DDump_EAM (psPolynomial2D *poly) {
    458279
    459280    for (int i = 0; i < poly->nX + 1; i++) {
    … …  
    464285}   
    465286
    466 psPolynomial2D* RobustFit2D_nomask(psPolynomial2D* poly,
    467                             const psVector* x,
    468                             const psVector* y,
    469                             const psVector* z,
    470                             const psVector* dz)
    471 {
    472     psVector *X;
    473     psVector *Y;
    474     psVector *Z;
    475     psVector *dZ;
    476 
    477     psVector *zFit   = NULL;
    478     psVector *zResid = NULL;
    479     psStats  *stats  = NULL;
    480 
    481     X  = psVectorCopy (NULL, x, PS_TYPE_F64);
    482     Y  = psVectorCopy (NULL, y, PS_TYPE_F64);
    483     Z  = psVectorCopy (NULL, z, PS_TYPE_F64);
    484     dZ = psVectorCopy (NULL, dz, PS_TYPE_F64);
    485 
    486     for (int N = 0; N < 3; N++) {
    487         // XXX EAM : this would be better defined with an element mask
    488         poly   = VectorFitPolynomial2DOrd_EAM (poly, NULL, X, Y, Z, dZ);
    489         zFit   = Polynomial2DEvalVectorD (poly, x, y);
    490         zResid = (psVector *) psBinaryOp (NULL, (void *) z, "-", (void *) zFit);
    491 
    492         stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
    493         stats  = psVectorStats (stats, zResid, NULL, NULL, 0);
    494         psTrace (".psphot.RobustFit", 4, "residual stats for robust fit:  %g +/- %g (%d pts)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues);
    495 
    496         // re-create X, Y, Z, dZ if pts are valid
    497         int n = 0;
    498         for (int i = 0; i < zResid->n; i++) {
    499             if (fabs(zResid->data.F64[i] - stats->clippedMean) > 3*stats->clippedStdev) {
    500               continue;
    501             }
    502             X->data.F64[n]  =  x->data.F64[i];
    503             Y->data.F64[n]  =  y->data.F64[i];
    504             Z->data.F64[n]  =  z->data.F64[i];
    505             dZ->data.F64[n] = dz->data.F64[i];
    506             n++;
    507         }
    508         X->n = n;
    509         Y->n = n;
    510         Z->n = n;
    511         dZ->n = n;
    512     }
    513     return (poly);
    514 }
    515 
    516 // XXX EAM : be careful here with F32 vs F64 vectors
    517 psPolynomial2D* RobustFit2D(psPolynomial2D* poly,
    518                             psVector* mask,
    519                             const psVector* x,
    520                             const psVector* y,
    521                             const psVector* z,
    522                             const psVector* dz)
    523 {
    524     // PS_VECTOR_CHECK_NULL(mask, NULL);
    525     // PS_VECTOR_CHECK_NULL(x, NULL);
    526     // PS_VECTOR_CHECK_NULL(y, NULL);
    527     // PS_VECTOR_CHECK_NULL(z, NULL);
    528     // PS_VECTOR_CHECK_NULL(dz, NULL);
    529 
    530     psVector *zFit   = NULL;
    531     psVector *zResid = psVectorAlloc (x->n, PS_TYPE_F64);
    532     psStats  *stats  = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
    533 
    534     for (int N = 0; N < 3; N++) {
    535         poly   = VectorFitPolynomial2DOrd_EAM (poly, mask, x, y, z, dz);
    536         zFit   = Polynomial2DEvalVectorD (poly, x, y);
    537         zResid = (psVector *) psBinaryOp (zResid, (void *) z, "-", (void *) zFit);
    538 
    539         stats  = psVectorStats (stats, zResid, NULL, mask, 1);
    540         psTrace (".psphot.RobustFit", 4, "residual stats for robust fit:  %g +/- %g\n",
    541                  stats->sampleMean, stats->sampleStdev);
    542 
    543         // set mask if pts are not valid
    544         // we are masking out any point which is out of range
    545         // recovery is not allowed with this scheme
    546         for (int i = 0; i < zResid->n; i++) {
    547           if (mask->data.U8[i]) continue;
    548           if (fabs(zResid->data.F64[i] - stats->sampleMean) > 3*stats->sampleStdev) {
    549             mask->data.U8[i] = 1;
    550             continue;
    551           }       
    552         }
    553         psFree (zFit);
    554     }
    555     psFree (zResid);
    556     psFree (stats);
    557     return (poly);
    558 }
    559 
    560 // XXX EAM : VectorFitPolynomial2DOrd and Polynomial2DEvalVector require different types (F32 vs F64)
    561 
    562 psF64 Polynomial2DEval(const psPolynomial2D* myPoly,
    563                        psF64 x,
    564                        psF64 y)
     287psF64 Polynomial2DEval_EAM(const psPolynomial2D* myPoly,
     288                           psF64 x,
     289                           psF64 y)
    565290{
    566291    // PS_POLY_CHECK_NULL(myPoly, NAN);
    … …  
    587312    return(polySum);
    588313}
     314
     315// XXX EAM : this version uses myPoly->nX as Norder, not Nterms
     316psVector *Polynomial2DEvalVector_EAM(const psPolynomial2D *myPoly,
     317                                     const psVector *x,
     318                                     const psVector *y)
     319
     320{
     321    // PS_POLY_CHECK_NULL(myPoly, NULL);
     322    // PS_VECTOR_CHECK_NULL(x, NULL);
     323    // PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
     324    // PS_VECTOR_CHECK_NULL(y, NULL);
     325    // PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
     326
     327    psVector *tmp;
     328    psS32 vecLen=x->n;
     329
     330    // Determine the length of the output vector to by the minimum of the x,y vectors
     331    if (y->n < vecLen) {
     332        vecLen = y->n;
     333    }
     334
     335    // Create output vector to return
     336    tmp = psVectorAlloc(vecLen, PS_TYPE_F32);
     337
     338    // Evaluate the polynomial at the specified points
     339    for (psS32 i=0; i<vecLen; i++) {
     340        tmp->data.F32[i] = Polynomial2DEval_EAM(myPoly, x->data.F32[i], y->data.F32[i]);
     341    }
     342
     343    // Return output vector
     344    return(tmp);
     345}
     346
     347// XXX EAM : this version uses the F64 vectors
     348psVector *Polynomial2DEvalVectorD_EAM(const psPolynomial2D *myPoly,
     349                                  const psVector *x,
     350                                  const psVector *y)
     351
     352{
     353    // PS_POLY_CHECK_NULL(myPoly, NULL);
     354    // PS_VECTOR_CHECK_NULL(x, NULL);
     355    // PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
     356    // PS_VECTOR_CHECK_NULL(y, NULL);
     357    // PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
     358
     359    psVector *tmp;
     360    psS32 vecLen=x->n;
     361
     362    // Determine the length of the output vector to by the minimum of the x,y vectors
     363    if (y->n < vecLen) {
     364        vecLen = y->n;
     365    }
     366
     367    // Create output vector to return
     368    tmp = psVectorAlloc(vecLen, PS_TYPE_F64);
     369
     370    // Evaluate the polynomial at the specified points
     371    for (psS32 i=0; i<vecLen; i++) {
     372        tmp->data.F64[i] = Polynomial2DEval_EAM(myPoly,x->data.F64[i],y->data.F64[i]);
     373    }
     374
     375    // Return output vector
     376    return(tmp);
     377}
     378
     379// XXX EAM : BuildSums2D in analogy with BuildSums1D
     380static psImage *BuildSums2D_EAM(psImage* sums,
     381                                psF64 x,      psF64 y,
     382                                psS32 nXterm, psS32 nYterm)
     383{
     384    psS32 nXsum = 0;
     385    psS32 nYsum = 0;
     386    psF64 xSum = 1.0;
     387    psF64 ySum = 1.0;
     388
     389    nXsum = 2*nXterm;
     390    nYsum = 2*nYterm;
     391    if (sums == NULL) {
     392        sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
     393    }
     394    if ((nXsum != sums->numCols) || (nYsum != sums->numRows)) {
     395        psFree (sums);
     396        sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
     397    }
     398
     399    ySum = 1.0;
     400    for (int j = 0; j < nYsum; j++) {
     401        xSum = ySum;
     402        for (int i = 0; i < nXsum; i++) {
     403            sums->data.F64[i][j] = xSum;
     404            xSum *= x;
     405        }
     406        ySum *= y;
     407    }
     408    return (sums);
     409}
     410
     411// XXX EAM : test version of 2d fitting
     412psPolynomial2D* VectorFitPolynomial2D_EAM(psPolynomial2D* myPoly,
     413                                             psVector* mask,
     414                                             psMaskType maskValue,
     415                                             const psVector* z,
     416                                             const psVector* zErr,
     417                                             const psVector* x,
     418                                             const psVector* y)
     419{
     420    // I think this is 1 dimension down
     421    psImage*      A = NULL;
     422    psVector*     B = NULL;
     423    psImage*   Sums = NULL;
     424    psF64 wt;
     425    psS32 nTerm, nXterm, nYterm;
     426
     427    nXterm = myPoly->nX + 1;
     428    nYterm = myPoly->nY + 1;
     429    nTerm = nXterm * nYterm;
     430
     431    A = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
     432    B = psVectorAlloc(nTerm, PS_TYPE_F64);
     433
     434    // Initialize data structures (why is this not a function!)
     435    for (int i = 0; i < nTerm; i++) {
     436        B->data.F64[i] = 0.0;
     437        for (int j = 0; j < nTerm; j++) {
     438            A->data.F64[i][j] = 0.0;
     439        }
     440    }
     441
     442    // Sums look like: 1, x, x^2, ... x^(2n+1), y, xy, x^2y, ... x^(2n+1)
     443
     444    // Build the B and A data structs.
     445    for (int k  = 0; k < x->n; k++) {
     446        if ((mask != NULL) && (mask->data.U8[k] & maskValue)) continue;
     447        Sums = BuildSums2D_EAM(Sums, x->data.F64[k], y->data.F64[k], nXterm, nYterm);
     448     
     449        if (zErr == NULL) {
     450            wt = 1.0;
     451        } else {
     452            // this should probably by zErr^2 !!
     453            // this should filter zErr == 0 values
     454            wt = 1.0 / zErr->data.F64[k];
     455        }
     456
     457        // we could skip half of the array and assign at the end
     458        // we must handle masked orders
     459        for (int n = 0; n < nXterm; n++) {
     460            for (int m = 0; m < nYterm; m++) {
     461                B->data.F64[n+m*nXterm] += z->data.F64[k] * Sums->data.F64[n][m] * wt;
     462            }
     463        }
     464
     465        for (int i = 0; i < nXterm; i++) {
     466            for (int j = 0; j < nYterm; j++) {
     467                for (int n = 0; n < nXterm; n++) {
     468                    for (int m = 0; m < nYterm; m++) {
     469                        A->data.F64[i+j*nXterm][n+m*nXterm] += Sums->data.F64[i+n][j+m] * wt;
     470                    }
     471                }
     472            }
     473        }
     474    }
     475 
     476    // does the solution in place
     477    psGaussJordan (A, B);
     478   
     479    // XXX: How do we know if these routines were successful?
     480    // ALUD = psMatrixLUD(ALUD, &outPerm, A);
     481    // coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
     482
     483    for (int n = 0; n < nXterm; n++) {
     484        for (int m = 0; m < nYterm; m++) {
     485            myPoly->coeff[n][m] = B->data.F64[n+m*nXterm];
     486        }
     487    }
     488
     489    psFree(A);
     490    psFree(B);
     491    psFree(Sums);
     492
     493    psTrace(".psLib.dataManip.VectorFitPolynomial2DOrd", 4,
     494            "---- VectorFitPolynomial2DOrd() begin ----\n");
     495    return (myPoly);
     496}
     497
     498// XXX EAM : be careful here with F32 vs F64 vectors
     499psPolynomial2D* VectorClipFitPolynomial2D_EAM(psPolynomial2D* poly,
     500                                              psStats *stats,
     501                                              psVector* mask,
     502                                              psMaskType maskValue,
     503                                              const psVector* z,
     504                                              const psVector* dz,
     505                                              const psVector* x,
     506                                              const psVector* y)
     507{
     508    psVector *zFit   = NULL;
     509    psVector *zResid = psVectorAlloc (x->n, PS_TYPE_F64);
     510
     511    // XXX EAM : use SAMPLE_MEAN and SAMPLE_STDEV for stats:
     512    stats->options |= (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
     513
     514    for (int N = 0; N < 3; N++) {
     515        poly   = VectorFitPolynomial2D_EAM (poly, mask, maskValue, z, dz, x, y);
     516        zFit   = Polynomial2DEvalVectorD_EAM (poly, x, y);
     517        zResid = (psVector *) psBinaryOp (zResid, (void *) z, "-", (void *) zFit);
     518
     519        stats  = psVectorStats (stats, zResid, NULL, mask, maskValue);
     520        psTrace (".psphot.RobustFit", 4, "residual stats for robust fit:  %g +/- %g\n",
     521                 stats->sampleMean, stats->sampleStdev);
     522
     523        // set mask if pts are not valid
     524        // we are masking out any point which is out of range
     525        // recovery is not allowed with this scheme
     526        for (int i = 0; i < zResid->n; i++) {
     527            if (mask->data.U8[i]) continue;
     528            if (fabs(zResid->data.F64[i] - stats->sampleMean) > 3*stats->sampleStdev) {
     529                mask->data.U8[i] &= 0x01;
     530                continue;
     531            }       
     532        }
     533        psFree (zFit);
     534    }
     535    psFree (zResid);
     536    return (poly);
     537}
     538
     539// XXX EAM : VectorFitPolynomial2D and Polynomial2DEvalVector require different types (F32 vs F64)
Note: See TracChangeset for help on using the changeset viewer.