Changeset 2655

Timestamp:

Apr 4, 2006, 3:49:42 PM (18 years ago)

Author:

duncan

Message:

Adding ability to remove points outside of the mesh when interpolating

Location:

inundation

Files:

• fit_interpolate/interpolate.py (modified) (8 diffs)
• fit_interpolate/test_interpolate.py (modified) (12 diffs)
• utilities/polygon.py (modified) (2 diffs)
• utilities/test_polygon.py (modified) (1 diff)

inundation/fit_interpolate/interpolate.py

@@ -23 +23 @@
 
 from Numeric import zeros, array, Float, Int, dot, transpose, concatenate, \
-     ArrayType, allclose, take, NewAxis
+     ArrayType, allclose, take, NewAxis, arange
 
 from pyvolution.mesh import Mesh
@@ -31 +31 @@
 from pyvolution.quad import build_quadtree
 from utilities.numerical_tools import ensure_numeric, mean, INF
-from utilities.polygon import inside_polygon
+from utilities.polygon import in_and_outside_polygon
 from geospatial_data.geospatial_data import Geospatial_data
 
@@ -128 +128 @@
         point_coordinates = ensure_numeric(point_coordinates, Float)
         
-        #Remove points falling outside mesh boundary
-        # do this bit later - that sorta means this becomes an object
-        # get a list of what indices are outside the boundary
-        # maybe fill these rows with n/a?
-
+        if verbose: print 'Getting indices inside mesh boundary'
+        self.inside_poly_indices, self.outside_poly_indices  = \
+                     in_and_outside_polygon(point_coordinates,
+                                            self.mesh.get_boundary_polygon(),
+                                            closed = True, verbose = verbose)
         
         #Build n x m interpolation matrix
@@ -143 +143 @@
         A = Sparse(n,m)
 
-        #Compute matrix elements
-        for i in range(n):
+        #Compute matrix elements for points inside the mesh
+        for i in self.inside_poly_indices:
             #For each data_coordinate point
             if verbose and i%((n+10)/10)==0: print 'Doing %d of %d' %(i, n)
@@ -164 +164 @@
                     A[i,j] = sigmas[j]
             else:
-                # Hack so this message only appears when the test is
-                #run in fit_interpoate directory.
-                #Error message reminds me to do something
-                #with points outside the polygon
-                if os.getcwd()[-15:] == "fit_interpolate":
-                    print 'Could not find triangle for point', x 
+                msg = 'Could not find triangle for point', x 
+                raise Exception(msg)
         return A
 
@@ -232 +228 @@
                 self._A_can_be_reused = False
                 start=0
-                z = self.interpolate_block(f, point_coordinates[0:0])
+                # creating a dummy array to concatenate to.
+                z = zeros((0,2)) # 2 for points in 2D space
                 for end in range(start_blocking_len
                                  ,len(point_coordinates)
@@ -251 +248 @@
         doesn't occur.
 
+        Return the point data, z.
+        
         See interpolate for doc info.
         """
@@ -259 +258 @@
 
     def _get_point_data_z(self, f, verbose=False):
-        return self._A * f
-
+        """
+        Return the point data, z.
+        
+        Precondition,
+        The _A matrix has been created
+        """
+        z = self._A * f
+        for i in self.outside_poly_indices: # Taking into account points outside the
+                                            # mesh. 
+            z[i,:] = INF
+        return z
 
 

inundation/fit_interpolate/test_interpolate.py

@@ -20 +20 @@
 from coordinate_transforms.geo_reference import Geo_reference
 from shallow_water import Domain, Transmissive_boundary
-from utilities.numerical_tools import mean
+from utilities.numerical_tools import mean, INF
 from data_manager import get_dataobject
 from geospatial_data.geospatial_data import Geospatial_data
@@ -237 +237 @@
         assert allclose(sum(results, axis=1), 1.0)
 
-        #FIXME - have to change this test to check default info
-    def NO_test_arbitrary_datapoints_some_outside(self):
+    def test_arbitrary_datapoints_some_outside(self):
         """Try arbitrary datapoints one outside the triangle.
         That one should be ignored
@@ -253 +252 @@
         data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44], [5.0, 7.0]]
 
-
-        interp = Interpolate(points, vertices, data, precrop = True)
-        
-        results = interp._build_interpolation_matrix_A(data).todense()
-        assert allclose(sum(results, axis=1), 1.0)
-
-        interp = Interpolate(points, vertices, data, precrop = False)
+        interp = Interpolate(points, vertices, data)
         results = interp._build_interpolation_matrix_A(data).todense()
         assert allclose(sum(results, axis=1), [1,1,1,0])
@@ -301 +294 @@
         assert allclose(A, answer)
 
-
+        #FIXME -  tests are hopefully failing due to points outsdie the
+        # mesh
     def test_interpolate_attributes_to_points(self):
         v0 = [0.0, 0.0]
@@ -442 +436 @@
 
         z = interp.interpolate(f, point_coords)
-        answer = [ [0., 0., 0., 0.]] # (-1,-1)
+        answer = array([ [INF, INF, INF, INF]]) # (-1,-1)
 
         #print "***********"
@@ -452 +446 @@
         # of the mesh? Not currently.  
 
-        assert allclose(z, answer)
-
+        for i in range(4):
+            self.failUnless( z[0,i] == answer[0,i], 'Fail!')
+        
     def test_interpolate_attributes_to_pointsIV(self):
         a = [-1.0, 0.0]
@@ -594 +589 @@
 
     def test_interpolation_interface_time_only(self):
-        """Test spatio-temporal interpolation
-        Test that spatio temporal function performs the correct
-        interpolations in both time and space
-        """
+
+        # Test spatio-temporal interpolation
+        # Test that spatio temporal function performs the correct
+        # interpolations in both time and space
+        
 
 
@@ -667 +663 @@
 
     def test_interpolation_interface_spatial_only(self):
-        """Test spatio-temporal interpolation with constant time
-        """
-
+        # Test spatio-temporal interpolation with constant time
+        
         #Three timesteps
         time = [1.0, 5.0, 6.0]
@@ -728 +723 @@
 
     def test_interpolation_interface(self):
-        """Test spatio-temporal interpolation
-        Test that spatio temporal function performs the correct
-        interpolations in both time and space
-        """
-
+        # Test spatio-temporal interpolation
+        # Test that spatio temporal function performs the correct
+        # interpolations in both time and space
+    
 
         #Three timesteps
@@ -793 +787 @@
 
 
-    def qtest_interpolation_interface(self):
+    def test_points_outside_the_polygon(self):
         a = [-1.0, 0.0]
         b = [3.0, 4.0]
@@ -815 +809 @@
                         [-0.9, -1.5],
                         [0.5, -1.9],
-                        [999999, 9999999]]
+                        [999999, 9999999]] # point Outside poly
         geo_data = Geospatial_data(data_points = point_coords)
 
@@ -827 +821 @@
                   2*linear_function(point_coords) ]
         answer = transpose(answer)
+        answer[11,:] = [INF, INF]
         #print "z",z
-        #print "answer",answer
-        assert allclose(z, answer)
-
+        #print "answer _ fixed",answer
+        assert allclose(z[0:10], answer[0:10])
+        self.failUnless( z[11,1] == answer[11,1], 'Fail!')
+        self.failUnless( z[11,0] == answer[11,0], 'Fail!')
+
+        
+    def test_points_outside_the_polygon(self):
+        a = [-1.0, 0.0]
+        b = [3.0, 4.0]
+        c = [4.0, 1.0]
+        d = [-3.0, 2.0] #3
+        e = [-1.0, -2.0]
+        f = [1.0, -2.0] #5
+
+        vertices = [a, b, c, d,e,f]
+        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
+
+        point_coords = [[-2.0, 2.0],
+                        [-1.0, 1.0],
+                        [9999.0, 9999.0], # point Outside poly
+                        [-9999.0, 1.0], # point Outside poly
+                        [2.0, 1.0],
+                        [0.0, 0.0],
+                        [1.0, 0.0],
+                        [0.0, -1.0],
+                        [-0.2, -0.5],
+                        [-0.9, -1.5],
+                        [0.5, -1.9],
+                        [999999, 9999999]] # point Outside poly
+        geo_data = Geospatial_data(data_points = point_coords)
+
+        interp = Interpolate(vertices, triangles)
+        f = array([linear_function(vertices),2*linear_function(vertices) ])
+        f = transpose(f)
+        #print "f",f
+        z = interp.interpolate(f, geo_data)
+        #z = interp.interpolate(f, point_coords)
+        answer = [linear_function(point_coords),
+                  2*linear_function(point_coords) ]
+        answer = transpose(answer)
+        answer[2,:] = [INF, INF]
+        answer[3,:] = [INF, INF]
+        answer[11,:] = [INF, INF]
+        #print "z",z
+        #print "answer _ fixed",answer
+        assert allclose(z[0:1], answer[0:1])
+        assert allclose(z[4:10], answer[4:10])
+        for i in [2,3,11]:
+            self.failUnless( z[i,1] == answer[11,1], 'Fail!')
+            self.failUnless( z[i,0] == answer[11,0], 'Fail!')
         
 #-------------------------------------------------------------
 if __name__ == "__main__":
+
+    #suite = unittest.makeSuite(Test_Interpolate,'test')
     suite = unittest.makeSuite(Test_Interpolate,'test')
+    #suite = unittest.makeSuite(Test_Interpolate,'test_points_outside_the_polygon')
     runner = unittest.TextTestRunner(verbosity=1)
     runner.run(suite)

inundation/utilities/polygon.py

@@ -127 +127 @@
                                                 verbose=verbose)
 
+    
     if one_point:
         return count != 1
@@ -135 +136 @@
         else:
             return indices[count:][::-1]  #return reversed
+       
+
+def in_and_outside_polygon(points, polygon, closed = True, verbose = False):
+    """Determine points inside and outside a polygon
+
+       See separate_points_by_polygon for documentation
+
+       Returns an array of points inside and an array of points outside the polygon
+    """
+
+    if verbose: print 'Checking input to outside_polygon'
+    try:
+        points = ensure_numeric(points, Float)
+    except NameError, e:
+        raise NameError, e
+    except:
+        msg = 'Points could not be converted to Numeric array'
+        raise msg
+
+    try:
+        polygon = ensure_numeric(polygon, Float)
+    except NameError, e:
+        raise NameError, e
+    except:
+        msg = 'Polygon could not be converted to Numeric array'
+        raise msg
+
+
+
+    if len(points.shape) == 1:
+        one_point = True
+        points = reshape(points, (1,2))
+    else:
+        one_point = False
+
+    indices, count = separate_points_by_polygon(points, polygon,
+                                                closed=closed,
+                                                verbose=verbose)
+    # Returns an array of points inside and an array of points outside
+    # the polygon
+  
+    if count == len(indices):
+        # No points are outside
+        return indices[:count],[]
+    else:
+        return  indices[:count], indices[count:][::-1]  #return reversed
 
 

inundation/utilities/test_polygon.py

@@ -477 +477 @@
         assert inside_polygon(point, polygon)
 
+    def test_in_and_outside_polygon_main(self):
+
+
+        #Simplest case: Polygon is the unit square
+        polygon = [[0,0], [1,0], [1,1], [0,1]]
+
+        inside, outside =  in_and_outside_polygon( (0.5, 0.5), polygon )
+        assert inside[0] == 0
+        assert len(outside) == 0
+        
+        inside, outside =  in_and_outside_polygon(  (1., 0.5), polygon, closed=True)
+        assert inside[0] == 0
+        assert len(outside) == 0
+        
+        inside, outside =  in_and_outside_polygon(  (1., 0.5), polygon, closed=False)
+        assert len(inside) == 0
+        assert outside[0] == 0
+
+        points =  [(1., 0.25),(1., 0.75) ]
+        inside, outside =  in_and_outside_polygon( points, polygon, closed=True)
+        assert (inside, [0,1])
+        assert len(outside) == 0
+        
+        inside, outside =  in_and_outside_polygon( points, polygon, closed=False)
+        assert len(inside) == 0
+        assert (outside, [0,1])
+
+       
+        points =  [(100., 0.25),(0.5, 0.5) ] 
+        inside, outside =  in_and_outside_polygon( points, polygon)
+        assert (inside, [1])
+        assert outside[0] == 0
+        
+        points =  [(100., 0.25),(0.5, 0.5), (39,20), (0.6,0.7),(56,43),(67,90) ] 
+        inside, outside =  in_and_outside_polygon( points, polygon)
+        assert (inside, [1,3])
+        assert (outside, [0,2,4,5])
+        
+    def zzztest_inside_polygon_main(self):  
+        print "inside",inside
+        print "outside",outside
+        
+        assert not inside_polygon( (0.5, 1.5), polygon )
+        assert not inside_polygon( (0.5, -0.5), polygon )
+        assert not inside_polygon( (-0.5, 0.5), polygon )
+        assert not inside_polygon( (1.5, 0.5), polygon )
+
+        #Try point on borders
+        assert inside_polygon( (1., 0.5), polygon, closed=True)
+        assert inside_polygon( (0.5, 1), polygon, closed=True)
+        assert inside_polygon( (0., 0.5), polygon, closed=True)
+        assert inside_polygon( (0.5, 0.), polygon, closed=True)
+
+        assert not inside_polygon( (0.5, 1), polygon, closed=False)
+        assert not inside_polygon( (0., 0.5), polygon, closed=False)
+        assert not inside_polygon( (0.5, 0.), polygon, closed=False)
+        assert not inside_polygon( (1., 0.5), polygon, closed=False)
+
+
+
+        #From real example (that failed)
+        polygon = [[20,20], [40,20], [40,40], [20,40]]
+        points = [ [40, 50] ]
+        res = inside_polygon(points, polygon)
+        assert len(res) == 0
+
+        polygon = [[20,20], [40,20], [40,40], [20,40]]
+        points = [ [25, 25], [30, 20], [40, 50], [90, 20], [40, 90] ]
+        res = inside_polygon(points, polygon)
+        assert len(res) == 2
+        assert allclose(res, [0,1])
+        
 #-------------------------------------------------------------
 if __name__ == "__main__":