--- python/plugins/fTools/tools/doGeometry.py.old	Thu Feb 26 20:12:00 2015
+++ python/plugins/fTools/tools/doGeometry.py	Thu Jul 16 09:31:20 2015
@@ -36,6 +36,11 @@
 import voronoi
 from sets import Set
 
+# For comparing equality in lat/long - 10^6 meters are < 1 degree
+# so this gives mm resolution at lat=0. For finer resolution tasks
+# a localised projection will likely be used, so this should be OK.
+FLT_EPSILON = 10**-9
+
 class GeometryDialog( QDialog, Ui_Dialog ):
   def __init__( self, iface, function ):
     QDialog.__init__( self, iface.mainWindow() )
@@ -729,41 +734,91 @@
     """ Clip voronoi function based on code written for Inkscape
         Copyright (C) 2010 Alvin Penner, penner@vaxxine.com
     """
+    def float_eq( a, b ):
+      # the initial == test is for efficiency
+      return a == b or abs( a - b ) <= FLT_EPSILON
+    
     def clip_line( x1, y1, x2, y2, w, h, x, y ):
+      # where both x's are left of the clipping region no line needed
       if x1 < 0 - x and x2 < 0 - x:
         return [ 0, 0, 0, 0 ]
+      # where both x's are right of the clipping region no line needed
       if x1 > w + x and x2 > w + x:
         return [ 0, 0, 0, 0 ]
-      if x1 < 0 - x:
-        y1 = ( y1 * x2 - y2 * x1 ) / ( x2 - x1 )
-        x1 = 0 - x
-      if x2 < 0 - x:
-        y2 = ( y1 * x2 - y2 * x1 ) / ( x2 - x1 )
-        x2 = 0 - x
-      if x1 > w + x:
-        y1 = y1 + ( w + x - x1 ) * ( y2 - y1 ) / ( x2 - x1 )
-        x1 = w + x
-      if x2 > w + x:
-        y2 = y1 + ( w + x - x1 ) *( y2 - y1 ) / ( x2 - x1 )
-        x2 = w + x
+      # where both y's are below the clipping region no line needed
       if y1 < 0 - y and y2 < 0 - y:
         return [ 0, 0, 0, 0 ]
+      # where both y's are above the clipping region no line needed
       if y1 > h + y and y2 > h + y:
         return [ 0, 0, 0, 0 ]
-      if x1 == x2 and y1 == y2:
+      # where end points are too close to tell apart no line needed
+      if float_eq( x1, x2 ) and float_eq( y1, y2 ):
         return [ 0, 0, 0, 0 ]
-      if y1 < 0 - y:
-        x1 = ( x1 * y2 - x2 * y1 ) / ( y2 - y1 )
+      # y = mx + b: calculate gradient but make sure is not a vertical line
+      if not float_eq( x1, x2 ):
+        m = (y2 - y1) / (x2 - x1)
+      else:
+        m = float("+inf") # vertical lines have infinite gradient
+      # y = mx + b: calculate intercept, if horizontal m==0 so use first Y value
+      if ( not float_eq( y1, y2 ) and not float_eq( x1, x2 ) ) or float_eq( y1, y2):
+        b = y1 - m*x1 # b = y1 when y1 == y2
+      else: # y1 != y2 and x1 == x2
+        b = float("nan") # vertical lines have no intercept
+      # for simplicity later define if the line is vertical
+      is_vertical = ( m == float("+inf") and b == float("nan") )
+      # where either of the x1,y1 or x2,y2 points are outside of the clipping region
+      # need to move them to be on the border of the region; this is done by recalculating
+      # the x or y value based on the equation of the line holding the other axis as the
+      # value of the clipping rectangle border that is violated.
+      if x1 < 0 - x: # first point is left of region
+        # impossible to have vertical line here, could be horizontal
+        # but that doesn't matter as m = 0, b = y1 in that case
+        x1 = 0 - x
+        y1 = m*x1 + b
+      if x2 < 0 - x: # second point is left of region
+        # impossible to have vertical line here, could be horizontal
+        # but that doesn't matter as m = 0, b = y1 in that case
+        x2 = 0 - x
+        y2 = m*x2 + b
+      if x1 > w + x: # first point is right of region
+        # impossible to have vertical line here, could be horizontal
+        # but that doesn't matter as m = 0, b = y1 in that case
+        x1 = w + x
+        y1 = m*x1 + b
+      if x2 > w + x: # second point is right of region
+        # impossible to have vertical line here, could be horizontal
+        # but that doesn't matter as m = 0, b = y1 in that case
+        x2 = w + x
+        y2 = m*x2 + b
+      if y1 < 0 - y: # first point is below the region
+        # impossible to have horizontal line here, could be vertical though
+        # if is vertical then x value does not need to change
         y1 = 0 - y
-      if y2 < 0 - y:
-        x2 = ( x1 * y2 - x2 * y1 ) / ( y2 - y1 )
+        if not is_vertical: # if veritical then X doesn't need to change
+          x1 = (y1 - b) / m
+      if y2 < 0 - y: # second point is below the region
+        # impossible to have horizontal line here, could be vertical though
+        # if is vertical then x value does not need to change
         y2 = 0 - y
-      if y1 > h + y:
-        x1 = x1 + ( h + y - y1 ) * ( x2 - x1 ) / ( y2 - y1 )
+        if not is_vertical: # if vertical then X doesn't need to change
+          x2 = (y2 - b) / m
+      if y1 > h + y: # first point is above the region
+        # impossible to have horizontal line here, could be vertical though
+        # if is vertical then x value does not need to change
         y1 = h + y
-      if y2 > h + y:
-        x2 = x1 + ( h + y - y1) * ( x2 - x1 ) / ( y2 - y1 )
+        if not is_vertical: # if vertical then X doesn't need to change
+          x1 = (y1 - b) / m
+      if y2 > h + y: # second point is above the region
+        # impossible to have horizontal line here, could be vertical though
+        # if is vertical then x value does not need to change
         y2 = h + y
+        if not is_vertical: # if vertical then X doesn't need to change
+          x2 = (y2 - b) / m
+      # recheck point approximate equivalance because both points may have been
+      # outside the region and thus were both clipped to a point on the boundary
+      if float_eq( x1, x2 ) and float_eq( y1, y2 ):
+        return [ 0, 0, 0, 0 ]
+      # return the clipped region
       return [ x1, y1, x2, y2 ]
     lines = []
     hasXMin = False
@@ -772,35 +827,35 @@
     hasYMax = False
     for edge in edges:
       if edge[ 1 ] >= 0 and edge[ 2 ] >= 0:       # two vertices
-          [ x1, y1, x2, y2 ] = clip_line( c.vertices[ edge[ 1 ] ][ 0 ], c.vertices[ edge[ 1 ] ][ 1 ], c.vertices[ edge[ 2 ] ][ 0 ], c.vertices[ edge[ 2 ] ][ 1 ], width, height, exX, exY )
+        [ x1, y1, x2, y2 ] = clip_line( c.vertices[ edge[ 1 ] ][ 0 ], c.vertices[ edge[ 1 ] ][ 1 ], c.vertices[ edge[ 2 ] ][ 0 ], c.vertices[ edge[ 2 ] ][ 1 ], width, height, exX, exY )
       elif edge[ 1 ] >= 0:                      # only one vertex
         if c.lines[ edge[ 0 ] ][ 1 ] == 0:      # vertical line
-          xtemp = c.lines[ edge[ 0 ] ][ 2 ] / c.lines[ edge[ 0 ] ][ 0 ]
+          xtemp = c.vertices[ edge[ 1 ] ][ 0 ]
           if c.vertices[ edge[ 1 ] ][ 1 ] > ( height + exY ) / 2:
             ytemp = height + exY
           else:
-            ytemp = 0 - exX
+            ytemp = 0.0 - exY
         else:
           xtemp = width + exX
-          ytemp = ( c.lines[ edge[ 0 ] ][ 2 ] - ( width + exX ) * c.lines[ edge[ 0 ] ][ 0 ] ) / c.lines[ edge[ 0 ] ][ 1 ]
+          ytemp = ( c.lines[ edge[ 0 ] ][ 2 ] - xtemp * c.lines[ edge[ 0 ] ][ 0 ] ) / c.lines[ edge[ 0 ] ][ 1 ]
         [ x1, y1, x2, y2 ] = clip_line( c.vertices[ edge[ 1 ] ][ 0 ], c.vertices[ edge[ 1 ] ][ 1 ], xtemp, ytemp, width, height, exX, exY )
       elif edge[ 2 ] >= 0:                       # only one vertex
         if c.lines[ edge[ 0 ] ][ 1 ] == 0:       # vertical line
-          xtemp = c.lines[ edge[ 0 ] ][ 2 ] / c.lines[ edge[ 0 ] ][ 0 ]
+          xtemp = c.vertices[ edge[ 2 ] ][ 0 ]
           if c.vertices[ edge[ 2 ] ][ 1 ] > ( height + exY ) / 2:
             ytemp = height + exY
           else:
             ytemp = 0.0 - exY
         else:
           xtemp = 0.0 - exX
-          ytemp = c.lines[ edge[ 0 ] ][ 2 ] / c.lines[ edge[ 0 ] ][ 1 ]
+          ytemp = ( c.lines[ edge[ 0 ] ][ 2 ] - xtemp * c.lines[ edge[ 0 ] ][ 0 ] ) / c.lines[ edge[ 0 ] ][ 1 ]
         [ x1, y1, x2, y2 ] = clip_line( xtemp, ytemp, c.vertices[ edge[ 2 ] ][ 0 ], c.vertices[ edge[ 2 ] ][ 1 ], width, height, exX, exY )
       if x1 or x2 or y1 or y2:
         lines.append( QgsPoint( x1 + extent.xMinimum(), y1 + extent.yMinimum() ) )
         lines.append( QgsPoint( x2 + extent.xMinimum(), y2 + extent.yMinimum() ) )
-        if 0 - exX in ( x1, x2 ):
+        if 0.0 - exX in ( x1, x2 ):
           hasXMin = True
-        if 0 - exY in ( y1, y2 ):
+        if 0.0 - exY in ( y1, y2 ):
           hasYMin = True
         if height + exY in ( y1, y2 ):
           hasYMax = True
