IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 3540


Ignore:
Timestamp:
Mar 29, 2005, 9:41:56 AM (21 years ago)
Author:
gusciora
Message:

Plane transform invert, combine, and fit.

Location:
trunk/psLib/src
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/astro/psCoord.c

    r3470 r3540  
    1010*  @author GLG, MHPCC
    1111*
    12 *  @version $Revision: 1.59 $ $Name: not supported by cvs2svn $
    13 *  @date $Date: 2005-03-22 01:35:49 $
     12*  @version $Revision: 1.60 $ $Name: not supported by cvs2svn $
     13*  @date $Date: 2005-03-29 19:41:56 $
    1414*
    1515*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    191191    return (out);
    192192}
    193 
    194 // XXX: Must code this.
    195 psPlaneTransform *psPlaneTransformInvert(
    196     psPlaneTransform *out,
    197     const psPlaneTransform *in,
    198     psRegion *region,
    199     int nSamples)
    200 {
    201     PS_PTR_CHECK_NULL(in, NULL);
    202     PS_PTR_CHECK_NULL(region, NULL);
    203 
    204     return(NULL);
    205 }
    206 
    207 
    208 // XXX: Must code this.
    209 psPlaneTransform *psPlaneTransformCombine(
    210     psPlaneTransform *out,
    211     const psPlaneTransform *trans1,
    212     const psPlaneTransform *trans2)
    213 {
    214     PS_PTR_CHECK_NULL(trans1, NULL);
    215     PS_PTR_CHECK_NULL(trans2, NULL);
    216 
    217     return(NULL);
    218 }
    219 
    220 
    221 // XXX: Must code this.
    222 // XXX: What about nRejIter?  Iterations?
    223 
    224 bool psPlaneTranformFit(
    225     psPlaneTransform *trans,
    226     const psArray *source,
    227     const psArray *dest,
    228     int nRejIter,
    229     float sigmaClip)
    230 {
    231     PS_PTR_CHECK_NULL(trans, NULL);
    232     PS_PTR_CHECK_NULL(source, NULL);
    233     PS_PTR_CHECK_NULL(dest, NULL);
    234 
    235     return(NULL);
    236 }
    237 
    238193
    239194/******************************************************************************

  • trunk/psLib/src/astro/psCoord.h

    r3450 r3540  
    1010*  @author GLG, MHPCC
    1111*
    12 *  @version $Revision: 1.28 $ $Name: not supported by cvs2svn $
    13 *  @date $Date: 2005-03-18 21:34:43 $
     12*  @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
     13*  @date $Date: 2005-03-29 19:41:56 $
    1414*
    1515*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    227227);
    228228
    229 // XXX: Doxygenate.
    230 psPlaneTransform *psPlaneTransformInvert(
    231     psPlaneTransform *out,
    232     const psPlaneTransform *in,
    233     psRegion *region,
    234     int nSamples
    235 );
    236 
    237 // XXX: Doxygenate.
    238 psPlaneTransform *psPlaneTransformCombine(
    239     psPlaneTransform *out,
    240     const psPlaneTransform *trans1,
    241     const psPlaneTransform *trans2
    242 );
    243 
    244 // XXX: Doxygenate.
    245 bool psPlaneTranformFit(
    246     psPlaneTransform *trans,
    247     const psArray *source,
    248     const psArray *dest,
    249     int nRejIter,
    250     float sigmaClip
    251 );
    252 
    253 
    254 
    255229/** Allocator for psSphereTransform
    256230 *

  • trunk/psLib/src/astronomy/psAstrometry.c

    r3497 r3540  
    88 *  @author GLG, MHPCC
    99 *
    10  *  @version $Revision: 1.60 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2005-03-24 22:36:16 $
     10 *  @version $Revision: 1.61 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2005-03-29 19:41:56 $
    1212 *
    1313 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    913913
    914914
    915 
    916 
    917 
    918 /*****************************************************************************
    919 psPlaneTransformInvertHEY(out, in, region, nSamples): this is an earlier
    920 version of this function which doesnot separate the code for building and
    921 solving the matrixequations.
    922  
    923 XXX: Delete this codeonce you know the other version works.
    924  *****************************************************************************/
    925 psPlaneTransform *psPlaneTransformInvertHEY(psPlaneTransform *out,
    926         const psPlaneTransform *in,
    927         psRegion *region,
    928         int nSamples)
    929 {
    930     PS_PTR_CHECK_NULL(in, NULL);
    931     if (isProjectionLinear((psPlaneTransform *) in)) {
    932         return(invertPlaneTransform((psPlaneTransform *) in));
    933     }
    934     psPlaneTransform *myPT = NULL;
    935     psPlane *inCoord = psPlaneAlloc();
    936     psPlane *outCoord = psPlaneAlloc();
    937 
    938     PS_PTR_CHECK_NULL(region, NULL);
    939     PS_INT_COMPARE(0, nSamples, NULL);
    940 
    941     // XXX: Is this correct?
    942     psS32 order = PS_MAX(in->x->nX, in->x->nY);
    943 
    944     //
    945     // Allocate a new psPlaneTransform if "out" is NULL, or has the wrong size.
    946     //
    947     if (out == NULL) {
    948         myPT = psPlaneTransformAlloc(order, order);
    949     } else {
    950         if (!((out->x->nX == order) &&
    951                 (out->x->nY == order) &&
    952                 (out->y->nX == order) &&
    953                 (out->y->nY == order))) {
    954             psFree(out);
    955             myPT = psPlaneTransformAlloc(order, order);
    956         } else {
    957             myPT = out;
    958         }
    959     }
    960     // XXX: Initialize myPT?
    961 
    962     //
    963     // Create fake polynomial to use in evaluation
    964     //
    965     psDPolynomial2D *fakePoly = psDPolynomial2DAlloc(order, order, PS_POLYNOMIAL_ORD);
    966     for (int i = 0; i < order; i++) {
    967         for (int j = 0; j < order; j++) {
    968             fakePoly->coeff[i][j] = 1.0; // Set all coeffecients to 1
    969             fakePoly->mask[i][j] = 1; // Mask all coefficients; unmask to evaluate
    970         }
    971     }
    972 
    973     //
    974     // Create a grid of xin,yin --> xout,yout
    975     //
    976     psVector *xIn = psVectorAlloc(nSamples * nSamples, PS_TYPE_F32);
    977     psVector *yIn = psVectorAlloc(nSamples * nSamples, PS_TYPE_F32);
    978     psVector *xOut = psVectorAlloc(nSamples * nSamples, PS_TYPE_F32);
    979     psVector *yOut = psVectorAlloc(nSamples * nSamples, PS_TYPE_F32);
    980 
    981     //
    982     // Initialize the grid of points
    983     //
    984     for (int yint = 0; yint < nSamples; yint++) {
    985         inCoord->y = region->y0 + ((psF32) yint) * ((region->y1 - region->y0) / ((psF32) nSamples));
    986         for (int xint = 0; xint < nSamples; xint++) {
    987             inCoord->x = region->x0 + ((psF32) xint) * ((region->x1 - region->x0) / ((psF32) nSamples));
    988 
    989             (void)psPlaneTransformApply(outCoord, in, inCoord);
    990             xOut->data.F32[yint*nSamples + xint] = inCoord->x;
    991             yOut->data.F32[yint*nSamples + xint] = inCoord->y;
    992             xIn->data.F32[yint*nSamples + xint] = outCoord->x;
    993             yIn->data.F32[yint*nSamples + xint] = outCoord->y;
    994         }
    995     }
    996 
    997     //
    998     // Initialise the matrix and vectors
    999     //
    1000     psS32 nCoeff = order * (order + 1) / 2; // Number of polynomial coefficients
    1001     psImage *matrix = psImageAlloc(nCoeff, nCoeff, PS_TYPE_F64); // Matrix for solution
    1002     psVector *xVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in x
    1003     psVector *yVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in y
    1004     for (psS32 i = 0; i < nCoeff; i++) {
    1005         for (psS32 j = 0; j < nCoeff; j++) {
    1006             matrix->data.F64[i][j] = 0.0;
    1007         }
    1008         xVector->data.F64[i] = 0.0;
    1009         yVector->data.F64[i] = 0.0;
    1010     }
    1011 
    1012     //
    1013     // Iterate over the grid points
    1014     //
    1015     for (psS32 g = 0; g < nSamples*nSamples; g++) {
    1016         // Iterate over the polynomial coefficients, accumulating the matrix and vectors
    1017 
    1018         for (psS32 i = 0, ijIndex = 0; i < order; i++) {
    1019             for (psS32 j = 0; j < order - i; j++, ijIndex++) {
    1020 
    1021                 fakePoly->mask[i][j] = 0;
    1022                 psF64 ijPoly = psDPolynomial2DEval(fakePoly, (psF64) xIn->data.F32[g], (psF64) yIn->data.F32[g]);
    1023                 fakePoly->mask[i][j] = 1;
    1024 
    1025                 for (psS32 m = 0, mnIndex = 0; m < order; m++) {
    1026                     for (psS32 n = 0; n < order - m; n++, mnIndex++) {
    1027                         fakePoly->mask[m][n] = 0;
    1028                         psF64 mnPoly = psDPolynomial2DEval(fakePoly, (psF64) xIn->data.F32[g], (psF64) yIn->data.F32[g]);
    1029                         fakePoly->mask[m][n] = 1;
    1030 
    1031                         matrix->data.F64[ijIndex][mnIndex] += ijPoly * mnPoly;
    1032                     }
    1033                 }
    1034 
    1035                 xVector->data.F64[ijIndex] += ijPoly * (psF64)xOut->data.F32[g];
    1036                 yVector->data.F64[ijIndex] += ijPoly * (psF64)yOut->data.F32[g];
    1037             }
    1038         }
    1039     }
    1040 
    1041     //
    1042     // Solution via LU Decomposition
    1043     //
    1044     psVector *permutation = psVectorAlloc(nCoeff, PS_TYPE_F64); // Permutation vector for LU Decomposition
    1045     psImage *luMatrix = psMatrixLUD(NULL, &permutation, matrix); // LU decomposed matrix
    1046     psVector *xSolution = psMatrixLUSolve(NULL, luMatrix, xVector, permutation); // Solution in x
    1047     psVector *ySolution = psMatrixLUSolve(NULL, luMatrix, yVector, permutation); // Solution in y
    1048 
    1049     //
    1050     // Stuff coefficients into transformation
    1051     //
    1052     for (psS32 i = 0, ijIndex = 0; i < order; i++) {
    1053         for (psS32 j = 0; j < order - i; j++, ijIndex++) {
    1054             myPT->x->coeff[i][j] = xSolution->data.F64[ijIndex];
    1055             myPT->y->coeff[i][j] = ySolution->data.F64[ijIndex];
    1056         }
    1057     }
    1058 
    1059     return(myPT);
    1060 }
    1061 
    1062 
    1063915/*****************************************************************************
    1064916multiplyCoeffs(trans1, trans2): Takes two 2-D polynomials as input and
    … …  
    1100952
    1101953/*****************************************************************************
    1102 psPlaneTransformCombineTmp(out, trans1, trans2)
     954psPlaneTransformCombine(out, trans1, trans2)
    1103955 
    1104956XXX: Much room for optimization.  Currently, we call the polyMultiply
    1105957routine far too many times.
    1106958 *****************************************************************************/
    1107 psPlaneTransform *psPlaneTransformCombineTmp(psPlaneTransform *out,
     959psPlaneTransform *psPlaneTransformCombine(psPlaneTransform *out,
    1108960        const psPlaneTransform *trans1,
    1109961        const psPlaneTransform *trans2)
    … …  
    12401092
    12411093/*****************************************************************************
    1242 psPlaneTranformFitTmp(trans, source, dest, nRejIter, sigmaClip)
     1094psPlaneTranformFit(trans, source, dest, nRejIter, sigmaClip)
    12431095 
    12441096XXX: What about nRejIter?  Iterations?
    12451097XXX: Use static vectors for internal data.
    12461098 *****************************************************************************/
    1247 bool psPlaneTranformFitTmp(psPlaneTransform *trans,
    1248                            const psArray *source,
    1249                            const psArray *dest,
    1250                            int nRejIter,
    1251                            float sigmaClip)
     1099bool psPlaneTranformFit(psPlaneTransform *trans,
     1100                        const psArray *source,
     1101                        const psArray *dest,
     1102                        int nRejIter,
     1103                        float sigmaClip)
    12521104{
    12531105    PS_PTR_CHECK_NULL(trans, NULL);
    … …  
    13601212
    13611213/*****************************************************************************
    1362 psPlaneTransformInvertTmp(out, in, region, nSamples)
     1214psPlaneTransformInvert(out, in, region, nSamples)
    13631215 
    13641216// XXX: Use static data structures.
    13651217 *****************************************************************************/
    1366 psPlaneTransform *psPlaneTransformInvertTmp(psPlaneTransform *out,
     1218psPlaneTransform *psPlaneTransformInvert(psPlaneTransform *out,
    13671219        const psPlaneTransform *in,
    13681220        psRegion *region,
    … …  
    14481300        }
    14491301    }
    1450     bool rc = psPlaneTranformFitTmp(myPT, inData, outData, 10, 100.0);
     1302    bool rc = psPlaneTranformFit(myPT, inData, outData, 10, 100.0);
    14511303
    14521304    psFree(inCoord);

  • trunk/psLib/src/astronomy/psAstrometry.h

    r3497 r3540  
    88*  @author GLG, MHPCC
    99*
    10 *  @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
    11 *  @date $Date: 2005-03-24 22:36:16 $
     10*  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
     11*  @date $Date: 2005-03-29 19:41:56 $
    1212*
    1313*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    537537// XXX: These functions don't belong here.  Will migrate to psCoords.c later.
    538538// XXX: Doxygenate.
    539 psPlaneTransform *psPlaneTransformInvertTmp(
     539psPlaneTransform *psPlaneTransformInvert(
    540540    psPlaneTransform *out,
    541541    const psPlaneTransform *in,
    … …  
    545545
    546546// XXX: Doxygenate.
    547 psPlaneTransform *psPlaneTransformCombineTmp(
     547psPlaneTransform *psPlaneTransformCombine(
    548548    psPlaneTransform *out,
    549549    const psPlaneTransform *trans1,
    … …  
    552552
    553553// XXX: Doxygenate.
    554 bool psPlaneTranformFitTmp(
     554bool psPlaneTranformFit(
    555555    psPlaneTransform *trans,
    556556    const psArray *source,

  • trunk/psLib/src/astronomy/psCoord.c

    r3470 r3540  
    1010*  @author GLG, MHPCC
    1111*
    12 *  @version $Revision: 1.59 $ $Name: not supported by cvs2svn $
    13 *  @date $Date: 2005-03-22 01:35:49 $
     12*  @version $Revision: 1.60 $ $Name: not supported by cvs2svn $
     13*  @date $Date: 2005-03-29 19:41:56 $
    1414*
    1515*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    191191    return (out);
    192192}
    193 
    194 // XXX: Must code this.
    195 psPlaneTransform *psPlaneTransformInvert(
    196     psPlaneTransform *out,
    197     const psPlaneTransform *in,
    198     psRegion *region,
    199     int nSamples)
    200 {
    201     PS_PTR_CHECK_NULL(in, NULL);
    202     PS_PTR_CHECK_NULL(region, NULL);
    203 
    204     return(NULL);
    205 }
    206 
    207 
    208 // XXX: Must code this.
    209 psPlaneTransform *psPlaneTransformCombine(
    210     psPlaneTransform *out,
    211     const psPlaneTransform *trans1,
    212     const psPlaneTransform *trans2)
    213 {
    214     PS_PTR_CHECK_NULL(trans1, NULL);
    215     PS_PTR_CHECK_NULL(trans2, NULL);
    216 
    217     return(NULL);
    218 }
    219 
    220 
    221 // XXX: Must code this.
    222 // XXX: What about nRejIter?  Iterations?
    223 
    224 bool psPlaneTranformFit(
    225     psPlaneTransform *trans,
    226     const psArray *source,
    227     const psArray *dest,
    228     int nRejIter,
    229     float sigmaClip)
    230 {
    231     PS_PTR_CHECK_NULL(trans, NULL);
    232     PS_PTR_CHECK_NULL(source, NULL);
    233     PS_PTR_CHECK_NULL(dest, NULL);
    234 
    235     return(NULL);
    236 }
    237 
    238193
    239194/******************************************************************************

  • trunk/psLib/src/astronomy/psCoord.h

    r3450 r3540  
    1010*  @author GLG, MHPCC
    1111*
    12 *  @version $Revision: 1.28 $ $Name: not supported by cvs2svn $
    13 *  @date $Date: 2005-03-18 21:34:43 $
     12*  @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
     13*  @date $Date: 2005-03-29 19:41:56 $
    1414*
    1515*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    227227);
    228228
    229 // XXX: Doxygenate.
    230 psPlaneTransform *psPlaneTransformInvert(
    231     psPlaneTransform *out,
    232     const psPlaneTransform *in,
    233     psRegion *region,
    234     int nSamples
    235 );
    236 
    237 // XXX: Doxygenate.
    238 psPlaneTransform *psPlaneTransformCombine(
    239     psPlaneTransform *out,
    240     const psPlaneTransform *trans1,
    241     const psPlaneTransform *trans2
    242 );
    243 
    244 // XXX: Doxygenate.
    245 bool psPlaneTranformFit(
    246     psPlaneTransform *trans,
    247     const psArray *source,
    248     const psArray *dest,
    249     int nRejIter,
    250     float sigmaClip
    251 );
    252 
    253 
    254 
    255229/** Allocator for psSphereTransform
    256230 *

  • trunk/psLib/src/dataManip/psMinimize.c

    r3476 r3540  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.107 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2005-03-22 21:52:49 $
     11 *  @version $Revision: 1.108 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2005-03-29 19:41:56 $
    1313 *
    1414 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    14141414XXX: Since this is an internal function, many of the parameter checks are
    14151415     redundant.
     1416XXX: Don't modify the psMinimization argument.
    14161417 *****************************************************************************/
    14171418#define PS_LINEMIN_MAX_ITERATIONS 30

  • trunk/psLib/src/dataManip/psStats.c

    r3529 r3540  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.121 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2005-03-28 22:49:33 $
     11 *  @version $Revision: 1.122 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2005-03-29 19:41:56 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii

  • trunk/psLib/src/math/psMinimize.c

    r3476 r3540  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.107 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2005-03-22 21:52:49 $
     11 *  @version $Revision: 1.108 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2005-03-29 19:41:56 $
    1313 *
    1414 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    14141414XXX: Since this is an internal function, many of the parameter checks are
    14151415     redundant.
     1416XXX: Don't modify the psMinimization argument.
    14161417 *****************************************************************************/
    14171418#define PS_LINEMIN_MAX_ITERATIONS 30

  • trunk/psLib/src/math/psStats.c

    r3529 r3540  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.121 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2005-03-28 22:49:33 $
     11 *  @version $Revision: 1.122 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2005-03-29 19:41:56 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
Note: See TracChangeset for help on using the changeset viewer.