IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 21025


Ignore:
Timestamp:
Dec 17, 2008, 9:12:34 AM (17 years ago)
Author:
Paul Price
Message:

Adding assertions to protect from SEGV.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/astrom/pmAstrometryWCS.c

    r20720 r21025  
    77 *  @author EAM, IfA
    88 *
    9  *  @version $Revision: 1.31 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2008-11-13 19:38:16 $
     9 *  @version $Revision: 1.32 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2008-12-17 19:12:34 $
    1111 *
    1212 *  Copyright 2006 Institute for Astronomy, University of Hawaii
    … …  
    8080
    8181// convert toFPA / toSky components to traditional WCS
    82 // we require the header to have NAXIS1,NAXIS2, the field of the FPA 
     82// we require the header to have NAXIS1,NAXIS2, the field of the FPA
    8383// the center of the TPA/Sky projection is 0.5*(NAXIS1,NAXIS2)
    8484bool pmAstromReadBilevelMosaic (pmFPA *fpa, const psMetadata *header)
    … …  
    9696
    9797    if (!status1 || !status2) {
    98         Nx = psMetadataLookupS32 (&status1, header, "ZNAXIS1");
    99         Ny = psMetadataLookupS32 (&status2, header, "ZNAXIS2");
     98        Nx = psMetadataLookupS32 (&status1, header, "ZNAXIS1");
     99        Ny = psMetadataLookupS32 (&status2, header, "ZNAXIS2");
    100100    }
    101101
    102102    if (!status1 || !status2) {
    103         psFree (wcs);
     103        psFree (wcs);
    104104        psError(PS_ERR_UNKNOWN, false, "missing required FPA size in header");
    105105        return false;
    … …  
    465465        }
    466466
    467         // apply the exiting fromTPA transformation to make the new toFPA consistent with the toTPA layter
    468         // XXX this only works if toTPA is at most a linear transformation
     467        // apply the exiting fromTPA transformation to make the new toFPA consistent with the toTPA layter
     468        // XXX this only works if toTPA is at most a linear transformation
    469469        psPlaneTransform *toFPAnew = psPlaneTransformAlloc(toFPA->x->nX, toFPA->x->nY);
    470         for (int i = 0; i <= toFPA->x->nX; i++) {
    471           for (int j = 0; j <= toFPA->x->nY; j++) {
    472             double f1 = toFPA->x->coeffMask[i][j] ? 0.0 : fpa->fromTPA->x->coeff[1][0]*toFPA->x->coeff[i][j];
    473             double f2 = toFPA->y->coeffMask[i][j] ? 0.0 : fpa->fromTPA->x->coeff[0][1]*toFPA->y->coeff[i][j];
    474             toFPAnew->x->coeff[i][j] = f1 + f2;
    475 
    476             double g1 = toFPA->x->coeffMask[i][j] ? 0.0 : fpa->fromTPA->y->coeff[1][0]*toFPA->x->coeff[i][j];
    477             double g2 = toFPA->y->coeffMask[i][j] ? 0.0 : fpa->fromTPA->y->coeff[0][1]*toFPA->y->coeff[i][j];
    478             toFPAnew->y->coeff[i][j] = g1 + g2;
    479           }
    480         }
    481         toFPAnew->x->coeff[0][0] += fpa->fromTPA->x->coeff[0][0];
    482         toFPAnew->y->coeff[0][0] += fpa->fromTPA->y->coeff[0][0];
    483 
    484         psFree (toFPA);
    485         toFPA = toFPAnew;
     470        for (int i = 0; i <= toFPA->x->nX; i++) {
     471          for (int j = 0; j <= toFPA->x->nY; j++) {
     472            double f1 = toFPA->x->coeffMask[i][j] ? 0.0 : fpa->fromTPA->x->coeff[1][0]*toFPA->x->coeff[i][j];
     473            double f2 = toFPA->y->coeffMask[i][j] ? 0.0 : fpa->fromTPA->x->coeff[0][1]*toFPA->y->coeff[i][j];
     474            toFPAnew->x->coeff[i][j] = f1 + f2;
     475
     476            double g1 = toFPA->x->coeffMask[i][j] ? 0.0 : fpa->fromTPA->y->coeff[1][0]*toFPA->x->coeff[i][j];
     477            double g2 = toFPA->y->coeffMask[i][j] ? 0.0 : fpa->fromTPA->y->coeff[0][1]*toFPA->y->coeff[i][j];
     478            toFPAnew->y->coeff[i][j] = g1 + g2;
     479          }
     480        }
     481        toFPAnew->x->coeff[0][0] += fpa->fromTPA->x->coeff[0][0];
     482        toFPAnew->y->coeff[0][0] += fpa->fromTPA->y->coeff[0][0];
     483
     484        psFree (toFPA);
     485        toFPA = toFPAnew;
    486486
    487487        // adjust reference pixel for new toSky reference coordinate
    … …  
    492492        psSphere *sky = psSphereAlloc();
    493493
    494         sky->r = toSky->R;
    495         sky->d = toSky->D;
     494        sky->r = toSky->R;
     495        sky->d = toSky->D;
    496496        psProject (tp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection
    497497        psPlaneTransformApply (fpOld, fpa->fromTPA, tp);
    498498
    499         sky->r = fpa->toSky->R;
    500         sky->d = fpa->toSky->D;
     499        sky->r = fpa->toSky->R;
     500        sky->d = fpa->toSky->D;
    501501        psProject (tp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection
    502502        psPlaneTransformApply (fpNew, fpa->fromTPA, tp);
    … …  
    595595pmAstromWCS *pmAstromWCSfromFPA (const pmFPA *fpa, const pmChip *chip, double tol)
    596596{
     597    PS_ASSERT_PTR_NON_NULL(fpa, NULL);
     598    PS_ASSERT_PTR_NON_NULL(chip, NULL);
     599    PS_ASSERT_PTR_NON_NULL(fpa->toFPA, NULL);
     600    PS_ASSERT_PTR_NON_NULL(fpa->toTPA, NULL);
     601
    597602    // XXX require chip->toFPA->x->nX == chip->toFPA->x->nY
    598603    // XXX require chip->toFPA->y->nX == chip->toFPA->y->nY
    … …  
    607612    // create a temporary transform which combines toTPA and toFPA.  allow toTPA to have 0th
    608613    // and 1st order terms
    609    
     614
    610615    // XXX require fpa->toTPA->x->nX == fpa->toTPA->x->nY
    611616    // XXX require fpa->toTPA->y->nX == fpa->toTPA->y->nY
    … …  
    618623
    619624    for (int i = 0; i <= toTPA->x->nX; i++) {
    620         for (int j = 0; j <= toTPA->x->nY; j++) {
    621             double f1 = chip->toFPA->x->coeffMask[i][j] ? 0.0 : fpa->toTPA->x->coeff[1][0]*chip->toFPA->x->coeff[i][j];
    622             double f2 = chip->toFPA->y->coeffMask[i][j] ? 0.0 : fpa->toTPA->x->coeff[0][1]*chip->toFPA->y->coeff[i][j];
    623             toTPA->x->coeff[i][j] = f1 + f2;
    624 
    625             double g1 = chip->toFPA->x->coeffMask[i][j] ? 0.0 : fpa->toTPA->y->coeff[1][0]*chip->toFPA->x->coeff[i][j];
    626             double g2 = chip->toFPA->y->coeffMask[i][j] ? 0.0 : fpa->toTPA->y->coeff[0][1]*chip->toFPA->y->coeff[i][j];
    627             toTPA->y->coeff[i][j] = g1 + g2;
    628         }
     625        for (int j = 0; j <= toTPA->x->nY; j++) {
     626            double f1 = chip->toFPA->x->coeffMask[i][j] ? 0.0 : fpa->toTPA->x->coeff[1][0]*chip->toFPA->x->coeff[i][j];
     627            double f2 = chip->toFPA->y->coeffMask[i][j] ? 0.0 : fpa->toTPA->x->coeff[0][1]*chip->toFPA->y->coeff[i][j];
     628            toTPA->x->coeff[i][j] = f1 + f2;
     629
     630            double g1 = chip->toFPA->x->coeffMask[i][j] ? 0.0 : fpa->toTPA->y->coeff[1][0]*chip->toFPA->x->coeff[i][j];
     631            double g2 = chip->toFPA->y->coeffMask[i][j] ? 0.0 : fpa->toTPA->y->coeff[0][1]*chip->toFPA->y->coeff[i][j];
     632            toTPA->y->coeff[i][j] = g1 + g2;
     633        }
    629634    }
    630635    toTPA->x->coeff[0][0] += fpa->toTPA->x->coeff[0][0];
    … …  
    641646    psPlane *center = psPlaneTransformGetCenter (toTPA, tol);
    642647    if (!center) {
    643         psError(PS_ERR_UNKNOWN, false, "Unable to solve for TPA center.");
    644         psFree (toTPA);
    645         psFree (wcs);
    646         return NULL;
     648        psError(PS_ERR_UNKNOWN, false, "Unable to solve for TPA center.");
     649        psFree (toTPA);
     650        psFree (wcs);
     651        return NULL;
    647652    }
    648653
    … …  
    703708    psPlane *center = psPlaneTransformGetCenter (chip->toFPA, tol);
    704709    if (!center) {
    705         psError(PS_ERR_UNKNOWN, false, "Unable to solve for TPA center.");
    706         psFree (wcs);
    707         return NULL;
     710        psError(PS_ERR_UNKNOWN, false, "Unable to solve for TPA center.");
     711        psFree (wcs);
     712        return NULL;
    708713    }
    709714
    … …  
    744749    psPlane *center = psPlaneTransformGetCenter (fpa->toTPA, tol);
    745750    if (!center) {
    746         psError(PS_ERR_UNKNOWN, false, "Unable to solve for TPA center.");
    747         psFree (wcs);
    748         return NULL;
     751        psError(PS_ERR_UNKNOWN, false, "Unable to solve for TPA center.");
     752        psFree (wcs);
     753        return NULL;
    749754    }
    750755
    … …  
    801806
    802807/*****
    803  
     808
    804809For mosaic astrometry, we need to have a starting set of projection terms in which the
    805810chip-to-FPA terms result in a fixed physical unit on the focal plane (eg, pixels or
    … …  
    814819the chip-to-FPA scaling (eg, pc11) to match the variations in the effective plate scale for
    815820each chip (eg, cdelt1).  Thus, we need to carry around both the
    816  
     821
    817822*****/
    818823
    … …  
    821826   L,M: coord on the focal plane (pixels)
    822827   P,Q: coord in the tangent plane (microns or mm?)
    823    R,D: coord on the sky 
    824  
     828   R,D: coord on the sky
     829
    825830   this function creates WCS terms which convert directly from chip to sky.
    826831   this function requires a linear, unrotated toTPA distortion term
    … …  
    830835
    831836/* at this point, we have extracted from the header the WCS terms in the form of a polynomial,
    832  * wcs->trans, which will convert X,Y in pixels to L,M in degrees.  we also have the following 
     837 * wcs->trans, which will convert X,Y in pixels to L,M in degrees.  we also have the following
    833838 * elements defined:
    834839 * type (projection type)
    835840 * crval1,2 (in RA,DEC degrees)
    836  * crpix1,2 
     841 * crpix1,2
    837842 * cdelt1,2 (in degrees / pixel)
    838843 * pixelScale (microns / pixel)
    839  * 
     844 *
    840845 * now we convert wcs->trans to toFPA, which is different from wcs->trans in 3 important ways:
    841846 * 1) the output is in microns (not degrees): divide by cdelt1,2
Note: See TracChangeset for help on using the changeset viewer.