IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 14983


Ignore:
Timestamp:
Sep 21, 2007, 5:05:50 PM (19 years ago)
Author:
Paul Price
Message:

Including psMemory.h where required.

Location:
trunk/psLib/src
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r14923 r14983  
    88 *  @author Eugene Magnier, IfA
    99 *
    10  *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2007-09-20 23:53:46 $
     10 *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2007-09-22 03:05:50 $
    1212 *
    1313 *  Copyright 2007 Institute for Astronomy, University of Hawaii
    … …  
    1919
    2020#include <stdio.h>
     21#include "psMemory.h"
    2122#include "psError.h"
    2223#include "psAbort.h"
    … …  
    3839
    3940void psImageBinningSetRuffSize(psImageBinning *binning, psImageBinningAlign align) {
    40    
     41
    4142    assert (binning->nXfine > 0);
    4243    assert (binning->nYfine > 0);
    … …  
    5455    switch (align) {
    5556      case PS_IMAGE_BINNING_LEFT:
    56         binning->nXoff = 0;
    57         binning->nYoff = 0;
    58         break;
     57        binning->nXoff = 0;
     58        binning->nYoff = 0;
     59        break;
    5960      case PS_IMAGE_BINNING_CENTER:
    60         binning->nXoff = (binning->nXruff * binning->nXbin - binning->nXfine) / 2;
    61         binning->nYoff = (binning->nYruff * binning->nYbin - binning->nYfine) / 2;
    62         break;
     61        binning->nXoff = (binning->nXruff * binning->nXbin - binning->nXfine) / 2;
     62        binning->nYoff = (binning->nYruff * binning->nYbin - binning->nYfine) / 2;
     63        break;
    6364      case PS_IMAGE_BINNING_RIGHT:
    64         binning->nXoff = (binning->nXruff * binning->nXbin - binning->nXfine);
    65         binning->nYoff = (binning->nYruff * binning->nYbin - binning->nYfine);
    66         break;
     65        binning->nXoff = (binning->nXruff * binning->nXbin - binning->nXfine);
     66        binning->nYoff = (binning->nYruff * binning->nYbin - binning->nYfine);
     67        break;
    6768      default:
    68         psAbort ("programming error in %s: impossible case\n", __func__);
     69        psAbort ("programming error in %s: impossible case\n", __func__);
    6970    }
    7071    return;
    … …  
    8182
    8283    if (image != NULL) {
    83         binning->nXskip = image->col0 - binning->nXoff;
    84         binning->nYskip = image->row0 - binning->nYoff;
     84        binning->nXskip = image->col0 - binning->nXoff;
     85        binning->nYskip = image->row0 - binning->nYoff;
    8586    } else {
    86         binning->nXskip = 0 - binning->nXoff;
    87         binning->nYskip = 0 - binning->nYoff;
     87        binning->nXskip = 0 - binning->nXoff;
     88        binning->nYskip = 0 - binning->nYoff;
    8889    }
    8990    return;
    … …  
    9192
    9293void psImageBinningSetScale(psImageBinning *binning, psImageBinningAlign align) {
    93    
     94
    9495    assert (binning->nXfine > 0);
    9596    assert (binning->nYfine > 0);
    … …  
    107108    switch (align) {
    108109      case PS_IMAGE_BINNING_LEFT:
    109         binning->nXoff = 0;
    110         binning->nYoff = 0;
    111         break;
     110        binning->nXoff = 0;
     111        binning->nYoff = 0;
     112        break;
    112113      case PS_IMAGE_BINNING_CENTER:
    113         binning->nXoff = (binning->nXruff * binning->nXbin - binning->nXfine) / 2;
    114         binning->nYoff = (binning->nYruff * binning->nYbin - binning->nYfine) / 2;
    115         break;
     114        binning->nXoff = (binning->nXruff * binning->nXbin - binning->nXfine) / 2;
     115        binning->nYoff = (binning->nYruff * binning->nYbin - binning->nYfine) / 2;
     116        break;
    116117      case PS_IMAGE_BINNING_RIGHT:
    117         binning->nXoff = (binning->nXruff * binning->nXbin - binning->nXfine);
    118         binning->nYoff = (binning->nYruff * binning->nYbin - binning->nYfine);
    119         break;
     118        binning->nXoff = (binning->nXruff * binning->nXbin - binning->nXfine);
     119        binning->nYoff = (binning->nYruff * binning->nYbin - binning->nYfine);
     120        break;
    120121      default:
    121         psAbort ("programming error in %s: impossible case\n", __func__);
     122        psAbort ("programming error in %s: impossible case\n", __func__);
    122123    }
    123124    return;

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

    r14924 r14983  
    77 *  @author Eugene Magnier, IfA
    88 *
    9  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2007-09-20 23:54:25 $
     9 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2007-09-22 03:05:50 $
    1111 *
    1212 *  Copyright 2007 Institute for Astronomy, University of Hawaii
    … …  
    1818
    1919#include <stdio.h>
     20
    2021#include "psError.h"
    2122#include "psAbort.h"
    … …  
    2627#include "psFitsImage.h"
    2728
     29#include "psMemory.h"
    2830#include "psVector.h"
    2931#include "psImage.h"
    … …  
    9698
    9799    // accumulate the values for each map pixel
    98    
     100
    99101    // we can do this by accumulating a vector of pixel indexes for each cell
    100102    psArray *pixelSets = psArrayAlloc (map->map->numCols*map->map->numRows);
    101103    for (int i = 0; i < pixelSets->n; i++) {
    102         pixelSets->data[i] = psVectorAllocEmpty (4, PS_TYPE_S32);
     104        pixelSets->data[i] = psVectorAllocEmpty (4, PS_TYPE_S32);
    103105    }
    104106    // associate each value with a cell
    105107    for (int i = 0; i < x->n; i++) {
    106         int xRuff = psImageBinningGetRuffX (map->binning, x->data.F32[i]);
    107         int yRuff = psImageBinningGetRuffY (map->binning, y->data.F32[i]);
    108 
    109         int bin = xRuff + yRuff*map->map->numCols;
    110         assert (bin >= 0);
    111         assert (bin < pixelSets->n);
    112        
    113         psVector *pixels = pixelSets->data[bin];
    114         pixels->data.S32[pixels->n] = i;
    115         psVectorExtend (pixels, 4, 1);
     108        int xRuff = psImageBinningGetRuffX (map->binning, x->data.F32[i]);
     109        int yRuff = psImageBinningGetRuffY (map->binning, y->data.F32[i]);
     110
     111        int bin = xRuff + yRuff*map->map->numCols;
     112        assert (bin >= 0);
     113        assert (bin < pixelSets->n);
     114
     115        psVector *pixels = pixelSets->data[bin];
     116        pixels->data.S32[pixels->n] = i;
     117        psVectorExtend (pixels, 4, 1);
    116118    }
    117119
    … …  
    123125    int Ny = map->map->numRows;
    124126    for (int iy = 0; iy < Ny; iy++) {
    125         for (int ix = 0; ix < Nx; ix++) {
    126    
    127             // pixel index for this cell
    128             psVector *pixels = pixelSets->data[ix + iy*Nx];
    129 
    130             // storage vectors
    131             psVector *xCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
    132             psVector *yCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
    133             psVector *fCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
    134 
    135             // error vector, if needed
    136             psVector *dfCell = NULL;
    137             if (df) {
    138                 dfCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
    139             }           
    140 
    141             // collect data for this cell
    142             for (int i = 0; i < pixels->n; i++) {
    143                 int bin = pixels->data.S32[i];
    144                 // convert x,y in the fine image to the ruff image
    145                 xCell->data.F32[i]  = psImageBinningGetRuffX (map->binning, x->data.F32[bin]);
    146                 yCell->data.F32[i]  = psImageBinningGetRuffY (map->binning, y->data.F32[bin]);
    147                 fCell->data.F32[i]  = f->data.F32[bin];
    148                 if (df) {
    149                     dfCell->data.F32[i] = df->data.F32[bin];
    150                 }
    151             }
    152 
    153             // reset the stats to avoid contamination from the previous loop
    154             psStatsInit (map->stats);
    155 
    156             // get the value
    157             // XXX need to supply a mask and skip the masked pixels when calculating the centroid
    158             // this will not in general be properly weighted...
    159             if (psVectorStats (map->stats, fCell, dfCell, NULL, 0)) {
    160                 mask->data.U8[iy][ix] = 0;
    161                 // XXX ensure only one option is selected, or save both position and width
    162                 map->map->data.F32[iy][ix] = psStatsGetValue (map->stats, map->stats->options);
    163 
    164                 // calculate the mean position and save:
    165                 psStatsInit (meanStat);
    166                 psVectorStats (meanStat, xCell, NULL, NULL, 0);
    167                 xCoord->data.F32[iy][ix] = psStatsGetValue (meanStat, meanStat->options);
    168                 psStatsInit (meanStat);
    169                 psVectorStats (meanStat, yCell, NULL, NULL, 0);
    170                 yCoord->data.F32[iy][ix] = psStatsGetValue (meanStat, meanStat->options);
    171             } else {
    172                 mask->data.U8[iy][ix] = 1;
    173             }
    174            
    175             psFree (xCell);
    176             psFree (yCell);
    177             psFree (fCell);
    178             psFree (dfCell);
    179         }
     127        for (int ix = 0; ix < Nx; ix++) {
     128
     129            // pixel index for this cell
     130            psVector *pixels = pixelSets->data[ix + iy*Nx];
     131
     132            // storage vectors
     133            psVector *xCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
     134            psVector *yCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
     135            psVector *fCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
     136
     137            // error vector, if needed
     138            psVector *dfCell = NULL;
     139            if (df) {
     140                dfCell = psVectorAlloc (pixels->n, PS_TYPE_F32);
     141            }
     142
     143            // collect data for this cell
     144            for (int i = 0; i < pixels->n; i++) {
     145                int bin = pixels->data.S32[i];
     146                // convert x,y in the fine image to the ruff image
     147                xCell->data.F32[i]  = psImageBinningGetRuffX (map->binning, x->data.F32[bin]);
     148                yCell->data.F32[i]  = psImageBinningGetRuffY (map->binning, y->data.F32[bin]);
     149                fCell->data.F32[i]  = f->data.F32[bin];
     150                if (df) {
     151                    dfCell->data.F32[i] = df->data.F32[bin];
     152                }
     153            }
     154
     155            // reset the stats to avoid contamination from the previous loop
     156            psStatsInit (map->stats);
     157
     158            // get the value
     159            // XXX need to supply a mask and skip the masked pixels when calculating the centroid
     160            // this will not in general be properly weighted...
     161            if (psVectorStats (map->stats, fCell, dfCell, NULL, 0)) {
     162                mask->data.U8[iy][ix] = 0;
     163                // XXX ensure only one option is selected, or save both position and width
     164                map->map->data.F32[iy][ix] = psStatsGetValue (map->stats, map->stats->options);
     165
     166                // calculate the mean position and save:
     167                psStatsInit (meanStat);
     168                psVectorStats (meanStat, xCell, NULL, NULL, 0);
     169                xCoord->data.F32[iy][ix] = psStatsGetValue (meanStat, meanStat->options);
     170                psStatsInit (meanStat);
     171                psVectorStats (meanStat, yCell, NULL, NULL, 0);
     172                yCoord->data.F32[iy][ix] = psStatsGetValue (meanStat, meanStat->options);
     173            } else {
     174                mask->data.U8[iy][ix] = 1;
     175            }
     176
     177            psFree (xCell);
     178            psFree (yCell);
     179            psFree (fCell);
     180            psFree (dfCell);
     181        }
    180182    }
    181183    psFree (pixelSets);
    … …  
    193195    map->nGood = mask->numCols * mask->numRows - map->nBad - map->nPoor;
    194196    if (map->nBad > badFrac * mask->numCols * mask->numRows) {
    195         psFree (xCoord);
    196         psFree (yCoord);
    197         psFree (mask);
    198         return false;
     197        psFree (xCoord);
     198        psFree (yCoord);
     199        psFree (mask);
     200        return false;
    199201    }
    200202
    … …  
    222224// using the points given, generate a map with maximum resolution that yields only good and ok pixels
    223225bool psImageMapGenerateScale (psImageMap *map, psVector *x, psVector *y, psVector *f, psVector *df, float badFrac) {
    224    
     226
    225227    int nXruff, nYruff;
    226        
     228
    227229    while (!psImageMapGenerate (map, x, y, f, df, badFrac)) {
    228         // if we failed to build an acceptable map, decrease nXruff, nYruff as appropriate, and
    229         // try again...  try to keep the aspect ratio.
    230 
    231         // if both axes are at 1, give up
    232         if ((map->binning->nXruff == 1) && (map->binning->nYruff == 1)) {
    233             return false;
    234         }
    235 
    236         // if one axis is at 1, decrement the other
    237         if (map->binning->nXruff == 1) {
    238             nXruff = map->binning->nXruff;
    239             nYruff = map->binning->nYruff - 1;
    240             psImageMapModifyScale (map, nXruff, nYruff);
    241             continue;
    242         }
    243         if (map->binning->nYruff == 1) {
    244             nYruff = map->binning->nYruff;
    245             nXruff = map->binning->nXruff - 1;
    246             psImageMapModifyScale (map, nXruff, nYruff);
    247             continue;
    248         }
    249 
    250         // otherwise, decrement the larger axis, and set the smaller based
    251         // on the aspect ratio
    252         float aRatio = map->binning->nXruff / map->binning->nYruff;
    253         if (map->binning->nXruff > map->binning->nYruff) {
    254             nXruff = map->binning->nXruff - 1;
    255             nYruff = (int)(0.5 + (nXruff / aRatio));
    256         } else {
    257             nYruff = map->binning->nYruff - 1;
    258             nXruff = (int)(0.5 + (nYruff * aRatio));
    259         }
    260 
    261         psImageMapModifyScale (map, nXruff, nYruff);
     230        // if we failed to build an acceptable map, decrease nXruff, nYruff as appropriate, and
     231        // try again...  try to keep the aspect ratio.
     232
     233        // if both axes are at 1, give up
     234        if ((map->binning->nXruff == 1) && (map->binning->nYruff == 1)) {
     235            return false;
     236        }
     237
     238        // if one axis is at 1, decrement the other
     239        if (map->binning->nXruff == 1) {
     240            nXruff = map->binning->nXruff;
     241            nYruff = map->binning->nYruff - 1;
     242            psImageMapModifyScale (map, nXruff, nYruff);
     243            continue;
     244        }
     245        if (map->binning->nYruff == 1) {
     246            nYruff = map->binning->nYruff;
     247            nXruff = map->binning->nXruff - 1;
     248            psImageMapModifyScale (map, nXruff, nYruff);
     249            continue;
     250        }
     251
     252        // otherwise, decrement the larger axis, and set the smaller based
     253        // on the aspect ratio
     254        float aRatio = map->binning->nXruff / map->binning->nYruff;
     255        if (map->binning->nXruff > map->binning->nYruff) {
     256            nXruff = map->binning->nXruff - 1;
     257            nYruff = (int)(0.5 + (nXruff / aRatio));
     258        } else {
     259            nYruff = map->binning->nYruff - 1;
     260            nXruff = (int)(0.5 + (nYruff * aRatio));
     261        }
     262
     263        psImageMapModifyScale (map, nXruff, nYruff);
    262264    }
    263265    return true;
    … …  
    271273    result = psImageUnbinPixel_V2(x, y, map->map, map->binning);
    272274
    273     return result; 
     275    return result;
    274276}
    275277
    … …  
    284286
    285287    for (int i = 0; i < x->n; i++) {
    286         result->data.F32[i] = psImageUnbinPixel_V2(x->data.F32[i], y->data.F32[i], map->map, map->binning);
    287     }
    288 
    289     return result; 
    290 }
     288        result->data.F32[i] = psImageUnbinPixel_V2(x->data.F32[i], y->data.F32[i], map->map, map->binning);
     289    }
     290
     291    return result;
     292}

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

    r14960 r14983  
    77 *  @author Eugene Magnier, IfA
    88 *
    9  *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2007-09-21 02:45:33 $
     9 *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2007-09-22 03:05:50 $
    1111 *
    1212 *  Copyright 2007 Institute for Astronomy, University of Hawaii
    … …  
    1818
    1919#include <stdio.h>
     20#include "psMemory.h"
    2021#include "psError.h"
    2122#include "psAbort.h"
    … …  
    5758    // no spatial information, just calculate mean & stdev
    5859    if ((Nx == 1) && (Ny == 1)) {
    59         psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    60        
    61         // XXX does ROBUST_MEDIAN work with weight?
    62         psVectorStats (stats, f, NULL, mask, maskValue);
    63 
    64         map->map->data.F32[0][0]   = stats->robustMedian;
    65         map->error->data.F32[0][0] = stats->robustStdev;
    66         psFree (stats);
    67         return true;
     60        psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     61
     62        // XXX does ROBUST_MEDIAN work with weight?
     63        psVectorStats (stats, f, NULL, mask, maskValue);
     64
     65        map->map->data.F32[0][0]   = stats->robustMedian;
     66        map->error->data.F32[0][0] = stats->robustStdev;
     67        psFree (stats);
     68        return true;
    6869    }
    6970
    7071    if (Nx == 1) {
    71         psAbort ("un-implemented edge case");
    72         goto insert;
     72        psAbort ("un-implemented edge case");
     73        goto insert;
    7374    }
    7475    if (Ny == 1) {
    75         psAbort ("un-implemented edge case");
    76         goto insert;
     76        psAbort ("un-implemented edge case");
     77        goto insert;
    7778    }
    7879
    … …  
    8283
    8384    for (int i = 0; i < 3; i++) {
    84         SAv[i] = SAm[i] + 1;
    85         // TAv[i] = TAm[i] + 1;
     85        SAv[i] = SAm[i] + 1;
     86        // TAv[i] = TAm[i] + 1;
    8687    }
    8788    sA = SAv + 1;
    8889    // tA = TAv + 1;
    8990
    90     // elements of the matrix equation Ax = B; we are solving for the vector x 
     91    // elements of the matrix equation Ax = B; we are solving for the vector x
    9192    psImage *A = psImageAlloc (Nx*Ny, Nx*Ny, PS_TYPE_F32);
    9293    psVector *B = psVectorAlloc (Nx*Ny, PS_TYPE_F32);
    … …  
    9596    psVectorInit (B, 0.0);
    9697
    97     // we are looping over the Nx,Ny image map elements; 
    98     // the matrix equation contains Nx*Ny rows and columns 
     98    // we are looping over the Nx,Ny image map elements;
     99    // the matrix equation contains Nx*Ny rows and columns
    99100    // for (int n = 1; n < Nx - 1; n++) {
    100101    // for (int m = 1; m < Ny - 1; m++) {
    101    
     102
    102103    // float Total = 0.0;
    103104    for (int n = 0; n < Nx; n++) {
    104         for (int m = 0; m < Ny; m++) {
    105             // define & init summing variables
    106             float rx_rx_ry_ry = 0;
    107             float rx_rx_dy_ry = 0;
    108             float dx_rx_ry_ry = 0;
    109             float dx_rx_dy_ry = 0;
    110             float fi_rx_ry    = 0;
    111             float rx_rx_py_py = 0;
    112             float rx_rx_qy_py = 0;
    113             float dx_rx_py_py = 0;
    114             float dx_rx_qy_py = 0;
    115             float fi_rx_py    = 0;
    116             float px_px_ry_ry = 0;
    117             float px_px_dy_ry = 0;
    118             float qx_px_ry_ry = 0;
    119             float qx_px_dy_ry = 0;
    120             float fi_px_ry    = 0;
    121             float px_px_py_py = 0;
    122             float px_px_qy_py = 0;
    123             float qx_px_py_py = 0;
    124             float qx_px_qy_py = 0;
    125             float fi_px_py    = 0;
    126 
    127             // generate the sums for the fitting matrix element I,J
    128             // I = n + nX*m
    129             // J = (n + jn) + nX*(m + jm)
    130             for (int i = 0; i < x->n; i++) {
    131 
    132                 if (mask && (mask->data.U8[i] & maskValue)) continue;
    133 
    134                 // base coordinate offset for this point (x,y) relative to this map element (n,m)
    135                 // float dx = x->data.F32[i] - psImageBinningGetFineX (map->binning, n + 0.5);
    136                 // float dy = y->data.F32[i] - psImageBinningGetFineY (map->binning, m + 0.5);
    137 
    138                 float dx = psImageBinningGetRuffX (map->binning, x->data.F32[i]) - (n + 0.5);
    139                 float dy = psImageBinningGetRuffY (map->binning, y->data.F32[i]) - (m + 0.5);
    140 
    141                 // edge cases to include:
    142                 bool edgeX = false;
    143                 edgeX |= ((n == 1) && (dx < -1.0));
    144                 edgeX |= ((n == Nx - 2) && (dx > +1.0));
    145 
    146                 bool edgeY = false;
    147                 edgeY |= ((m == 1) && (dy < -1.0));
    148                 edgeY |= ((m == Ny - 2) && (dy > +1.0));
    149                
    150                 // skip points outside of 2x2 grid centered on n,m:
    151                 if (!edgeX && (fabs(dx) > 1.0)) continue;
    152                 if (!edgeY && (fabs(dy) > 1.0)) continue;
    153 
    154                 // related offset values
    155                 float rx = 1.0 - dx;
    156                 float ry = 1.0 - dy;
    157                 float px = 1.0 + dx;
    158                 float py = 1.0 + dy;
    159                 float qx = -dx;
    160                 float qy = -dy;
    161 
    162                 // data value & weight for this point
    163                 float fi = f->data.F32[i];
    164                 float wt = 1.0;
    165                 if (df != NULL) {
    166                     if (df->data.F32[i] == 0.0) {
    167                         wt = 0.0;
    168                     } else {
    169                         wt = 1.0 / PS_SQR(df->data.F32[i]); // XXX test for dz == NULL or dz_i = 0
    170                     }
    171                 }
    172 
    173                 // sum the appropriate elements for the different quadrants
    174 
    175                 int Qx = (dx >= 0) ? 1 : 0;
    176                 if (n ==      0) Qx = 1;
    177                 if (n == Nx - 1) Qx = 0;
    178 
    179                 int Qy = (dy >= 0) ? 1 : 0;
    180                 if (m ==      0) Qy = 1;
    181                 if (m == Ny - 1) Qy = 0;
    182 
    183                 // points at offset 1,1
    184                 if ((Qx == 1) && (Qy == 1)) {
    185                     rx_rx_ry_ry += rx*rx*ry*ry*wt;
    186                     rx_rx_dy_ry += rx*rx*dy*ry*wt;
    187                     dx_rx_ry_ry += dx*rx*ry*ry*wt;
    188                     dx_rx_dy_ry += dx*rx*dy*ry*wt;
    189                     fi_rx_ry    += fi*rx*ry*wt;
    190                 }
    191                 // points at offset 1,0
    192                 if ((Qx == 1) && (Qy == 0)) {
    193                     rx_rx_py_py += rx*rx*py*py*wt;
    194                     rx_rx_qy_py += rx*rx*qy*py*wt;
    195                     dx_rx_py_py += dx*rx*py*py*wt;
    196                     dx_rx_qy_py += dx*rx*qy*py*wt;
    197                     fi_rx_py    += fi*rx*py*wt;
    198                 }
    199                 // points at offset 0,1
    200                 if ((Qx == 0) && (Qy == 1)) {
    201                     px_px_ry_ry += px*px*ry*ry*wt;
    202                     px_px_dy_ry += px*px*dy*ry*wt;
    203                     qx_px_ry_ry += qx*px*ry*ry*wt;
    204                     qx_px_dy_ry += qx*px*dy*ry*wt;
    205                     fi_px_ry    += fi*px*ry*wt;
    206                 }
    207                 // points at offset 0,0
    208                 if ((Qx == 0) && (Qy == 0)) {
    209                     px_px_py_py += px*px*py*py*wt;
    210                     px_px_qy_py += px*px*qy*py*wt;
    211                     qx_px_py_py += qx*px*py*py*wt;
    212                     qx_px_qy_py += qx*px*qy*py*wt;
    213                     fi_px_py    += fi*px*py*wt;
    214                 }
    215             }           
    216 
    217             // the chi-square derivatives have elements of the form g(n+jn,m+jm)*A(jn,jm),
    218             // jn,jm = -1 to +1. Convert the sums above into the correct coefficients
    219             sA[-1][-1] = qx_px_qy_py;
    220             sA[-1][ 0] = qx_px_ry_ry + qx_px_py_py;
    221             sA[-1][+1] = qx_px_dy_ry;
    222             sA[ 0][-1] = rx_rx_qy_py + px_px_qy_py;
    223             sA[ 0][ 0] = rx_rx_ry_ry + px_px_ry_ry + rx_rx_py_py + px_px_py_py;
    224             sA[ 0][+1] = rx_rx_dy_ry + px_px_dy_ry;
    225             sA[+1][-1] = dx_rx_qy_py;
    226             sA[+1][ 0] = dx_rx_ry_ry + dx_rx_py_py;
    227             sA[+1][+1] = dx_rx_dy_ry;
    228 
    229         insert:
    230             // I[ 0][ 0] = index for this n,m element:
    231             I = n + Nx * m;
    232             B->data.F32[I] = fi_rx_ry + fi_rx_py + fi_px_ry + fi_px_py;
    233 
    234             // insert these values into their corresponding locations in A, B
    235             // float Sum = 0.0;
    236             for (int jn = -1; jn <= +1; jn++) {
    237                 if (n + jn <   0) continue;
    238                 if (n + jn >= Nx) continue;
    239                 for (int jm = -1; jm <= +1; jm++) {
    240                     if (m + jm <   0) continue;
    241                     if (m + jm >= Ny) continue;
    242                     J = (n + jn) + Nx * (m + jm);
    243                     A->data.F32[J][I] = sA[jn][jm];
    244                     // fprintf (stderr, "A %d %d (%d %d : %d %d): %f\n", I, J, n, m, n + jn, m + jm, sA[jn][jm]);
    245                     // Sum += sA[jn][jm];
    246                 }
    247             }
    248             // fprintf (stderr, "B %d (%d %d) : %f  :  %f\n", I, n, m, B->data.F32[I], Sum);
    249             // Total += Sum;
    250         }
     105        for (int m = 0; m < Ny; m++) {
     106            // define & init summing variables
     107            float rx_rx_ry_ry = 0;
     108            float rx_rx_dy_ry = 0;
     109            float dx_rx_ry_ry = 0;
     110            float dx_rx_dy_ry = 0;
     111            float fi_rx_ry    = 0;
     112            float rx_rx_py_py = 0;
     113            float rx_rx_qy_py = 0;
     114            float dx_rx_py_py = 0;
     115            float dx_rx_qy_py = 0;
     116            float fi_rx_py    = 0;
     117            float px_px_ry_ry = 0;
     118            float px_px_dy_ry = 0;
     119            float qx_px_ry_ry = 0;
     120            float qx_px_dy_ry = 0;
     121            float fi_px_ry    = 0;
     122            float px_px_py_py = 0;
     123            float px_px_qy_py = 0;
     124            float qx_px_py_py = 0;
     125            float qx_px_qy_py = 0;
     126            float fi_px_py    = 0;
     127
     128            // generate the sums for the fitting matrix element I,J
     129            // I = n + nX*m
     130            // J = (n + jn) + nX*(m + jm)
     131            for (int i = 0; i < x->n; i++) {
     132
     133                if (mask && (mask->data.U8[i] & maskValue)) continue;
     134
     135                // base coordinate offset for this point (x,y) relative to this map element (n,m)
     136                // float dx = x->data.F32[i] - psImageBinningGetFineX (map->binning, n + 0.5);
     137                // float dy = y->data.F32[i] - psImageBinningGetFineY (map->binning, m + 0.5);
     138
     139                float dx = psImageBinningGetRuffX (map->binning, x->data.F32[i]) - (n + 0.5);
     140                float dy = psImageBinningGetRuffY (map->binning, y->data.F32[i]) - (m + 0.5);
     141
     142                // edge cases to include:
     143                bool edgeX = false;
     144                edgeX |= ((n == 1) && (dx < -1.0));
     145                edgeX |= ((n == Nx - 2) && (dx > +1.0));
     146
     147                bool edgeY = false;
     148                edgeY |= ((m == 1) && (dy < -1.0));
     149                edgeY |= ((m == Ny - 2) && (dy > +1.0));
     150
     151                // skip points outside of 2x2 grid centered on n,m:
     152                if (!edgeX && (fabs(dx) > 1.0)) continue;
     153                if (!edgeY && (fabs(dy) > 1.0)) continue;
     154
     155                // related offset values
     156                float rx = 1.0 - dx;
     157                float ry = 1.0 - dy;
     158                float px = 1.0 + dx;
     159                float py = 1.0 + dy;
     160                float qx = -dx;
     161                float qy = -dy;
     162
     163                // data value & weight for this point
     164                float fi = f->data.F32[i];
     165                float wt = 1.0;
     166                if (df != NULL) {
     167                    if (df->data.F32[i] == 0.0) {
     168                        wt = 0.0;
     169                    } else {
     170                        wt = 1.0 / PS_SQR(df->data.F32[i]); // XXX test for dz == NULL or dz_i = 0
     171                    }
     172                }
     173
     174                // sum the appropriate elements for the different quadrants
     175
     176                int Qx = (dx >= 0) ? 1 : 0;
     177                if (n ==      0) Qx = 1;
     178                if (n == Nx - 1) Qx = 0;
     179
     180                int Qy = (dy >= 0) ? 1 : 0;
     181                if (m ==      0) Qy = 1;
     182                if (m == Ny - 1) Qy = 0;
     183
     184                // points at offset 1,1
     185                if ((Qx == 1) && (Qy == 1)) {
     186                    rx_rx_ry_ry += rx*rx*ry*ry*wt;
     187                    rx_rx_dy_ry += rx*rx*dy*ry*wt;
     188                    dx_rx_ry_ry += dx*rx*ry*ry*wt;
     189                    dx_rx_dy_ry += dx*rx*dy*ry*wt;
     190                    fi_rx_ry    += fi*rx*ry*wt;
     191                }
     192                // points at offset 1,0
     193                if ((Qx == 1) && (Qy == 0)) {
     194                    rx_rx_py_py += rx*rx*py*py*wt;
     195                    rx_rx_qy_py += rx*rx*qy*py*wt;
     196                    dx_rx_py_py += dx*rx*py*py*wt;
     197                    dx_rx_qy_py += dx*rx*qy*py*wt;
     198                    fi_rx_py    += fi*rx*py*wt;
     199                }
     200                // points at offset 0,1
     201                if ((Qx == 0) && (Qy == 1)) {
     202                    px_px_ry_ry += px*px*ry*ry*wt;
     203                    px_px_dy_ry += px*px*dy*ry*wt;
     204                    qx_px_ry_ry += qx*px*ry*ry*wt;
     205                    qx_px_dy_ry += qx*px*dy*ry*wt;
     206                    fi_px_ry    += fi*px*ry*wt;
     207                }
     208                // points at offset 0,0
     209                if ((Qx == 0) && (Qy == 0)) {
     210                    px_px_py_py += px*px*py*py*wt;
     211                    px_px_qy_py += px*px*qy*py*wt;
     212                    qx_px_py_py += qx*px*py*py*wt;
     213                    qx_px_qy_py += qx*px*qy*py*wt;
     214                    fi_px_py    += fi*px*py*wt;
     215                }
     216            }
     217
     218            // the chi-square derivatives have elements of the form g(n+jn,m+jm)*A(jn,jm),
     219            // jn,jm = -1 to +1. Convert the sums above into the correct coefficients
     220            sA[-1][-1] = qx_px_qy_py;
     221            sA[-1][ 0] = qx_px_ry_ry + qx_px_py_py;
     222            sA[-1][+1] = qx_px_dy_ry;
     223            sA[ 0][-1] = rx_rx_qy_py + px_px_qy_py;
     224            sA[ 0][ 0] = rx_rx_ry_ry + px_px_ry_ry + rx_rx_py_py + px_px_py_py;
     225            sA[ 0][+1] = rx_rx_dy_ry + px_px_dy_ry;
     226            sA[+1][-1] = dx_rx_qy_py;
     227            sA[+1][ 0] = dx_rx_ry_ry + dx_rx_py_py;
     228            sA[+1][+1] = dx_rx_dy_ry;
     229
     230        insert:
     231            // I[ 0][ 0] = index for this n,m element:
     232            I = n + Nx * m;
     233            B->data.F32[I] = fi_rx_ry + fi_rx_py + fi_px_ry + fi_px_py;
     234
     235            // insert these values into their corresponding locations in A, B
     236            // float Sum = 0.0;
     237            for (int jn = -1; jn <= +1; jn++) {
     238                if (n + jn <   0) continue;
     239                if (n + jn >= Nx) continue;
     240                for (int jm = -1; jm <= +1; jm++) {
     241                    if (m + jm <   0) continue;
     242                    if (m + jm >= Ny) continue;
     243                    J = (n + jn) + Nx * (m + jm);
     244                    A->data.F32[J][I] = sA[jn][jm];
     245                    // fprintf (stderr, "A %d %d (%d %d : %d %d): %f\n", I, J, n, m, n + jn, m + jm, sA[jn][jm]);
     246                    // Sum += sA[jn][jm];
     247                }
     248            }
     249            // fprintf (stderr, "B %d (%d %d) : %f  :  %f\n", I, n, m, B->data.F32[I], Sum);
     250            // Total += Sum;
     251        }
    251252    }
    252253    // fprintf (stderr, "Total: %f\n", Total);
    … …  
    256257    psVectorInit (Empty, 0);
    257258    for (int i = 0; i < Nx*Ny; i++) {
    258         if (A->data.F32[i][i] == 0.0) {
    259             Empty->data.S8[i] = 1;
    260             for (int j = 0; j < Nx*Ny; j++) {
    261                 A->data.F32[i][j] = 0.0;
    262                 A->data.F32[j][i] = 0.0;
    263             }
    264             A->data.F32[i][i] = 1.0;
    265             B->data.F32[i] = 0.0;
    266         }
     259        if (A->data.F32[i][i] == 0.0) {
     260            Empty->data.S8[i] = 1;
     261            for (int j = 0; j < Nx*Ny; j++) {
     262                A->data.F32[i][j] = 0.0;
     263                A->data.F32[j][i] = 0.0;
     264            }
     265            A->data.F32[i][i] = 1.0;
     266            B->data.F32[i] = 0.0;
     267        }
    267268    }
    268269
    … …  
    274275    psImage *vector = psImageAlloc (1, B->n, PS_TYPE_F32);
    275276    for (int n = 0; n < B->n; n++) {
    276         vector->data.F32[0][n] = B->data.F32[n];
     277        vector->data.F32[0][n] = B->data.F32[n];
    277278    }
    278279
    … …  
    284285
    285286    if (!psMatrixGJSolveF32(A, B)) {
    286         psAbort ("failed on linear equations");
    287         psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
    288         psFree (A);
    289         psFree (B);
    290         return false;
    291     }
    292    
     287        psAbort ("failed on linear equations");
     288        psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
     289        psFree (A);
     290        psFree (B);
     291        return false;
     292    }
     293
    293294    // set bad values to NaN
    294295    for (int i = 0; i < Nx*Ny; i++) {
    295         if (Empty->data.S8[i]) {
    296             B->data.F32[i] = NAN;
    297         }
     296        if (Empty->data.S8[i]) {
     297            B->data.F32[i] = NAN;
     298        }
    298299    }
    299300
    300301
    301302    for (int n = 0; n < Nx; n++) {
    302         for (int m = 0; m < Ny; m++) {
    303             I = n + Nx * m;
    304             map->map->data.F32[m][n] = B->data.F32[I];
    305             map->error->data.F32[m][n] = sqrt(A->data.F32[I][I]);
    306         }
     303        for (int m = 0; m < Ny; m++) {
     304            I = n + Nx * m;
     305            map->map->data.F32[m][n] = B->data.F32[I];
     306            map->error->data.F32[m][n] = sqrt(A->data.F32[I][I]);
     307        }
    307308    }
    308309
    … …  
    353354    psVector *mask = inMask;
    354355    if (!inMask) {
    355         mask = psVectorAlloc (x->n, PS_TYPE_U8);
    356         psVectorInit (mask, 0);
     356        mask = psVectorAlloc (x->n, PS_TYPE_U8);
     357        psVectorInit (mask, 0);
    357358    }
    358359
    … …  
    368369        if (!psImageMapFit(map, mask, maskValue, x, y, f, df)) {
    369370            psError(PS_ERR_UNKNOWN, false, "Could not fit image map.\n");
    370             psFree(resid);
    371             if (!inMask) psFree (mask);
     371            psFree(resid);
     372            if (!inMask) psFree (mask);
    372373            return false;
    373374        }
    … …  
    377378            psError(PS_ERR_UNKNOWN, false, "Failure in psImageMapEvalVector().\n");
    378379            psFree(resid);
    379             if (!inMask) psFree (mask);
     380            if (!inMask) psFree (mask);
    380381            return false;
    381382        }
    382383        for (int i = 0 ; i < f->n ; i++) {
    383             resid->data.F32[i] = (f->data.F32[i] - fit->data.F32[i]);
     384            resid->data.F32[i] = (f->data.F32[i] - fit->data.F32[i]);
    384385        }
    385386
    … …  
    388389            psFree(resid);
    389390            psFree(fit);
    390             if (!inMask) psFree (mask);
     391            if (!inMask) psFree (mask);
    391392            return false;
    392393        }
    … …  
    404405        // recovery is not allowed with this scheme
    405406        for (psS32 i = 0; i < resid->n; i++) {
    406             // XXX this prevents recovery of previously masked values
     407            // XXX this prevents recovery of previously masked values
    407408            if (mask->data.U8[i] & maskValue) {
    408409                continue;
    … …  
    410411
    411412            if ((resid->data.F32[i] - meanValue > maxClipValue) || (resid->data.F32[i] - meanValue < minClipValue)) {
    412                 psTrace("psLib.imageops", 6, "Masking element %d  : %f vs %f : resid is %f\n", i, f->data.F32[i], fit->data.F32[i], resid->data.F32[i]);
    413                 mask->data.U8[i] |= 0x01;
     413                psTrace("psLib.imageops", 6, "Masking element %d  : %f vs %f : resid is %f\n", i, f->data.F32[i], fit->data.F32[i], resid->data.F32[i]);
     414                mask->data.U8[i] |= 0x01;
    414415                continue;
    415416            }

  • trunk/psLib/src/types/psMetadataConfig.c

    r14686 r14983  
    1111*  @author Joshua Hoblitt, University of Hawaii 2006-2007
    1212*
    13 *  @version $Revision: 1.139 $ $Name: not supported by cvs2svn $
    14 *  @date $Date: 2007-08-29 00:27:13 $
     13*  @version $Revision: 1.140 $ $Name: not supported by cvs2svn $
     14*  @date $Date: 2007-09-22 03:05:50 $
    1515*
    1616*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    3131#include <unistd.h>
    3232
     33#include "psMemory.h"
    3334#include "psMetadataConfig.h"
    3435#include "psAssert.h"
Note: See TracChangeset for help on using the changeset viewer.