Changeset 24382

Timestamp:

Jun 11, 2009, 6:09:52 PM (17 years ago)

Author:

eugene

Message:

modifications to support a two-pass check on the line range, limiting the initial ra,dec range to the region of the chip of interest using just the projection portion

Location:

trunk/magic/remove/src

Files:

• Line.c (modified) (1 diff)
• Line.h (modified) (1 diff)
• streaksastrom.c (modified) (1 diff)
• streaksastrom.h (modified) (1 diff)
• streaksextern.c (modified) (2 diffs)

trunk/magic/remove/src/Line.c

@@ -225 +225 @@
     vertices[2].x = numCols; vertices[2].y = numRows;
     vertices[3].x = 0;       vertices[3].y = numRows;
+
+    for (i = 0; i < 4 && found < 2; ++i)
+    {
+        boundLine.begin = vertices[i];
+        boundLine.end   = vertices[(i + 1) % 4];
+        if (LineIntercept (line, &boundLine, &tuple1, &tuple2, false, true))
+        {
+            if (found == 0)
+            {
+                clipLine.begin = tuple1;
+                ++found;
+            }
+            else if (tuple1.x != clipLine.begin.x || 
+                     tuple1.y != clipLine.begin.y)
+            {
+                clipLine.end = tuple1;
+                ++found;
+            }
+        }
+    }
+    
+    // If two endpoints are found, clip the line
+    
+    if (found > 1)
+    {
+        if (clipLine.begin.x <= clipLine.end.x)
+        {
+            line->begin = clipLine.begin;
+            line->end   = clipLine.end;
+        }
+        else
+        {
+            line->begin = clipLine.end;
+            line->end   = clipLine.begin;
+        }
+    }
+    return found > 1;
+}
+
+/** Clip the line between (minX,minY) and (maxX,maxY)
+
+    @param[in,out] line line to be clipped within the bounds
+    @param[in] minX minimum X (columns) for the line
+    @param[in] minY minimum Y (rows) for the line
+    @param[in] maxX maximum X (columns) for the line
+    @param[in] maxY maximum Y (rows) for the line
+    @return true if line overlaps the clip boundaries           */
+
+bool LineClipFull (Line *line, int minX, int minY, int maxX, int maxY)
+{
+    unsigned int i, found = 0;
+    Line boundLine, clipLine;
+    strkPt tuple1, tuple2, vertices[4];
+    vertices[0].x = minX; vertices[0].y = minY;
+    vertices[1].x = maxX; vertices[1].y = minY;
+    vertices[2].x = maxX; vertices[2].y = maxY;
+    vertices[3].x = minX; vertices[3].y = maxY;
 
     for (i = 0; i < 4 && found < 2; ++i)

trunk/magic/remove/src/Line.h

@@ -21 +21 @@
 extern bool LineClip (Line *line, int numCols, int numRows);
 
+/** Clip the line between (minX,minY) and (maxX,maxY)
+
+    @param[in,out] line line to be clipped within the bounds
+    @param[in] minX minimum X (columns) for the line
+    @param[in] minY minimum Y (rows) for the line
+    @param[in] maxX maximum X (columns) for the line
+    @param[in] maxY maximum Y (rows) for the line
+    @return true if line overlaps the clip boundaries           */
+
+extern bool LineClipFull (Line *line, int minX, int minY, int maxX, int maxY);
+
 /** Map a line to an image for its specified width and append as
     a list of pixel positions

trunk/magic/remove/src/streaksastrom.c

@@ -150 +150 @@
  
 bool
+SkyToLocal(strkPt *outPt, strkAstrom *astrom, double ra, double dec)
+{
+    // generate a local project using the RA, DEC of the 0,0 pixel of the chip as the
+    // projection center with the same plate scale as the nominal TP->Sky astrometry.
+
+    pmFPA *fpa   = (pmFPA *) astrom->fpa;
+    pmChip *chip = (pmChip *) astrom->chip;
+
+    // find the RA,DEC coords of the 0,0 pixel for this chip:
+
+    psPlane ptTP;
+    psSphere ptSky;
+
+    ptSky.r = ra;
+    ptSky.d = dec;
+    ptSky.rErr = 0.0;
+    ptSky.dErr = 0.0;
+
+    psProject(&ptTP, &ptSky, fpa->toSky);
+
+    outPt->x = ptTP.x;
+    outPt->y = ptTP.y;
+
+    return true;
+}
+
+bool
+LocalToSky(strkPt *outPt, strkAstrom *astrom, strkPt *inPt)
+{
+    // generate a local project using the RA, DEC of the 0,0 pixel of the chip as the
+    // projection center with the same plate scale as the nominal TP->Sky astrometry.
+
+    pmFPA *fpa   = (pmFPA *) astrom->fpa;
+    pmChip *chip = (pmChip *) astrom->chip;
+
+    // find the RA,DEC coords of the 0,0 pixel for this chip:
+
+    psPlane ptTP;
+    psSphere ptSky;
+
+    ptTP.x = inPt->x;
+    ptTP.y = inPt->y;
+    ptTP.xErr = 0.0;
+    ptTP.yErr = 0.0;
+
+    psDeproject(&ptSky, &ptTP, fpa->toSky);
+
+    outPt->x = ptSky.r;
+    outPt->y = ptSky.d;
+
+    return true;
+}
+
+bool
+componentBounds(int *minX, int *minY, int *maxX, int *maxY, strkAstrom *astrom, int numCols, int numRows)
+{
+    // find the bounds of the (padded) chip region in tangent-plane coordinates
+
+    pmFPA *fpa   = (pmFPA *) astrom->fpa;
+    pmChip *chip = (pmChip *) astrom->chip;
+
+    psPlane ptCH, ptFP, TPo, TPx, TPy;
+
+    // coordinate of the chip center:
+    ptCH.x = 0.5*numCols;
+    ptCH.y = 0.5*numRows;
+    psPlaneTransformApply(&ptFP, chip->toFPA, &ptCH);
+    psPlaneTransformApply(&TPo, fpa->toTPA,  &ptFP);
+
+    // coordinate of the chip center + dX/2:
+    ptCH.x = numCols;
+    ptCH.y = 0.5*numRows;
+    psPlaneTransformApply(&ptFP, chip->toFPA, &ptCH);
+    psPlaneTransformApply(&TPx, fpa->toTPA,  &ptFP);
+
+    // coordinate of the chip center + dY/2:
+    ptCH.x = 0.5*numCols;
+    ptCH.y = numRows;
+    psPlaneTransformApply(&ptFP, chip->toFPA, &ptCH);
+    psPlaneTransformApply(&TPy, fpa->toTPA,  &ptFP);
+
+    // half-lengths of the two sides in tangent-plane coords:
+    double xSize = hypot (TPx.x - TPo.x, TPx.y - TPo.y);
+    double ySize = hypot (TPy.x - TPo.x, TPy.y - TPo.y);
+    double radius = hypot (xSize, ySize);
+
+    // define the region encompassed by the radius with some padding:
+    *minX = TPo.x - 1.1*radius;
+    *minY = TPo.y - 1.1*radius;
+    *maxX = TPo.x + 1.1*radius;
+    *maxY = TPo.y + 1.1*radius;
+
+    return true;
+}
+
+bool
 skyToCell(strkPt *outPt, strkAstrom *astrom, double ra, double dec)
 {

trunk/magic/remove/src/streaksastrom.h

@@ -32 +32 @@
 extern void linearizeTransforms(strkAstrom *astrom);
 
+extern bool SkyToLocal(strkPt *outPt, strkAstrom *astrom, double ra, double dec);
+extern bool LocalToSky(strkPt *outPt, strkAstrom *astrom, strkPt *inPt);
+extern bool componentBounds(int *minX, int *minY, int *maxX, int *maxY, strkAstrom *astrom, int numCols, int numRows);
+
 #endif // STREAKS_ASTROM_H

trunk/magic/remove/src/streaksextern.c

@@ -36 +36 @@
     StreakPixels *pixels = psArrayAllocEmpty (1024);
     int streaksOnComponent = 0;
+
+    int minX, minY, maxX, maxY;
+
+    // find the chip dimensions in the tangent-plane coordinates (length of hypotenuse)
+    componentBounds (&minX, &minY, &maxX, &maxY, astrom, numCols, numRows);
+
     for (i = 0; i != streaks->size; ++i)
     {
@@ -41 +47 @@
 
         line.width = streaks->list[i].width;
-        if (skyToCell (&line.begin, astrom,
-                       streaks->list[i].ra1, streaks->list[i].dec1) &&
-            skyToCell (&line.end, astrom,
-                       streaks->list[i].ra2, streaks->list[i].dec2) &&
+
+        /* Use tangent plane coordinates to narrow down the ra,dec range of the line closer to
+         * the chip boundaries.  Use these new ra,dec positions to generate the line on the
+         * chip using the full non-linear astrometry */
+           
+        // project the ends of the line using a linear projection centered on the chip center:
+        Line full;
+        SkyToLocal (&full.begin, astrom, streaks->list[i].ra1, streaks->list[i].dec1);
+        SkyToLocal (&full.end,   astrom, streaks->list[i].ra2, streaks->list[i].dec2);
+
+        // clip the line to a square box with diameter = hypotenuse of the chip image centerd
+        // on the chip center in tangent-plane coordinates.  skip the rest of this streak if
+        // the line does not intersect this region
+        if (!LineClipFull (&full, minX, minY, maxX, maxY)) {
+            continue;
+        }
+
+        // convert the end points back into ra, dec pairs:
+        strkPt sky1, sky2;
+        LocalToSky (&sky1, astrom, &full.begin);
+        LocalToSky (&sky2, astrom, &full.end);
+
+        if (skyToCell (&line.begin, astrom, sky1.x, sky1.y) &&
+            skyToCell (&line.end,   astrom, sky2.x, sky2.y) &&
             LineClip (&line, numCols, numRows))
         {