IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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

Including psMemory.h where required.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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            }
Note: See TracChangeset for help on using the changeset viewer.