IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5542


Ignore:
Timestamp:
Nov 18, 2005, 9:39:29 AM (20 years ago)
Author:
gusciora
Message:

Asserts for psPolynomial, andvarious psCoord mods.

Location:
trunk/psLib/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r5294 r5542  
    1010*  @author GLG, MHPCC
    1111*
    12 *  @version $Revision: 1.88 $ $Name: not supported by cvs2svn $
    13 *  @date $Date: 2005-10-12 21:02:20 $
     12*  @version $Revision: 1.89 $ $Name: not supported by cvs2svn $
     13*  @date $Date: 2005-11-18 19:39:29 $
    1414*
    1515*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    5757"transform" is linear.
    5858 
    59 This program assumes that the inverse of the following linear equations:
    60         X2 = A + (B * X1) + (C * Y1);
    61         Y2 = D + (E * X1) + (F * Y1);
    62 is
    63         Y1 = (Y2 - ((E/B) * X2) - D + ((E*A)/B)) / (F - ((C*E)/B));
    64         X1 = (Y2 - ((F/C) * X2) - D + ((F*A)/C)) / (E - ((F*B)/C));
    65 or
    66  X1 = (-D + ((F*A)/C)) / (E - ((F*B)/C)) +
    67       (X2 * -((F/C) / (E - ((F*B)/C)))) +
    68       (Y2 * (1.0 / (E - ((F*B)/C))));
    69  Y1 = (-D + ((E*A)/B))/(F - ((C*E)/B)) +
    70       (X2 * -((E/B) / (F - ((C*E)/B)))) +
    71       (Y2 * (1.0 / (F - ((C*E)/B))));
    72  
    7359XXX: Since there is now a general psPlaneTransformInvert() function, we
    7460should rename this.
    7561 
    76 XXX: Use the ADD version which is based on determinants.
     62To derive this transformation, start with the simple 2-by-2 matrix inversion
     63based on discriminants.  This will invert the
     64    (x_2, y_2) = matrix(a b c d) * vector (x, y)
     65where there are no constant terms.  Then you substitute
     66    x_2 = x_1 - e
     67    y_2 = y_1 - f
     68for (x_2, y_2) to get the desired inverse.
    7769 *****************************************************************************/
    7870psPlaneTransform *p_psPlaneTransformLinearInvert(psPlaneTransform *transform)
    7971{
    80     PS_ASSERT_PTR_NON_NULL(transform, 0);
    81     PS_ASSERT_PTR_NON_NULL(transform->x, 0);
    82     PS_ASSERT_PTR_NON_NULL(transform->y, 0);
    83 
    84     psF64 A = 0.0;
    85     psF64 B = 0.0;
    86     psF64 C = 0.0;
    87     psF64 D = 0.0;
    88     psF64 E = 0.0;
    89     psF64 F = 0.0;
    90 
    91     A = transform->x->coeff[0][0];
    92     if (transform->x->nX >= 1) {
    93         B = transform->x->coeff[1][0];
    94     }
    95     if (transform->x->nY >= 1) {
    96         C = transform->x->coeff[0][1];
    97     }
    98     D = transform->y->coeff[0][0];
    99     if (transform->y->nX >= 1) {
    100         E = transform->y->coeff[1][0];
    101     }
    102     if (transform->y->nY >= 1) {
    103         F = transform->y->coeff[0][1];
    104     }
    105 
    106     psPlaneTransform *out = psPlaneTransformAlloc(2, 2);
    107 
    108     /* This is sample code from IfA.  It didn't work initially, and I did not
    109        spend any time debugging it.
    110 
     72    PS_ASSERT_PTR_NON_NULL(transform, NULL);
     73    PS_ASSERT_PTR_NON_NULL(transform->x, NULL);
     74    PS_ASSERT_PTR_NON_NULL(transform->y, NULL);
     75    if ((transform->x->nX < 1) ||
     76            (transform->x->nY < 1) ||
     77            (transform->y->nX < 1) ||
     78            (transform->y->nY < 1)) {
     79        psLogMsg(__func__, PS_LOG_WARN, "WARNING: input transform is not invertible.");
     80        return(NULL);
     81    }
     82
     83    psPlaneTransform *out = psPlaneTransformAlloc(1, 1);
     84
     85    if (0) {
     86        //This is sample code from IfA.  It didn't work initially, and I did not
     87        //spend any time debugging it.
    11188        psF64 a = transform->x->coeff[1][0];
    11289        psF64 b = transform->x->coeff[0][1];
    … …  
    11895        psF64 invDet = 1.0 / (a * d - b * c); // Inverse of the determinant
    11996
    120         // Not entirely sure why this works, but it appears to do so....................................!
     97        // Not entirely sure why this works, but it appears to do so
    12198        out->x->coeff[1][0] = invDet * a;
    12299        out->x->coeff[0][1] = - invDet * b;
    … …  
    126103        out->x->coeff[0][0] = - invDet * (d * e + c * f);
    127104        out->y->coeff[0][0] = - invDet * (b * e + a * f);
    128     */
    129     out->x->coeff[0][0] = (-D + ((F*A)/C)) / (E - ((F*B)/C));
    130     out->x->coeff[1][0] = -(F/C) / (E - ((F*B)/C));
    131     out->x->coeff[0][1] =  1.0 / (E - ((F*B)/C));
    132     out->y->coeff[0][0] = (-D + ((E*A)/B)) / (F - ((C*E)/B));
    133     out->y->coeff[1][0] = -(E/B) / (F - ((C*E)/B));
    134     out->y->coeff[0][1] =  1.0 / (F - ((C*E)/B));
     105    }
     106
     107    if (1) {
     108        psF64 A = transform->x->coeff[1][0];
     109        psF64 B = transform->x->coeff[0][1];
     110        psF64 C = transform->y->coeff[1][0];
     111        psF64 D = transform->y->coeff[0][1];
     112        psF64 E = transform->x->coeff[0][0];
     113        psF64 F = transform->y->coeff[0][0];
     114        psF64 invDet = 1.0 / (A * D - B * C);
     115
     116        out->x->coeff[1][0] = invDet *  D;
     117        out->x->coeff[0][1] = invDet * -B;
     118        out->y->coeff[1][0] = invDet * -C;
     119        out->y->coeff[0][1] = invDet *  A;
     120
     121        //
     122        // This comes from substituting the (x_2 = x_1 - e) and (y_2 = y_1 - f) equations:
     123        //
     124        out->x->coeff[0][0] = invDet * (B * F - D * E);
     125        out->y->coeff[0][0] = invDet * (E * C - A * F);
     126
     127    }
     128    //    printf("HMMM: out->x: (%f %f %f)\n", out->x->coeff[0][0], out->x->coeff[1][0], out->x->coeff[0][1]);
     129    //    printf("HMMM: out->y: (%f %f %f)\n", out->y->coeff[0][0], out->y->coeff[1][0], out->y->coeff[0][1]);
    135130
    136131    return(out);
    … …  
    143138 
    144139Returns:
    145     1: if linear
    146     0: otherwise
    147  
    148 XXX: This should be a psBool
     140    true: if linear
     141    false: otherwise
    149142 *****************************************************************************/
    150 psS32 p_psIsProjectionLinear(psPlaneTransform *transform)
    151 {
    152     PS_ASSERT_PTR_NON_NULL(transform, 0);
    153     PS_ASSERT_PTR_NON_NULL(transform->x, 0);
    154     PS_ASSERT_PTR_NON_NULL(transform->y, 0);
     143psBool p_psIsProjectionLinear(psPlaneTransform *transform)
     144{
     145    PS_ASSERT_PTR_NON_NULL(transform, false);
     146    PS_ASSERT_PTR_NON_NULL(transform->x, false);
     147    PS_ASSERT_PTR_NON_NULL(transform->y, false);
    155148
    156149    for (psS32 i=0;i<(1 + transform->x->nX);i++) {
    … …  
    160153                        ((i == 0) && (j == 1)) ||
    161154                        ((i == 1) && (j == 0)))) {
    162                     return(0);
     155                    return(false);
    163156                }
    164157            }
    … …  
    172165                        ((i == 0) && (j == 1)) ||
    173166                        ((i == 1) && (j == 0)))) {
    174                     return(0);
     167                    return(false);
    175168                }
    176169            }
    … …  
    178171    }
    179172
    180     return(1);
     173    return(true);
    181174}
    182175
    … …  
    215208
    216209    psPlaneTransform *pt = psAlloc(sizeof(psPlaneTransform));
    217     pt->x = psPolynomial2DAlloc(n1-1, n2-1, PS_POLYNOMIAL_ORD);
    218     pt->y = psPolynomial2DAlloc(n1-1, n2-1, PS_POLYNOMIAL_ORD);
     210    pt->x = psPolynomial2DAlloc(n1, n2, PS_POLYNOMIAL_ORD);
     211    pt->y = psPolynomial2DAlloc(n1, n2, PS_POLYNOMIAL_ORD);
    219212
    220213    psMemSetDeallocator(pt, (psFreeFunc) planeTransformFree);
    … …  
    243236}
    244237
    245 psPlane* psPlaneTransformApply(psPlane* out,
    246                                const psPlaneTransform* transform,
    247                                const psPlane* coords)
     238psPlane* psPlaneTransformApply(
     239    psPlane* out,
     240    const psPlaneTransform* transform,
     241    const psPlane* coords)
    248242{
    249243    PS_ASSERT_PTR_NON_NULL(transform, NULL);
    … …  
    256250    }
    257251
    258     out->x = psPolynomial2DEval(
    259                  transform->x,
    260                  coords->x,
    261                  coords->y
    262              );
    263 
    264     out->y = psPolynomial2DEval(
    265                  transform->y,
    266                  coords->x,
    267                  coords->y
    268              );
     252    out->x = psPolynomial2DEval(transform->x, coords->x, coords->y);
     253    out->y = psPolynomial2DEval(transform->y, coords->x, coords->y);
    269254
    270255    return (out);
    … …  
    383368
    384369
    385 psPlane* psProject(const psSphere* coord,
    386                    const psProjection* projection)
     370psPlane* p_psProject(
     371    psPlane *outPlane,
     372    const psSphere* coord,
     373    const psProjection* projection)
    387374{
    388375    PS_ASSERT_PTR_NON_NULL(coord, NULL);
    … …  
    393380
    394381    // Allocate return value
    395     psPlane* out = psPlaneAlloc();
     382    psPlane* out = NULL;
     383    if (outPlane == NULL) {
     384        out = psPlaneAlloc();
     385    } else {
     386        out = outPlane;
     387    }
    396388
    397389    // Convert to projection spherical coordinate system
    … …  
    435427}
    436428
    437 psSphere* psDeproject(const psPlane* coord,
    438                       const psProjection* projection)
     429psPlane* psProject(const psSphere* coord,
     430                   const psProjection* projection)
     431{
     432    return(p_psProject(NULL, coord, projection));
     433}
     434
     435psSphere* p_psDeproject(
     436    psSphere *outSphere,
     437    const psPlane* coord,
     438    const psProjection* projection)
    439439{
    440440    PS_ASSERT_PTR_NON_NULL(coord, NULL);
    … …  
    445445
    446446    // Allocate return sphere structure
    447     psSphere* out = psSphereAlloc();
     447    psSphere *out = NULL;
     448    if (outSphere == NULL) {
     449        out = psSphereAlloc();
     450    } else {
     451        out = outSphere;
     452        // XXX: Do a memory checkpointer.
     453    }
    448454
    449455    // Remove plate scales
    … …  
    490496}
    491497
     498psSphere* psDeproject(const psPlane* coord,
     499                      const psProjection* projection)
     500{
     501    return(p_psDeproject(NULL, coord, projection));
     502}
     503
     504
    492505/*****************************************************************************
    493506multiplyDPoly2D(trans1, trans2): Takes two 2-D polynomials as input and
    … …  
    813826{
    814827    PS_ASSERT_PTR_NON_NULL(in, NULL);
     828    PS_ASSERT_PTR_NON_NULL(in->x, NULL);
     829    PS_ASSERT_PTR_NON_NULL(in->y, NULL);
     830    if ((in->x->nX < 1) ||
     831            (in->x->nY < 1) ||
     832            (in->y->nX < 1) ||
     833            (in->y->nY < 1)) {
     834        psLogMsg(__func__, PS_LOG_WARN, "WARNING: input transform is not invertible.");
     835        return(NULL);
     836    }
    815837    //
    816838    // If the transform is linear, then invert it exactly and return.

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

    r4966 r5542  
    1111*  @author GLG, MHPCC
    1212*
    13 *  @version $Revision: 1.46 $ $Name: not supported by cvs2svn $
    14 *  @date $Date: 2005-09-07 23:51:20 $
     13*  @version $Revision: 1.47 $ $Name: not supported by cvs2svn $
     14*  @date $Date: 2005-11-18 19:39:29 $
    1515*
    1616*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    290290);
    291291
     292psPlane* p_psProject(
     293    psPlane *out,
     294    const psSphere* coord,             ///< coordinate to project
     295    const psProjection* projection     ///< parameters of the projection
     296);
     297
    292298/** Reverse projection of a linear coordinate to a spherical coordinate system
    293299 *
    … …  
    295301 */
    296302psSphere* psDeproject(
     303    const psPlane* coord,              ///< coordinate to project
     304    const psProjection* projection     ///< parameters of the projection
     305);
     306
     307/** Private reverse projection of a linear coordinate to a spherical coordinate system
     308 *
     309 *  @return psPlane*    projected coordinate
     310 */
     311psSphere* p_psDeproject(
     312    psSphere *outSphere,
    297313    const psPlane* coord,              ///< coordinate to project
    298314    const psProjection* projection     ///< parameters of the projection
    … …  
    314330 *  the order of the projection
    315331*/
    316 psS32 p_psIsProjectionLinear(
     332psBool p_psIsProjectionLinear(
    317333    psPlaneTransform *transform        ///< transform to test for linearity
    318334);
     335
    319336
    320337/** inverts a given transformation.

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

    r5294 r5542  
    77*  polynomials.  It also contains a Gaussian functions.
    88*
    9 *  @version $Revision: 1.130 $ $Name: not supported by cvs2svn $
    10 *  @date $Date: 2005-10-12 21:02:20 $
     9*  @version $Revision: 1.131 $ $Name: not supported by cvs2svn $
     10*  @date $Date: 2005-11-18 19:39:29 $
    1111*
    1212*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    644644                                     psPolynomialType type)
    645645{
    646     PS_ASSERT_INT_POSITIVE(nX, NULL);
    647     PS_ASSERT_INT_POSITIVE(nY, NULL);
     646    PS_ASSERT_INT_NONNEGATIVE(nX, NULL);
     647    PS_ASSERT_INT_NONNEGATIVE(nY, NULL);
    648648
    649649    unsigned int x = 0;
    … …  
    682682                                     psPolynomialType type)
    683683{
    684     PS_ASSERT_INT_POSITIVE(nX, NULL);
    685     PS_ASSERT_INT_POSITIVE(nY, NULL);
    686     PS_ASSERT_INT_POSITIVE(nZ, NULL);
     684    PS_ASSERT_INT_NONNEGATIVE(nX, NULL);
     685    PS_ASSERT_INT_NONNEGATIVE(nY, NULL);
     686    PS_ASSERT_INT_NONNEGATIVE(nZ, NULL);
    687687
    688688    unsigned int x = 0;
    … …  
    731731                                     psPolynomialType type)
    732732{
    733     PS_ASSERT_INT_POSITIVE(nX, NULL);
    734     PS_ASSERT_INT_POSITIVE(nY, NULL);
    735     PS_ASSERT_INT_POSITIVE(nZ, NULL);
    736     PS_ASSERT_INT_POSITIVE(nT, NULL);
     733    PS_ASSERT_INT_NONNEGATIVE(nX, NULL);
     734    PS_ASSERT_INT_NONNEGATIVE(nY, NULL);
     735    PS_ASSERT_INT_NONNEGATIVE(nZ, NULL);
     736    PS_ASSERT_INT_NONNEGATIVE(nT, NULL);
    737737
    738738    unsigned int x = 0;
Note: See TracChangeset for help on using the changeset viewer.