Index: trunk/magic/remove/src/Line.c
===================================================================
--- trunk/magic/remove/src/Line.c	(revision 21156)
+++ trunk/magic/remove/src/Line.c	(revision 24344)
@@ -16,4 +16,13 @@
 {
     double temp = *first;
+    *first = *second;
+    *second = temp;
+}
+
+/** Internal routine to swap integer values */
+
+void SwapInt (int* first, int* second)
+{
+    int temp = *first;
     *first = *second;
     *second = temp;
@@ -255,15 +264,53 @@
 }
 
-/** Map a line to an image for its specified width and store as a list
-    of pixel positions
+/** Move a line by the specified X and Y offsets
+    @param[in,out] line Line to move by X and Y offsets
+    @param[in] xOffset X shift applied to both endpoints
+    @param[in] yOffset Y shift applied to both endpoints        */
+
+void LineMove (Line *line, double xOffset, double yOffset)
+{
+    line->begin.x += xOffset;
+    line->begin.y += yOffset;
+    line->end.x   += xOffset;
+    line->end.y   += yOffset;
+}
+
+/** Return the maximum bounds between the line endpoints and
+    current bounds
+    @param[in] line Line endpoints to compare
+    @param[in,out] xMin minimum X value to update
+    @param[in,out] xMax maximum X value to update
+    @param[in,out] yMin minimum Y value to update
+    @param[in,out] yMax maximum Y value to update               */
+
+void MaxBounds (Line *line, int *xMin, int *xMax,  int *yMin, int *yMax)
+{
+    if (line->begin.x < *xMin) *xMin = (int) floor (line->begin.x);
+    if (line->end.x   < *xMin) *xMin = (int) floor (line->end.x);
+    if (line->begin.y < *yMin) *yMin = (int) floor (line->begin.y);
+    if (line->end.y   < *yMin) *yMin = (int) floor (line->end.y);
+
+    if (line->begin.x > *xMax) *xMax = (int) ceil (line->begin.x);
+    if (line->end.x   > *xMax) *xMax = (int) ceil (line->end.x);
+    if (line->begin.y > *yMax) *yMax = (int) ceil (line->begin.y);
+    if (line->end.y   > *yMax) *yMax = (int) ceil (line->end.y);
+}
+
+/** Map a line to an image for its specified width and store as
+    a list of pixel positions
 
     @param[out] pixels list of PixelPos pointers corresponding
                        based on the line settings
-    @param[in] line Line to map to pixels                       */
-
-void PixelsFromLine (StreakPixels* pixels, Line *line)
-{
+    @param[in] line Line to map to pixels   
+    @param[in] numCols maximum X (columns) for the line
+    @param[in] numRows maximum Y (rows) for the line            */
+
+void PixelsFromLine (StreakPixels* pixels, Line *line, int numCols, int numRows)
+{
+    Line offsetLine;
     PixelPos *pixel;
-    double slope, xOffset, yOffset, xMid, yMid, xBegin, yBegin, xEnd, yEnd, x, y;
+    double slope, xOffset, yOffset, xMid, yMid;
+    int x, y, xBegin = numCols, yBegin = numRows, xEnd = 0, yEnd = 0;
 
     // Extract the endpoints
@@ -280,7 +327,31 @@
     double dr = sqrt (dx * dx + dy * dy);
     double halfWidth  = line->width / 2.0;
-    double halfWidth2 = halfWidth * halfWidth;
     if (!dr) return;
-    
+
+    // Compute the intercepts of line width bounds and determine maximum
+    // bounds in each axis
+
+    xOffset = -halfWidth * dy / dr;
+    yOffset =  halfWidth * dx / dr;
+
+    offsetLine = *line;
+    LineMove (&offsetLine, xOffset, yOffset);
+    if (LineClip (&offsetLine, numCols, numRows))
+        MaxBounds (&offsetLine, &xBegin, &xEnd, &yBegin, &yEnd);
+
+    offsetLine = *line;
+    LineMove (&offsetLine, -xOffset, -yOffset);
+    if (LineClip (&offsetLine, numCols, numRows))
+        MaxBounds (&offsetLine, &xBegin, &xEnd, &yBegin, &yEnd);
+
+    if (xBegin > xEnd) 
+        SwapInt (&xBegin, &xEnd);
+    else
+        ++xEnd;
+    if (yBegin > yEnd)
+        SwapInt (&yBegin, &yEnd);
+    else
+        ++yEnd;
+
     // Step point by point based on the dominate axis
     
@@ -307,25 +378,20 @@
         // Compute the x and y offsets for the line width extent
 
-        xOffset = halfWidth * dy / dr;
-        yOffset = halfWidth * dr / dx;
-        yMid   = y1 + slope * (floor (x1 - xOffset) - x1);
-        xBegin = floor (x1 - xOffset);
-        xEnd   = ceil  (x2 + xOffset) + 1.0;
-
-        for (x = xBegin; x < xEnd; ++x)
-        {
-            yBegin = floor (yMid - yOffset);
-            yEnd   = ceil  (yMid + yOffset) + 1.0;
-            for (y = yBegin; y < yEnd; ++y)
+        yMid = y1 + slope * (xBegin - x1);
+        yOffset = fabs (halfWidth * dr / dx);
+
+        for (x = xBegin; x != xEnd; ++x)
+        {
+            yBegin = (int) floor (yMid - yOffset);
+            yEnd   = (int) ceil  (yMid + yOffset) + 1;
+            for (y = yBegin; y != yEnd; ++y)
             {
-                if (DistanceSquared (line, x, y) <= halfWidth2)
+                if (x >=0 && x < numCols && y >= 0 && y < numRows)
                 {
-                    if (x >=0 && y >= 0) {
-                        pixel = psAlloc (sizeof(PixelPos));
-                        pixel->x = (unsigned int) x;
-                        pixel->y = (unsigned int) y;
-                        psArrayAdd (pixels, 1024, pixel);
-                        psFree (pixel);
-                    }
+                    pixel = psAlloc (sizeof(PixelPos));
+                    pixel->x = (unsigned int) x;
+                    pixel->y = (unsigned int) y;
+                    psArrayAdd (pixels, 1024, pixel);
+                    psFree (pixel);
                 }
             }
@@ -354,27 +420,21 @@
         
         // Compute the x and y offsets for the line width extent
-        
-        xOffset = halfWidth * dr / dy;
-        yOffset = halfWidth * dx / dr;
-
-        xMid   = x1 + slope * (floor (y1 - yOffset) - y1);
-        yBegin = floor (y1 - yOffset);
-        yEnd   = ceil  (y2 + yOffset) + 1.0;
-        
-        for (y = yBegin; y < yEnd; ++y)
-        {
-            xBegin = floor (xMid - xOffset);
-            xEnd   = ceil  (xMid + xOffset) + 1.0;
-            for (x = xBegin; x < xEnd; ++x)
+
+        xMid = x1 + slope * (yBegin - y1);
+        xOffset = fabs (halfWidth * dr / dy);
+        
+        for (y = yBegin; y != yEnd; ++y)
+        {
+            xBegin = (int) floor (xMid - xOffset);
+            xEnd   = (int) ceil  (xMid + xOffset) + 1;
+            for (x = xBegin; x != xEnd; ++x)
             {
-                if (DistanceSquared (line, x, y) <= halfWidth2)
+                if (x >=0 && x < numCols && y >= 0 && y < numRows)
                 {
-                    if (x >=0 && y >= 0) {
-                        pixel = psAlloc (sizeof(PixelPos));
-                        pixel->x = (unsigned int) x;
-                        pixel->y = (unsigned int) y;
-                        psArrayAdd (pixels, 1024, pixel);
-                        psFree (pixel);
-                    }
+                    pixel = psAlloc (sizeof(PixelPos));
+                    pixel->x = (unsigned int) x;
+                    pixel->y = (unsigned int) y;
+                    psArrayAdd (pixels, 1024, pixel);
+                    psFree (pixel);
                 }
             }