IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 15841


Ignore:
Timestamp:
Dec 14, 2007, 3:20:03 PM (18 years ago)
Author:
Paul Price
Message:

A little bit of cleaning up, and making it so that the statistic to be used may be specified.

Location:
trunk/psLib/src
Files:
1 added
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/imageops/Makefile.am

    r14929 r15841  
    3434        psImageBinning.h \
    3535        psImageMap.h \
     36        psImageMapFit.h \
    3637        psImagePixelInterpolate.h \
    3738        psImageUnbin.h

  • trunk/psLib/src/imageops/psImageMap.c

    r15599 r15841  
    77 *  @author Eugene Magnier, IfA
    88 *
    9  *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2007-11-13 18:28:02 $
     9 *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2007-12-15 01:20:03 $
    1111 *
    1212 *  Copyright 2007 Institute for Astronomy, University of Hawaii
    … …  
    3232#include "psStats.h"
    3333#include "psImageBinning.h"
    34 #include "psImageMap.h"
    3534#include "psImagePixelInterpolate.h"
    3635#include "psImageUnbin.h"
    3736
     37#include "psImageMap.h"
     38
    3839static void psImageMapFree(psImageMap *map) {
    3940
    4041    if (!map) return;
    4142
    42     psFree (map->map);
    43     psFree (map->error);
    44     psFree (map->stats);
    45     psFree (map->binning);
     43    psFree(map->map);
     44    psFree(map->error);
     45    psFree(map->stats);
     46    psFree(map->binning);
    4647
    4748    return;
    … …  
    103104}
    104105
    105 bool psImageMapModifyScale(psImageMap *map, int nXruff, int nYruff) {
    106 
     106bool psImageMapModifyScale(psImageMap *map, int nXruff, int nYruff)
     107{
    107108    assert (map);
    108109
    … …  
    121122// generate a psImageMap (or NULL) with the given number of superpixels in X and Y
    122123// this function returns an error if the output map has impossible holes
    123 bool psImageMapGenerate (psImageMap *map, psVector *x, psVector *y, psVector *f, psVector *df, float badFrac) {
     124bool psImageMapGenerate(psImageMap *map, const psVector *x, const psVector *y,
     125                        const psVector *f, const psVector *df, float badFrac)
     126{
     127    // XXX asserts
    124128
    125129    psImage *mask = psImageAlloc (map->map->numCols, map->map->numRows, PS_TYPE_MASK);
    … …  
    253257
    254258// using the points given, generate a map with maximum resolution that yields only good and ok pixels
    255 bool psImageMapGenerateScale (psImageMap *map, psVector *x, psVector *y, psVector *f, psVector *df, float badFrac) {
     259bool psImageMapGenerateScale(psImageMap *map, const psVector *x, const psVector *y,
     260                             const psVector *f, const psVector *df, float badFrac)
     261{
     262    // XXX Asserts
    256263
    257264    int nXruff, nYruff;
    … …  
    297304
    298305// x,y are in fractional pixel coords of the fine image (pixel center: 0.5)
    299 double psImageMapEval (psImageMap *map, float x, float y) {
     306double psImageMapEval(const psImageMap *map, float x, float y) {
    300307
    301308    double result;
    … …  
    306313}
    307314
    308 psVector *psImageMapEvalVector (psImageMap *map, psVector *x, psVector *y) {
     315psVector *psImageMapEvalVector(const psImageMap *map, const psVector *x, const psVector *y) {
    309316
    310317    assert (x);

  • trunk/psLib/src/imageops/psImageMap.h

    r15598 r15841  
    77 *  @author Eugene Magnier, IfA
    88 *
    9  *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2007-11-13 18:23:17 $
     9 *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2007-12-15 01:20:03 $
    1111 *
    1212 *  Copyright 2007 Institute for Astronomy, University of Hawaii
    … …  
    1818/// @addtogroup ImageOps Image Operations
    1919/// @{
     20
     21#include <psStats.h>
     22#include <psImage.h>
     23#include <psImageBinning.h>
     24#include <psVector.h>
    2025
    2126// a structure to describe the 2D variations of some quantity as a function of position the
    … …  
    3338    int nPoor;
    3439    int nGood;
     40    psStatsOptions singleMean, singleStdev;  // Statistics for mean and stdev when there's a single pixel
    3541} psImageMap;
    3642
    … …  
    4349
    4450// generate a psImageMap (or NULL) with the given number of superpixels in X and Y
    45 bool psImageMapGenerate (psImageMap *map, psVector *x, psVector *y, psVector *f, psVector *df, float badFrac);
     51bool psImageMapGenerate (psImageMap *map, const psVector *x, const psVector *y, const psVector *f, const psVector *df, float badFrac);
    4652
    47 bool psImageMapGenerateScale (psImageMap *map, psVector *x, psVector *y, psVector *f, psVector *df, float badFrac);
     53bool psImageMapGenerateScale (psImageMap *map, const psVector *x, const psVector *y, const psVector *f, const psVector *df, float badFrac);
    4854
    4955// apply the psImageMap to the given coordinate (fine image pixels)
    50 double psImageMapEval (psImageMap *map, float x, float y);
     56double psImageMapEval (const psImageMap *map, float x, float y);
    5157
    5258// apply the psImageMap to the given coordinate vectors (fine image pixels)
    53 psVector *psImageMapEvalVector (psImageMap *map, psVector *x, psVector *y);
    54 
    55 // fit the image map to a set of points
    56 bool psImageMapFit (psImageMap *map, psVector *mask, psMaskType maskValue, psVector *x, psVector *y, psVector *f, psVector *df);
    57 
    58 // fit the image map to a set of points
    59 bool psImageMapClipFit (psImageMap *map, psStats *stats, psVector *mask, psMaskType maskValue, psVector *x, psVector *y, psVector *f, psVector *df);
    60 
    61 bool psImageMapFit1DinY (psImageMap *map, psVector *mask, psMaskType maskValue, psVector *x, psVector *y, psVector *f, psVector *df);
    62 bool psImageMapFit1DinX (psImageMap *map, psVector *mask, psMaskType maskValue, psVector *x, psVector *y, psVector *f, psVector *df);
     59psVector *psImageMapEvalVector (const psImageMap *map, const psVector *x, const psVector *y);
    6360
    6461/// @}

  • trunk/psLib/src/imageops/psImageMapFit.c

    r15041 r15841  
    77 *  @author Eugene Magnier, IfA
    88 *
    9  *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2007-09-27 04:27:03 $
     9 *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2007-12-15 01:20:03 $
    1111 *
    1212 *  Copyright 2007 Institute for Astronomy, University of Hawaii
    … …  
    3737// #include "psImageUnbin.h"
    3838
     39#include "psImageMapFit.h"
     40
     41
    3942// given a randomly-sampled field of values & weights at points: (f, df) @ (x, y), find the
    4043// best fit image from which Bilinear interpolation yields the input field.  The fitted image
    … …  
    4447
    4548// map defines the output image dimensions and scaling.
    46 bool psImageMapFit (psImageMap *map, psVector *mask, psMaskType maskValue, psVector *x, psVector *y, psVector *f, psVector *df) {
    47 
    48     int I, J;
    49 
     49bool psImageMapFit(psImageMap *map, const psVector *mask, psMaskType maskValue,
     50                   const psVector *x, const psVector *y, const psVector *f, const psVector *df)
     51{
    5052    // XXX Add Asserts
    5153
    … …  
    5658    // no spatial information, just calculate mean & stdev
    5759    if ((Nx == 1) && (Ny == 1)) {
    58         psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     60        psStatsInit(map->stats);
     61
     62        // the user has supplied one of various stats option pairs,
     63        psStatsOptions mean = psStatsMeanOption(map->stats->options);
     64        psStatsOptions stdev = psStatsStdevOption(map->stats->options);
     65        if (!psStatsSingleOption(mean)) {
     66            psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected");
     67            return false;
     68        }
     69        if (!psStatsSingleOption(stdev)) {
     70            psError(PS_ERR_UNKNOWN, true, "no valid stdev stats option selected");
     71            return false;
     72        }
    5973
    6074        // XXX does ROBUST_MEDIAN work with weight?
    61         psVectorStats (stats, f, NULL, mask, maskValue);
    62 
    63         map->map->data.F32[0][0]   = stats->robustMedian;
    64         map->error->data.F32[0][0] = stats->robustStdev;
    65         psFree (stats);
     75        psVectorStats(map->stats, f, NULL, mask, maskValue);
     76
     77        map->map->data.F32[0][0]   = psStatsGetValue(map->stats, mean);
     78        map->error->data.F32[0][0] = psStatsGetValue(map->stats, stdev);
    6679        return true;
    6780    }
    6881
    6982    if (Nx == 1) {
    70         bool status;
    71         status = psImageMapFit1DinY (map, mask, maskValue, x, y, f, df);
    72         return status;
     83        bool status;
     84        status = psImageMapFit1DinY (map, mask, maskValue, x, y, f, df);
     85        return status;
    7386    }
    7487    if (Ny == 1) {
    75         bool status;
    76         status = psImageMapFit1DinX (map, mask, maskValue, x, y, f, df);
    77         return status;
     88        bool status;
     89        status = psImageMapFit1DinX (map, mask, maskValue, x, y, f, df);
     90        return status;
    7891    }
    7992
    8093    // set up the redirection table so we can use sA[-1][-1], etc
     94    // XXX psKernel does this for you --- PAP.
    8195    float SAm[3][3], *SAv[3], **sA;
    8296    // float TAm[3][3], *TAv[3], **tA;
    … …  
    182196                if (m == Ny - 1) Qy = 0;
    183197
    184                 assert (isfinite(fi));
    185                 assert (isfinite(wt));
    186                 assert (isfinite(rx));
    187                 assert (isfinite(ry));
     198                assert (isfinite(fi));
     199                assert (isfinite(wt));
     200                assert (isfinite(rx));
     201                assert (isfinite(ry));
    188202
    189203                // points at offset 1,1
    … …  
    234248
    235249            // I[ 0][ 0] = index for this n,m element:
    236             I = n + Nx * m;
     250            int I = n + Nx * m;
    237251            B->data.F32[I] = fi_rx_ry + fi_rx_py + fi_px_ry + fi_px_py;
    238252
    … …  
    245259                    if (m + jm <   0) continue;
    246260                    if (m + jm >= Ny) continue;
    247                     J = (n + jn) + Nx * (m + jm);
     261                    int J = (n + jn) + Nx * (m + jm);
    248262                    A->data.F32[J][I] = sA[jn][jm];
    249263                    // fprintf (stderr, "A %d %d (%d %d : %d %d): %f\n", I, J, n, m, n + jn, m + jm, sA[jn][jm]);
    … …  
    306320    for (int n = 0; n < Nx; n++) {
    307321        for (int m = 0; m < Ny; m++) {
    308             I = n + Nx * m;
     322            int I = n + Nx * m;
    309323            map->map->data.F32[m][n] = B->data.F32[I];
    310324            map->error->data.F32[m][n] = sqrt(A->data.F32[I][I]);
    … …  
    320334
    321335// measure residuals on each pass and clip outliers based on stats
    322 bool psImageMapClipFit (psImageMap *map, psStats *stats, psVector *inMask, psMaskType maskValue, psVector *x, psVector *y, psVector *f, psVector *df) {
    323 
     336bool psImageMapClipFit(psImageMap *map, psStats *stats, psVector *inMask, psMaskType maskValue,
     337                       const psVector *x, const psVector *y, const psVector *f, const psVector *df)
     338{
    324339    // XXX add in full PS_ASSERTS
    325     assert (map);
    326     assert (stats);
    327     assert (x);
    328     assert (y);
    329     assert (f);
     340    assert(map);
     341    assert(stats);
     342    assert(x);
     343    assert(y);
     344    assert(f);
    330345
    331346    // the user supplies one of various stats option pairs,
    332347    // determine the desired mean and stdev STATS options:
    333     psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_MEAN_V3 | PS_STAT_FITTED_MEAN_V4);
    334     psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV_V2 | PS_STAT_FITTED_STDEV_V3 | PS_STAT_FITTED_STDEV_V4);
    335     if (!meanOption) {
     348    psStatsOptions meanOption = psStatsMeanOption(stats->options);
     349    psStatsOptions stdevOption = psStatsStdevOption(stats->options);
     350    if (!psStatsSingleOption(meanOption)) {
    336351        psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected");
    337352        return false;
    338353    }
    339     if (!stdevOption) {
     354    if (!psStatsSingleOption(stdevOption)) {
    340355        psError(PS_ERR_UNKNOWN, true, "no valid stdev stats option selected");
    341356        return false;
    … …  
    436451
    437452// map defines the output image dimensions and scaling.
    438 bool psImageMapFit1DinY (psImageMap *map, psVector *mask, psMaskType maskValue, psVector *x, psVector *y, psVector *f, psVector *df) {
    439 
    440     int I, J;
    441 
     453bool psImageMapFit1DinY(psImageMap *map, const psVector *mask, psMaskType maskValue,
     454                        const psVector *x, const psVector *y, const psVector *f, const psVector *df)
     455{
    442456    // XXX Add Asserts
    443457    assert (map->binning->nXruff == 1);
    … …  
    459473
    460474    for (int m = 0; m < Ny; m++) {
    461         // define & init summing variables
    462         float ry_ry = 0;
    463         float dy_ry = 0;
    464         float fi_ry = 0;
    465         float py_py = 0;
    466         float qy_py = 0;
    467         float fi_py = 0;
    468 
    469         // generate the sums for the fitting matrix element I,J
    470         // I = m
    471         // J = m + jm
    472         for (int i = 0; i < y->n; i++) {
    473 
    474             if (mask && (mask->data.U8[i] & maskValue)) continue;
    475 
    476             float dy = psImageBinningGetRuffY (map->binning, y->data.F32[i]) - (m + 0.5);
    477 
    478             bool edgeY = false;
    479             edgeY |= ((m == 1) && (dy < -1.0));
    480             edgeY |= ((m == Ny - 2) && (dy > +1.0));
    481 
    482             // skip points outside of 2x2 grid centered on n,m:
    483             if (!edgeY && (fabs(dy) > 1.0)) continue;
    484 
    485             // related offset values
    486             float ry = 1.0 - dy;
    487             float py = 1.0 + dy;
    488             float qy = -dy;
    489 
    490             // data value & weight for this point
    491             float fi = f->data.F32[i];
    492             float wt = 1.0;
    493             if (df != NULL) {
    494                 if (df->data.F32[i] == 0.0) {
    495                     wt = 0.0;
    496                 } else {
    497                     wt = 1.0 / PS_SQR(df->data.F32[i]); // XXX test for dz == NULL or dz_i = 0
    498                 }
    499             }
    500 
    501             // sum the appropriate elements for the different quadrants
    502             int Qy = (dy >= 0) ? 1 : 0;
    503             if (m ==      0) Qy = 1;
    504             if (m == Ny - 1) Qy = 0;
    505 
    506             assert (isfinite(fi));
    507             assert (isfinite(wt));
    508             assert (isfinite(ry));
    509 
    510             // points at offset 1,1
    511             if (Qy == 1) {
    512                 ry_ry += ry*ry*wt;
    513                 dy_ry += dy*ry*wt;
    514                 fi_ry += fi*ry*wt;
    515             }
    516             // points at offset 1,0
    517             if (Qy == 0) {
    518                 py_py += py*py*wt;
    519                 qy_py += qy*py*wt;
    520                 fi_py += fi*py*wt;
    521             }
    522         }
    523 
    524         // the chi-square derivatives have elements of the form g(n+jn,m+jm)*A(jn,jm),
    525         // jn,jm = -1 to +1. Convert the sums above into the correct coefficients
    526         sA[-1] = qy_py;
    527         sA[ 0] = ry_ry + py_py;
    528         sA[+1] = dy_ry;
    529 
    530         // I[ 0][ 0] = index for this n,m element:
    531         I = m;
    532         B->data.F32[I] = fi_ry + fi_py;
    533 
    534         // insert these values into their corresponding locations in A, B
    535         for (int jm = -1; jm <= +1; jm++) {
    536             if (m + jm <   0) continue;
    537             if (m + jm >= Ny) continue;
    538             J = (m + jm);
    539             A->data.F32[J][I] = sA[jm];
    540         }
     475        // define & init summing variables
     476        float ry_ry = 0;
     477        float dy_ry = 0;
     478        float fi_ry = 0;
     479        float py_py = 0;
     480        float qy_py = 0;
     481        float fi_py = 0;
     482
     483        // generate the sums for the fitting matrix element I,J
     484        // I = m
     485        // J = m + jm
     486        for (int i = 0; i < y->n; i++) {
     487
     488            if (mask && (mask->data.U8[i] & maskValue)) continue;
     489
     490            float dy = psImageBinningGetRuffY (map->binning, y->data.F32[i]) - (m + 0.5);
     491
     492            bool edgeY = false;
     493            edgeY |= ((m == 1) && (dy < -1.0));
     494            edgeY |= ((m == Ny - 2) && (dy > +1.0));
     495
     496            // skip points outside of 2x2 grid centered on n,m:
     497            if (!edgeY && (fabs(dy) > 1.0)) continue;
     498
     499            // related offset values
     500            float ry = 1.0 - dy;
     501            float py = 1.0 + dy;
     502            float qy = -dy;
     503
     504            // data value & weight for this point
     505            float fi = f->data.F32[i];
     506            float wt = 1.0;
     507            if (df != NULL) {
     508                if (df->data.F32[i] == 0.0) {
     509                    wt = 0.0;
     510                } else {
     511                    wt = 1.0 / PS_SQR(df->data.F32[i]); // XXX test for dz == NULL or dz_i = 0
     512                }
     513            }
     514
     515            // sum the appropriate elements for the different quadrants
     516            int Qy = (dy >= 0) ? 1 : 0;
     517            if (m ==      0) Qy = 1;
     518            if (m == Ny - 1) Qy = 0;
     519
     520            assert (isfinite(fi));
     521            assert (isfinite(wt));
     522            assert (isfinite(ry));
     523
     524            // points at offset 1,1
     525            if (Qy == 1) {
     526                ry_ry += ry*ry*wt;
     527                dy_ry += dy*ry*wt;
     528                fi_ry += fi*ry*wt;
     529            }
     530            // points at offset 1,0
     531            if (Qy == 0) {
     532                py_py += py*py*wt;
     533                qy_py += qy*py*wt;
     534                fi_py += fi*py*wt;
     535            }
     536        }
     537
     538        // the chi-square derivatives have elements of the form g(n+jn,m+jm)*A(jn,jm),
     539        // jn,jm = -1 to +1. Convert the sums above into the correct coefficients
     540        sA[-1] = qy_py;
     541        sA[ 0] = ry_ry + py_py;
     542        sA[+1] = dy_ry;
     543
     544        // I[ 0][ 0] = index for this n,m element:
     545        int I = m;
     546        B->data.F32[I] = fi_ry + fi_py;
     547
     548        // insert these values into their corresponding locations in A, B
     549        for (int jm = -1; jm <= +1; jm++) {
     550            if (m + jm <   0) continue;
     551            if (m + jm >= Ny) continue;
     552            int J = (m + jm);
     553            A->data.F32[J][I] = sA[jm];
     554        }
    541555    }
    542556
    … …  
    545559    psVectorInit (Empty, 0);
    546560    for (int i = 0; i < Ny; i++) {
    547         if (A->data.F32[i][i] == 0.0) {
    548             Empty->data.S8[i] = 1;
    549             for (int j = 0; j < Ny; j++) {
    550                 A->data.F32[i][j] = 0.0;
    551                 A->data.F32[j][i] = 0.0;
    552             }
    553             A->data.F32[i][i] = 1.0;
    554             B->data.F32[i] = 0.0;
    555         }
     561        if (A->data.F32[i][i] == 0.0) {
     562            Empty->data.S8[i] = 1;
     563            for (int j = 0; j < Ny; j++) {
     564                A->data.F32[i][j] = 0.0;
     565                A->data.F32[j][i] = 0.0;
     566            }
     567            A->data.F32[i][i] = 1.0;
     568            B->data.F32[i] = 0.0;
     569        }
    556570    }
    557571
    558572    if (!psMatrixGJSolveF32(A, B)) {
    559         psAbort ("failed on linear equations");
    560         psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
    561         psFree (A);
    562         psFree (B);
    563         return false;
     573        psAbort ("failed on linear equations");
     574        psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
     575        psFree (A);
     576        psFree (B);
     577        return false;
    564578    }
    565579
    566580    // set bad values to NaN
    567581    for (int i = 0; i < Ny; i++) {
    568         if (Empty->data.S8[i]) {
    569             B->data.F32[i] = NAN;
    570         }
     582        if (Empty->data.S8[i]) {
     583            B->data.F32[i] = NAN;
     584        }
    571585    }
    572586
    573587    for (int m = 0; m < Ny; m++) {
    574         map->map->data.F32[m][0] = B->data.F32[m];
    575         map->error->data.F32[m][0] = sqrt(A->data.F32[m][m]);
     588        map->map->data.F32[m][0] = B->data.F32[m];
     589        map->error->data.F32[m][0] = sqrt(A->data.F32[m][m]);
    576590    }
    577591
    … …  
    584598
    585599// map defines the output image dimensions and scaling.
    586 bool psImageMapFit1DinX (psImageMap *map, psVector *mask, psMaskType maskValue, psVector *x, psVector *y, psVector *f, psVector *df) {
    587 
    588     int I, J;
    589 
     600bool psImageMapFit1DinX(psImageMap *map, const psVector *mask, psMaskType maskValue,
     601                        const psVector *x, const psVector *y, const psVector *f, const psVector *df)
     602{
    590603    // XXX Add Asserts
    591604    assert (map->binning->nXruff == 1);
    … …  
    607620
    608621    for (int m = 0; m < Nx; m++) {
    609         // define & init summing variables
    610         float rx_rx = 0;
    611         float dx_rx = 0;
    612         float fi_rx = 0;
    613         float px_px = 0;
    614         float qx_px = 0;
    615         float fi_px = 0;
    616 
    617         // generate the sums for the fitting matrix element I,J
    618         // I = m
    619         // J = m + jm
    620         for (int i = 0; i < x->n; i++) {
    621 
    622             if (mask && (mask->data.U8[i] & maskValue)) continue;
    623 
    624             float dx = psImageBinningGetRuffX (map->binning, x->data.F32[i]) - (m + 0.5);
    625 
    626             bool edgeX = false;
    627             edgeX |= ((m == 1) && (dx < -1.0));
    628             edgeX |= ((m == Nx - 2) && (dx > +1.0));
    629 
    630             // skip points outside of 2x2 grid centered on n,m:
    631             if (!edgeX && (fabs(dx) > 1.0)) continue;
    632 
    633             // related offset values
    634             float rx = 1.0 - dx;
    635             float px = 1.0 + dx;
    636             float qx = -dx;
    637 
    638             // data value & weight for this point
    639             float fi = f->data.F32[i];
    640             float wt = 1.0;
    641             if (df != NULL) {
    642                 if (df->data.F32[i] == 0.0) {
    643                     wt = 0.0;
    644                 } else {
    645                     wt = 1.0 / PS_SQR(df->data.F32[i]); // XXX test for dz == NULL or dz_i = 0
    646                 }
    647             }
    648 
    649             // sum the appropriate elements for the different quadrants
    650             int Qx = (dx >= 0) ? 1 : 0;
    651             if (m ==      0) Qx = 1;
    652             if (m == Nx - 1) Qx = 0;
    653 
    654             assert (isfinite(fi));
    655             assert (isfinite(wt));
    656             assert (isfinite(rx));
    657 
    658             // points at offset 1,1
    659             if (Qx == 1) {
    660                 rx_rx += rx*rx*wt;
    661                 dx_rx += dx*rx*wt;
    662                 fi_rx += fi*rx*wt;
    663             }
    664             // points at offset 1,0
    665             if (Qx == 0) {
    666                 px_px += px*px*wt;
    667                 qx_px += qx*px*wt;
    668                 fi_px += fi*px*wt;
    669             }
    670         }
    671 
    672         // the chi-square derivatives have elements of the form g(n+jn,m+jm)*A(jn,jm),
    673         // jn,jm = -1 to +1. Convert the sums above into the correct coefficients
    674         sA[-1] = qx_px;
    675         sA[ 0] = rx_rx + px_px;
    676         sA[+1] = dx_rx;
    677 
    678         // I[ 0][ 0] = index for this n,m element:
    679         I = m;
    680         B->data.F32[I] = fi_rx + fi_px;
    681 
    682         // insert these values into their corresponding locations in A, B
    683         for (int jm = -1; jm <= +1; jm++) {
    684             if (m + jm <   0) continue;
    685             if (m + jm >= Nx) continue;
    686             J = (m + jm);
    687             A->data.F32[J][I] = sA[jm];
    688         }
     622        // define & init summing variables
     623        float rx_rx = 0;
     624        float dx_rx = 0;
     625        float fi_rx = 0;
     626        float px_px = 0;
     627        float qx_px = 0;
     628        float fi_px = 0;
     629
     630        // generate the sums for the fitting matrix element I,J
     631        // I = m
     632        // J = m + jm
     633        for (int i = 0; i < x->n; i++) {
     634
     635            if (mask && (mask->data.U8[i] & maskValue)) continue;
     636
     637            float dx = psImageBinningGetRuffX (map->binning, x->data.F32[i]) - (m + 0.5);
     638
     639            bool edgeX = false;
     640            edgeX |= ((m == 1) && (dx < -1.0));
     641            edgeX |= ((m == Nx - 2) && (dx > +1.0));
     642
     643            // skip points outside of 2x2 grid centered on n,m:
     644            if (!edgeX && (fabs(dx) > 1.0)) continue;
     645
     646            // related offset values
     647            float rx = 1.0 - dx;
     648            float px = 1.0 + dx;
     649            float qx = -dx;
     650
     651            // data value & weight for this point
     652            float fi = f->data.F32[i];
     653            float wt = 1.0;
     654            if (df != NULL) {
     655                if (df->data.F32[i] == 0.0) {
     656                    wt = 0.0;
     657                } else {
     658                    wt = 1.0 / PS_SQR(df->data.F32[i]); // XXX test for dz == NULL or dz_i = 0
     659                }
     660            }
     661
     662            // sum the appropriate elements for the different quadrants
     663            int Qx = (dx >= 0) ? 1 : 0;
     664            if (m ==      0) Qx = 1;
     665            if (m == Nx - 1) Qx = 0;
     666
     667            assert (isfinite(fi));
     668            assert (isfinite(wt));
     669            assert (isfinite(rx));
     670
     671            // points at offset 1,1
     672            if (Qx == 1) {
     673                rx_rx += rx*rx*wt;
     674                dx_rx += dx*rx*wt;
     675                fi_rx += fi*rx*wt;
     676            }
     677            // points at offset 1,0
     678            if (Qx == 0) {
     679                px_px += px*px*wt;
     680                qx_px += qx*px*wt;
     681                fi_px += fi*px*wt;
     682            }
     683        }
     684
     685        // the chi-square derivatives have elements of the form g(n+jn,m+jm)*A(jn,jm),
     686        // jn,jm = -1 to +1. Convert the sums above into the correct coefficients
     687        sA[-1] = qx_px;
     688        sA[ 0] = rx_rx + px_px;
     689        sA[+1] = dx_rx;
     690
     691        // I[ 0][ 0] = index for this n,m element:
     692        int I = m;
     693        B->data.F32[I] = fi_rx + fi_px;
     694
     695        // insert these values into their corresponding locations in A, B
     696        for (int jm = -1; jm <= +1; jm++) {
     697            if (m + jm <   0) continue;
     698            if (m + jm >= Nx) continue;
     699            int J = (m + jm);
     700            A->data.F32[J][I] = sA[jm];
     701        }
    689702    }
    690703
    … …  
    693706    psVectorInit (Empty, 0);
    694707    for (int i = 0; i < Nx; i++) {
    695         if (A->data.F32[i][i] == 0.0) {
    696             Empty->data.S8[i] = 1;
    697             for (int j = 0; j < Nx; j++) {
    698                 A->data.F32[i][j] = 0.0;
    699                 A->data.F32[j][i] = 0.0;
    700             }
    701             A->data.F32[i][i] = 1.0;
    702             B->data.F32[i] = 0.0;
    703         }
     708        if (A->data.F32[i][i] == 0.0) {
     709            Empty->data.S8[i] = 1;
     710            for (int j = 0; j < Nx; j++) {
     711                A->data.F32[i][j] = 0.0;
     712                A->data.F32[j][i] = 0.0;
     713            }
     714            A->data.F32[i][i] = 1.0;
     715            B->data.F32[i] = 0.0;
     716        }
    704717    }
    705718
    706719    if (!psMatrixGJSolveF32(A, B)) {
    707         psAbort ("failed on linear equations");
    708         psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
    709         psFree (A);
    710         psFree (B);
    711         return false;
     720        psAbort ("failed on linear equations");
     721        psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
     722        psFree (A);
     723        psFree (B);
     724        return false;
    712725    }
    713726
    714727    // set bad values to NaN
    715728    for (int i = 0; i < Nx; i++) {
    716         if (Empty->data.S8[i]) {
    717             B->data.F32[i] = NAN;
    718         }
     729        if (Empty->data.S8[i]) {
     730            B->data.F32[i] = NAN;
     731        }
    719732    }
    720733
    721734    for (int m = 0; m < Nx; m++) {
    722         map->map->data.F32[0][m] = B->data.F32[I];
    723         map->error->data.F32[0][m] = sqrt(A->data.F32[I][I]);
     735        int I = m; // XXX I'm not entirely sure about this; it wasn't set for this scope --- PAP.
     736        map->map->data.F32[0][m] = B->data.F32[I];
     737        map->error->data.F32[0][m] = sqrt(A->data.F32[I][I]);
    724738    }
    725739

  • trunk/psLib/src/pslib_strict.h

    r15818 r15841  
    99*  @author Eric Van Alst, MHPCC
    1010*
    11 *  @version $Revision: 1.32 $ $Name: not supported by cvs2svn $
    12 *  @date $Date: 2007-12-14 00:40:15 $
     11*  @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
     12*  @date $Date: 2007-12-15 01:20:03 $
    1313*
    1414*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    5858#include "psImageUnbin.h"
    5959#include "psImageMap.h"
     60#include "psImageMapFit.h"
    6061
    6162#include "psImageJpeg.h"
Note: See TracChangeset for help on using the changeset viewer.