Merge bf0ce25 into 6a770ac

docs/usage.rst

@@ -12,19 +12,25 @@ In order to load a dataset (a single collection of images)::
     nadirs, obliques = evtech.load_dataset('/path/to/dataset')
 
     # Get the image data for the first nadir camera as a numpy array
-    img = nadirs[0].load_image()
-
-From where you can also look at operations with the camera, such as projecting a ray from the camera at a given pixel. Also projecting a latitude, longitude, elevation point into the camera to get the pixel location::
-
-    # Fill in 
+    nadir_cam = nadirs[0]
+    img = nadir_cam.load_image()
 
 The camera object also stores the average elevation of the image and the geographical bounds of the image, which can be retreived as a Shapely polgon::
 
     # Fill in
 
-We have provided a simple single-image height measurement function that uses the camera and two points to compute the height of an object::
+From here you can also look at operations with the camera, such as projecting a ray from the camera at a given pixel. Also projecting a latitude, longitude, elevation point into the camera to get the pixel location::
 
-    # Fill in
+    image_pt = nadir_cam.project_to_camera(lon, lat, elevation)
+    ray = nadir_cam.project_from_camera(co, row)
+
+Rays can be used with ground elevation to find the intersection of the pixel and the ground::
+
+    ground_point_at_pxiel = ray.intersect_at_elevation(nadir_cam.elevation)
+
+We have provided a simple single-image height measurement function that uses the camera and two points to compute the height of an object at a given elevation::
+
+    height = nadir_cam.height_between_points(base_img_pt, peak_image_pt, nadir_cam.elevation)
 
 Lastly, we have proivded a function to take multiple cameras and associated image points, and triangulate a three dimensional point::
 

evtech/camera.py

@@ -3,6 +3,7 @@
 import numpy as np
 import utm
 import cv2
+import math
 
 from pyproj import CRS, Transformer
 
@@ -22,7 +23,7 @@ class Camera():
     :type geo_bounds: list
     :param elev: The average elevation of the image
     :type elev: float
-    :param crs: The coordinate system for the projection matrix
+    :param crs: The coordinate system for the projection matrix (Must be linear such as UTM)
     :type crs: class: `pyproj.CRS`
     :param image_path: The filepath to the image data
     :type path: str
@@ -62,6 +63,9 @@ def project_from_camera(self, col, row):
         # subset the projection matrix
         m = self.projection_matrix[0:3,0:3]
 
+        # Offset for crop
+        col, row = self.to_full_image(col,row)
+
         # Get point and project to to normalized plane
         pt = np.transpose(np.array([[col,row,1.0]]))
         norm_pt = np.linalg.inv(m) @ pt
@@ -97,7 +101,63 @@ def project_to_camera(self, lon, lat, elevation):
         img_pt[1] -= self.image_bounds[1]
         img_pt = np.transpose(img_pt)
         return img_pt[0][0:2]
-
+
+    def height_between_points(self, base_point, peak_point, elev=None):
+        """ Compute the height between two image points, given the elevation of the base point. 
+        If no elevation is passed the stored elevation will be used.
+
+        :param base_point: The image point at the given elevation
+        :type base_point: list
+        :param peak_point: The image point to compute the height at
+        :type peak_point: list
+        :param elev: The associated elevation of the base point, defaults to None
+        :type elev: float, optional
+        """
+        if not elev:
+            elev = self.elevation
+
+        # Compute rays from each point
+        base_ray = self.project_from_camera(base_point[0],base_point[1])
+        peak_ray = self.project_from_camera(peak_point[0],peak_point[1])
+
+        # Compute the cosine of the angle between the two rays
+        [base_dir] = np.transpose(base_ray.direction)
+        [peak_dir] = np.transpose(peak_ray.direction)
+
+        dt = np.dot(base_dir, peak_dir)
+        c = dt/np.linalg.norm(base_dir)/np.linalg.norm(peak_dir)
+
+        # Compute depth at given elevation
+        depth = base_ray.depth_at_elevation(elev)
+
+        # Get depth at midpoint between points
+        depth_mid = depth*c
+
+        # Extract focal length from proj matrix
+        camera_matrix,_,_,_,_,_,_ = cv2.decomposeProjectionMatrix(self.projection_matrix)
+        focal_x = camera_matrix[0][0]
+        focal_y = camera_matrix[1][1]
+        focal = (focal_x + focal_y) / 2
+
+        # Compute height using simlar triangles
+        dist = np.linalg.norm(np.array(base_point) - np.array(peak_point))
+        height = dist / focal * depth_mid
+        return height[0]
+
+    def to_full_image(self, col, row):
+        """ Convert an image point from the subset image to the full image
+        
+        :param col: The column to offset
+        :type col: float
+        :param row: The row to offset
+        :type row: float
+        :return: The offset row,col
+        :rtype: tuple
+        """
+        r_col = col + self.image_bounds[1]
+        r_row = row + self.image_bounds[0]
+        return r_row, r_col
+
     def load_image(self, loader=cv2.imread):
         """ Load the image for this camera
         
@@ -116,12 +176,56 @@ def camera_from_json(json_data, image_path = ""):
     :param image_path: The path to the associated image data
     :type image_path: str, optional
     :return: A camera object
-    :rtype: classs
+    :rtype: evtech.Camera
     """
     # Determine proper UTM zone
     crs = utm_crs_from_latlon(json_data["geo_bounds"][1], 
                                 json_data["geo_bounds"][0])
 
     return Camera(np.array(json_data["projection"]), json_data["bounds"], 
                     json_data["camera_center"],json_data["geo_bounds"], 
-                    json_data["elevation"], crs, image_path)
+                    json_data["elevation"], crs, image_path)
+
+# def triangulate_point_from_cameras(cameras, points, to_latlng = False):
+# #     """ Triangulates a 3D point from two cameras and two image points
+
+# #     :param camera_1: The first camera
+# #     :type camera_1: evtech.Camera
+# #     :param point_1: The first image point as a 2 element list [col,row]
+# #     :type point_1: list
+# #     :param camera_2: The second camera
+# #     :type camera_2: evtech.Camera
+# #     :param point_2: The second image point as a 2 element list [col,row]
+# #     :type point_2: list
+# #     :param to_latlng: Flag to return the point as lat/lng/elevation, otherwise will return in the camera's CRS
+# #     :type to_latlng: bool
+# #     :return: A 3d point
+# #     :rtype: numpy.Array
+# #     """
+
+
+#     weights = [1.0 for i in range(len(cameras))]
+#     for i in range(10):
+#         A = []
+#         for cam, pt, weight in zip(cameras, points, weights):
+#             x,y = cam.to_full_image(float(pt[0]),float(pt[1]))
+#             A.append((x * cam.projection_matrix[2,:] - cam.projection_matrix[0,:])/weight)
+#             A.append((y * cam.projection_matrix[2,:] - cam.projection_matrix[1,:])/weight)
+
+#         # Solve for X
+#         u,d,vt=np.linalg.svd(A)
+#         X = vt[-1,0:3]/vt[-1,3] # normalize
+#         wX = np.transpose(np.append(X, [1.0]))
+
+#         # update weights
+#         for idx, (cam, weight) in enumerate(zip(cameras, weights)):
+#             weights[idx] = cam.projection_matrix[2,:] @ wX
+
+#         print(weights)
+
+#     # Convert if needed
+#     if to_latlng:
+#         transformer = Transformer.from_crs(cameras[0].crs, CRS.from_user_input(4326), always_xy=True)
+#         X[0],X[1],X[2] = transformer.transform(X[0], X[1], X[2])
+
+#     return X

evtech/ray.py

@@ -2,6 +2,7 @@
 
 import numpy as np
 from sklearn import preprocessing
+from pyproj import CRS, Transformer
 
 class Ray():
     """ A class to represent rays in three dimensional space
@@ -19,4 +20,43 @@ def __init__(self, origin, direction, crs):
         """
         self.origin = np.transpose(np.array([origin]))
         self.direction = np.transpose(preprocessing.normalize(np.array([direction]), norm='l2'))
-        self.crs = crs
+        self.crs = crs
+
+    def point_at_depth(self, depth):
+        """ Return a 3D point at a given depth along the ray
+        
+        :param depth: The depth along the ray
+        :type depth: float
+        :return: A 3d point
+        :rtype: numpy.array
+        """
+        return self.origin + depth * self.direction
+
+    def depth_at_elevation(self, elevation):
+        """ Return the depth at a given elevation
+        
+        :param elevation: The elevation to get the depth of
+        :type elevation: float
+        """
+        return (elevation - self.origin[2])/self.direction[2]
+
+    def intersect_at_elevation(self, elevation, latlng=True):
+        """Return a three dimensional point intersected at a given elevation
+        
+        :param elevation: The elevation to intersect
+        :type elevation: float
+        :param latlng: Return the point as lat,lng, elevation, defaults to True
+        :type latlng: bool, optional
+        :return: A 3d point
+        :rtype: numpy.array
+        """
+
+        depth = self.depth_at_elevation(elevation)
+        pt = self.point_at_depth(depth)
+
+        if latlng:
+            transformer = Transformer.from_crs(self.crs, CRS.from_user_input(4326), always_xy=True)
+            pt[0],pt[1],pt[2] = transformer.transform(pt[0], pt[1], pt[2])
+
+        [pt] = np.transpose(pt)
+        return pt

tests/test_camera.py

@@ -9,6 +9,7 @@
 from pyproj import CRS
 from evtech import Camera
 from evtech import camera_from_json
+#from evtech import triangulate_point_from_cameras
 
 class TestCamera(unittest.TestCase):
     """Tests for `evtech.camera` package."""
@@ -48,12 +49,16 @@ def test_project_to_camera(self):
         self.assertTrue(pt[1] >= y_idx and pt[1] < y_idx+1)
 
     def test_project_from_camera(self):
-
         ray = self.cam.project_from_camera(0,0)
         self.assertEqual(ray.origin[0], self.cen[0])
         self.assertEqual(ray.origin[1], self.cen[1])
         self.assertEqual(ray.origin[2], self.cen[2])
 
+    def test_to_full_image(self):
+        x, y = self.cam.to_full_image(0,0)
+        self.assertEqual(x, self.bounds[0])
+        self.assertEqual(y, self.bounds[1])
+
     def test_fromjson(self):
         data = {
             "id": "18274434", 
@@ -76,3 +81,82 @@ def test_loadimage(self):
         def loader(img):
             return self.path == img
         self.assertTrue(self.cam.load_image(loader))
+
+    def test_hightbetweenpoints(self):
+        cam1_json = {
+            "id": "13470217", 
+            "shot_id": "13470217_E_ILZLAK020024NeighObliq7313E_150610.jpg", 
+            "warehouse_id": 26946, 
+            "projection": [[1525.5867281279347, -15512.91424561646, -2311.8378111550846, 72183965325.08594], 
+            [-7573.84425712711, 803.6272922226443, -13570.519962708786, -650749980.3668021], 
+            [0.7925646349229568, -0.045523902505363464, -0.608173942069206, -109441.04682805175]], 
+            "bounds": [125, 267, 966, 550], 
+            "camera_center": [408968.8416940464, 4693116.473847266, 1716.97110001749], 
+            "geo_bounds": [-88.07612165733431, 42.38789365082783, -88.07494134030249, 42.389211554600976], 
+            "elevation": 254.16879272460938, 
+            "orientation": "E"
+        }
+        cam1 = camera_from_json(cam1_json)
+        base_pt = [41,118]
+        peak_pt = [35,90]
+
+        height = cam1.height_between_points(base_pt, peak_pt)
+        self.assertAlmostEqual(height, 5.51879092088098)
+
+        height = cam1.height_between_points(base_pt, peak_pt, cam1.elevation)
+        self.assertAlmostEqual(height, 5.51879092088098)
+
+
+
+    # def test_triangualte(self):
+    #     cam1_json = {
+    #         "id": "13470217", 
+    #         "shot_id": "13470217_E_ILZLAK020024NeighObliq7313E_150610.jpg", 
+    #         "warehouse_id": 26946, 
+    #         "projection": [[1525.5867281279347, -15512.91424561646, -2311.8378111550846, 72183965325.08594], 
+    #         [-7573.84425712711, 803.6272922226443, -13570.519962708786, -650749980.3668021], 
+    #         [0.7925646349229568, -0.045523902505363464, -0.608173942069206, -109441.04682805175]], 
+    #         "bounds": [125, 267, 966, 550], 
+    #         "camera_center": [408968.8416940464, 4693116.473847266, 1716.97110001749], 
+    #         "geo_bounds": [-88.07612165733431, 42.38789365082783, -88.07494134030249, 42.389211554600976], 
+    #         "elevation": 254.16879272460938, 
+    #         "orientation": "E"
+    #     }
+
+    #     cam2_json = {
+    #             "id": "13470217", 
+    #             "shot_id": "13470217_N_ILZLAK020024NeighObliq4265N_150610.jpg", 
+    #             "warehouse_id": 26946,
+    #             "projection": [[15116.757193147516, 3196.5054851458067, -2909.892827161719, -21213711732.266273], 
+    #             [-333.16135284503827, -8152.243893734243, -13224.287669937909, 38408498248.73735], 
+    #             [-0.06777729452508616, 0.7652991191741246, -0.6401701362908264, -3561739.5617974782]], 
+    #             "bounds": [3569, 2298, 4299, 3029], 
+    #             "camera_center": [411524.2286809936, 4691886.086275864, 1663.19605194418], 
+    #             "geo_bounds": [-88.07612165733431, 42.38789365082783, -88.07494134030249, 42.389211554600976], 
+    #             "elevation": 254.16879272460938, 
+    #             "orientation": "N"
+    #         }
+
+    #     cam3_json = {
+    #         "id": "13470217", 
+    #         "warehouse_id": "26946", 
+    #         "projection": [[-234.48497951320869, -11689.146112537686, -3420.9549093694854, 54967162069.77626], 
+    #         [-11527.74509904331, 527.9966478964207, -3108.9307732776556, 2267432568.205459], 
+    #         [0.07731721986909759, 0.01342309733163904, -0.996916676327768, -93150.24955090503]],
+    #          "bounds": [4405, 655, 5587, 1420], 
+    #          "camera_center": [411228.51669897616, 4693677.177776167, 1653.5802147550032], 
+    #          "geo_bounds": [-88.07607063663191, 42.387928513288855, -88.07499236028416, 42.38917669615173], 
+    #          "elevation": 250.522
+    #     }
+
+    #     cam1 = camera_from_json(cam1_json)
+    #     cam2 = camera_from_json(cam2_json)
+    #     cam3 = camera_from_json(cam3_json)
+    #     pt1 = [605,171]
+    #     pt2 = [304,536]
+    #     pt3 = [879,441]
+    #     cams = [cam1, cam2]
+    #     pts = [pt1, pt2]
+
+    #     world_pt = triangulate_point_from_cameras(cams, pts, True)
+    #     print(world_pt)

tests/test_ray.py

@@ -7,14 +7,28 @@
 import numpy as np
 import math
 
+from pyproj import CRS
 from evtech import Ray
+from evtech import Camera
 
 class TestRay(unittest.TestCase):
     """Tests for `evtech.ray` package."""
 
     def setUp(self):
+        self.proj = np.array([[-234.48497951320869, -11689.146112537686, -3420.9549093694854, 54967162069.77626], 
+             [-11527.74509904331, 527.9966478964207, -3108.9307732776556, 2267432568.205459], 
+             [0.07731721986909759, 0.01342309733163904, -0.996916676327768, -93150.24955090503]
+             ])
+        self.bounds = [4405, 655, 5587, 1420]
+        self.cen = [411228.51669897616, 4693677.177776167, 1653.5802147550032]
+        self.geo_bounds = [-88.07607063663191, 42.387928513288855, -88.07499236028416, 42.38917669615173]
+        self.elev = 250.522
+        self.crs = CRS.from_user_input(32616)
+        self.path = "foo.jpg"
+        self.cam = Camera(self.proj, self.bounds, self.cen, self.geo_bounds, self.elev, self.crs, self.path)
         pass
 
+
     def test_construct(self):
         ray = Ray([0,0,0],[1,1,1], None)
 
@@ -27,4 +41,7 @@ def test_construct(self):
         self.assertEqual(ray.direction[1],1.0/math.sqrt(3))
         self.assertEqual(ray.direction[2],1.0/math.sqrt(3))
 
-
+    def test_elevation_intersect(self):
+        ray = self.cam.project_from_camera(0,0)
+        pt = ray.intersect_at_elevation(self.elev)
+        self.assertAlmostEqual(pt[2], self.elev)