IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 14862


Ignore:
Timestamp:
Sep 16, 2007, 3:14:52 PM (19 years ago)
Author:
magnier
Message:

fixed algorithim for psImageMapFit

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branch_20070830/psLib/src/imageops/psImageMapFit.c

    r14859 r14862  
    77 *  @author Eugene Magnier, IfA
    88 *
    9  *  @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2007-09-15 19:43:43 $
     9 *  @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2007-09-17 01:14:52 $
    1111 *
    1212 *  Copyright 2007 Institute for Astronomy, University of Hawaii
    … …  
    2222
    2323// XXX for testing
    24 #include "psFits.h"
    25 #include "psFitsImage.h"
     24// #include "psFits.h"
     25// #include "psFitsImage.h"
    2626
    2727#include "psAssert.h"
    … …  
    5454    // set up the redirection table so we can use sA[-1][-1], etc
    5555    float SAm[3][3], *SAv[3], **sA;
    56     float TAm[3][3], *TAv[3], **tA;
     56    // float TAm[3][3], *TAv[3], **tA;
    5757
    5858    for (int i = 0; i < 3; i++) {
    5959        SAv[i] = SAm[i] + 1;
    60         TAv[i] = TAm[i] + 1;
     60        // TAv[i] = TAm[i] + 1;
    6161    }
    6262    sA = SAv + 1;
    63     tA = TAv + 1;
     63    // tA = TAv + 1;
    6464
    6565    // elements of the matrix equation Ax = B; we are solving for the vector x
    6666    psImage *A = psImageAlloc (Nx*Ny, Nx*Ny, PS_TYPE_F32);
    6767    psVector *B = psVectorAlloc (Nx*Ny, PS_TYPE_F32);
     68
     69    psImageInit (A, 0.0);
     70    psVectorInit (B, 0.0);
    6871
    6972    // we are looping over the Nx,Ny image map elements;
    … …  
    7275    // for (int m = 1; m < Ny - 1; m++) {
    7376   
    74     float Total = 0.0;
     77    // float Total = 0.0;
    7578    for (int n = 0; n < Nx; n++) {
    7679        for (int m = 0; m < Ny; m++) {
    … …  
    108111                float dy = psImageBinningGetRuffY (map->binning, y->data.F32[i]) - (m + 0.5);
    109112
     113                // edge cases to include:
     114                bool edgeX = false;
     115                edgeX |= ((n == 1) && (dx < -1.0));
     116                edgeX |= ((n == Nx - 2) && (dx > +1.0));
     117
     118                bool edgeY = false;
     119                edgeY |= ((m == 1) && (dy < -1.0));
     120                edgeY |= ((m == Ny - 2) && (dy > +1.0));
     121               
    110122                // skip points outside of 2x2 grid centered on n,m:
    111                 if (fabs(dx) > 1.0) continue;
    112                 if (fabs(dy) > 1.0) continue;
     123                if (!edgeX && (fabs(dx) > 1.0)) continue;
     124                if (!edgeY && (fabs(dy) > 1.0)) continue;
    113125
    114126                // related offset values
    … …  
    133145                // sum the appropriate elements for the different quadrants
    134146
     147                int Qx = (dx >= 0) ? 1 : 0;
     148                if (n ==      0) Qx = 1;
     149                if (n == Nx - 1) Qx = 0;
     150
     151                int Qy = (dy >= 0) ? 1 : 0;
     152                if (m ==      0) Qy = 1;
     153                if (m == Ny - 1) Qy = 0;
     154
    135155                // points at offset 1,1
    136                 if ((dx >= 0) && (dy >= 0)) {
     156                if ((Qx == 1) && (Qy == 1)) {
    137157                    rx_rx_ry_ry += rx*rx*ry*ry*wt;
    138158                    rx_rx_dy_ry += rx*rx*dy*ry*wt;
    … …  
    142162                }
    143163                // points at offset 1,0
    144                 if ((dx >= 0) && (dy < 0)) {
     164                if ((Qx == 1) && (Qy == 0)) {
    145165                    rx_rx_py_py += rx*rx*py*py*wt;
    146166                    rx_rx_qy_py += rx*rx*qy*py*wt;
    … …  
    150170                }
    151171                // points at offset 0,1
    152                 if ((dx <  0) && (dy >= 0)) {
     172                if ((Qx == 0) && (Qy == 1)) {
    153173                    px_px_ry_ry += px*px*ry*ry*wt;
    154174                    px_px_dy_ry += px*px*dy*ry*wt;
    … …  
    158178                }
    159179                // points at offset 0,0
    160                 if ((dx <  0) && (dy < 0)) {
     180                if ((Qx == 0) && (Qy == 0)) {
    161181                    px_px_py_py += px*px*py*py*wt;
    162182                    px_px_qy_py += px*px*qy*py*wt;
    … …  
    178198            sA[+1][ 0] = dx_rx_ry_ry + dx_rx_py_py;
    179199            sA[+1][+1] = dx_rx_dy_ry;
    180            
    181             // in for edge cases, we only have certain grid points available.  Depending on location
    182             // in the plane, we need to re-arrange the terms to account for this.  we
    183             // implicitly extrapolate across the g[n][m] grid point, assuming, eg g[n+1][m] =
    184             // 2*g[n][m] - g[n-1][m].  generate the true output coefficients for the different
    185             // cases.  start by zeroing out this array.
    186             for (int jn = -1; jn <= +1; jn++) {
    187                 for (int jm = -1; jm <= +1; jm++) {
    188                     tA[jn][jm] = 0;
    189                 }
    190             }
    191 
    192             // XXX test:
    193             // goto skip_ahead;
    194200
    195201            // XXX need to add three additional edge cases:
    … …  
    207213            }
    208214
    209             // UPPER-RIGHT Corner
    210             if ((n == Nx - 1) && (m == Ny - 1)) {
    211                 tA[-1][-1] = sA[-1][-1] -   sA[+1][+1] -   sA[+1][-1] -   sA[-1][+1];
    212                 tA[ 0][-1] = sA[ 0][-1] + 2*sA[+1][-1] -   sA[ 0][+1];
    213                 tA[-1][ 0] = sA[-1][ 0] + 2*sA[-1][+1] -   sA[+1][ 0];
    214                 tA[ 0][ 0] = sA[ 0][ 0] + 2*sA[+1][ 0] + 2*sA[ 0][+1] + 2*sA[+1][+1];
    215                 goto insert;
    216             }
    217             // LOWER-RIGHT Corner
    218             if ((n == Nx - 1) && (m == 0)) {
    219                 tA[-1][+1] = sA[-1][+1] -   sA[+1][-1] -   sA[+1][+1] -   sA[-1][-1];
    220                 tA[ 0][+1] = sA[ 0][+1] + 2*sA[+1][+1] -   sA[ 0][-1];
    221                 tA[-1][ 0] = sA[-1][ 0] + 2*sA[-1][-1] -   sA[+1][ 0];
    222                 tA[ 0][ 0] = sA[ 0][ 0] + 2*sA[+1][ 0] + 2*sA[ 0][-1] + 2*sA[+1][-1];
    223                 goto insert;
    224             }
    225             // LOWER-LEFT Corner
    226             if ((n == 0) && (m == 0)) {
    227                 tA[+1][+1] = sA[+1][+1] -   sA[-1][-1] -   sA[-1][+1] -   sA[+1][-1];
    228                 tA[ 0][+1] = sA[ 0][+1] + 2*sA[-1][+1] -   sA[ 0][-1];
    229                 tA[+1][ 0] = sA[+1][ 0] + 2*sA[+1][-1] -   sA[-1][ 0];
    230                 tA[ 0][ 0] = sA[ 0][ 0] + 2*sA[-1][ 0] + 2*sA[ 0][-1] + 2*sA[-1][-1];
    231                 goto insert;
    232             }
    233             // UPPER-LEFT Corner
    234             if ((n == 0) && (m == Ny - 1)) {
    235                 tA[+1][-1] = sA[+1][-1] -   sA[-1][+1] -   sA[-1][-1] -   sA[+1][+1];
    236                 tA[ 0][-1] = sA[ 0][-1] + 2*sA[-1][-1] -   sA[ 0][+1];
    237                 tA[+1][ 0] = sA[+1][ 0] + 2*sA[+1][+1] -   sA[-1][ 0];
    238                 tA[ 0][ 0] = sA[ 0][ 0] + 2*sA[-1][ 0] + 2*sA[ 0][+1] + 2*sA[-1][+1];
    239                 goto insert;
    240             }
    241             // LEFT Edge
    242             if (n == 0) {
    243                 tA[+1][-1] = sA[+1][-1] - sA[-1][-1];
    244                 tA[+1][ 0] = sA[+1][ 0] - sA[-1][ 0];
    245                 tA[+1][+1] = sA[+1][+1] - sA[-1][+1];
    246                 tA[ 0][-1] = sA[ 0][-1] + 2.0*sA[-1][-1];
    247                 tA[ 0][ 0] = sA[ 0][ 0] + 2.0*sA[-1][ 0];
    248                 tA[ 0][+1] = sA[ 0][+1] + 2.0*sA[-1][+1];
    249                 goto insert;
    250             }
    251             // RIGHT Edge
    252             if (n == Nx - 1) {
    253                 tA[-1][-1] = sA[-1][-1] - sA[+1][-1];
    254                 tA[-1][ 0] = sA[-1][ 0] - sA[+1][ 0];
    255                 tA[-1][+1] = sA[-1][+1] - sA[+1][+1];
    256                 tA[ 0][-1] = sA[ 0][-1] + 2.0*sA[+1][-1];
    257                 tA[ 0][ 0] = sA[ 0][ 0] + 2.0*sA[+1][ 0];
    258                 tA[ 0][+1] = sA[ 0][+1] + 2.0*sA[+1][+1];
    259                 goto insert;
    260             }
    261             // BOTTOM Edge
    262             if (m == 0) {
    263                 tA[-1][+1] = sA[-1][+1] - sA[-1][-1];
    264                 tA[ 0][+1] = sA[ 0][+1] - sA[ 0][-1];
    265                 tA[+1][+1] = sA[+1][+1] - sA[+1][-1];
    266                 tA[-1][ 0] = sA[-1][ 0] + 2.0*sA[-1][-1];
    267                 tA[ 0][ 0] = sA[ 0][ 0] + 2.0*sA[ 0][-1];
    268                 tA[+1][ 0] = sA[+1][ 0] + 2.0*sA[+1][-1];
    269                 goto insert;
    270             }
    271             // TOP Edge
    272             if (m == Ny - 1) {
    273                 tA[-1][-1] = sA[-1][-1] - sA[-1][+1];
    274                 tA[ 0][-1] = sA[ 0][-1] - sA[ 0][+1];
    275                 tA[+1][-1] = sA[+1][-1] - sA[+1][+1];
    276                 tA[-1][ 0] = sA[-1][ 0] + 2.0*sA[-1][+1];
    277                 tA[ 0][ 0] = sA[ 0][ 0] + 2.0*sA[ 0][+1];
    278                 tA[+1][ 0] = sA[+1][ 0] + 2.0*sA[+1][+1];
    279                 goto insert;
    280             }
    281 
    282             // skip_ahead:
    283             // center pixels:
    284             for (int jn = -1; jn <= +1; jn++) {
    285                 for (int jm = -1; jm <= +1; jm++) {
    286                     tA[jn][jm] = sA[jn][jm];
    287                 }
    288             }
    289 
    290215        insert:
    291216            // I[ 0][ 0] = index for this n,m element:
    … …  
    294219
    295220            // insert these values into their corresponding locations in A, B
    296             float Sum = 0.0;
     221            // float Sum = 0.0;
    297222            for (int jn = -1; jn <= +1; jn++) {
    298223                if (n + jn <   0) continue;
    … …  
    302227                    if (m + jm >= Ny) continue;
    303228                    J = (n + jn) + Nx * (m + jm);
    304                     A->data.F32[J][I] = tA[jn][jm];
    305                     fprintf (stderr, "A %d %d (%d %d : %d %d): %f\n", I, J, n, m, n + jn, m + jm, tA[jn][jm]);
    306                     Sum += tA[jn][jm];
    307                 }
    308             }
    309             fprintf (stderr, "B %d (%d %d) : %f  :  %f\n", I, n, m, B->data.F32[I], Sum);
    310             Total += Sum;
    311         }
    312     }
    313 
    314     fprintf (stderr, "Total: %f\n", Total);
     229                    A->data.F32[J][I] = sA[jn][jm];
     230                    // fprintf (stderr, "A %d %d (%d %d : %d %d): %f\n", I, J, n, m, n + jn, m + jm, sA[jn][jm]);
     231                    // Sum += sA[jn][jm];
     232                }
     233            }
     234            // fprintf (stderr, "B %d (%d %d) : %f  :  %f\n", I, n, m, B->data.F32[I], Sum);
     235            // Total += Sum;
     236        }
     237    }
     238    // fprintf (stderr, "Total: %f\n", Total);
    315239
    316240    // write 1s in unused diagonal elements
    … …  
    331255    # endif
    332256
    333     psFits *fits = psFitsOpen ("Agj.fits", "w");
    334     psFitsWriteImage (fits, NULL, A, 0, NULL);
    335     psFitsClose (fits);
     257    // psFits *fits = psFitsOpen ("Agj.fits", "w");
     258    // psFitsWriteImage (fits, NULL, A, 0, NULL);
     259    // psFitsClose (fits);
    336260
    337261    if (!psMatrixGJSolveF32(A, B)) {
Note: See TracChangeset for help on using the changeset viewer.