Feature perspective

.github/workflows/run-unit-tests.yml

@@ -15,7 +15,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: ['3.9', '3.11', '3.12', '3.x']
+        python-version: ['3.9', '3.11', '3.12', '3.13']
         os: [ macos-latest, ubuntu-latest, windows-latest ]
 
     steps:

README.md

@@ -16,14 +16,11 @@ pip install plan-rect
 
 ## Usage
 
-Rectification is performed with the ``plan-rect`` command.  It requires an image, camera interior parameters and marker locations as inputs, and creates a rectified image and rectification data file as outputs.  Its options are described below:
+Rectification is performed with the ``plan-rect`` command.  It requires an image and marker locations as inputs, and creates a rectified image and rectification data file as outputs.  Its options are described below:
 
 | Option                        | Value                                                          | Description                                                                                                                                        |
 |-------------------------------|----------------------------------------------------------------|----------------------------------------------------------------------------------------------------------------------------------------------------|
 | ``-im`` / ``--image``         | FILE                                                           | Path / URI of the source image (required).                                                                                                         |
-| ``-fl`` / ``--focal-len``     | FLOAT                                                          | Camera focal length (any units).                                                                                                                   |
-| ``-ss`` / ``--sensor-size``   | WIDTH HEIGHT                                                   | Camera sensor size in the same units as ``-fl`` / ``--focal-len``.                                                                                 |
-| ``-ip`` / ``--int-param``     | FILE                                                           | Path / URI of an Orthority format interior parameter file.                                                                                         |
 | ``-m`` / ``--marker``         | ID X Y COL ROW                                                 | Marker ID and location in world and pixel coordinates, with pixel coordinate origin at the bottom left image corner.                               | 
 | ``-g`` / ``--gcp``            | FILE                                                           | Path / URI of an Orthority GCP file defining marker locations.                                                                                     |
 | ``-r`` / ``--res``            | FLOAT                                                          | Rectified pixel size in meters.  Can be used twice for non-square pixels: ``--res WIDTH --res HEIGHT``.  Defaults to the ground sampling distance. |
@@ -35,35 +32,33 @@ Rectification is performed with the ``plan-rect`` command.  It requires an image
 | ``--version``                 |                                                                | Show the version and exit.                                                                                                                         |
 | ``--help``                    |                                                                | Show the help and exit.                                                                                                                            |
 
-Camera interior parameters are required with either ``-fl`` / ``--focal-len`` and ``-ss`` / ``--sensor-size``, or ``-ip`` / ``--int-param``.  
-
-Marker locations are required with either ``-m`` / ``--marker`` or ``-g`` / ``--gcp``.  The ``-m`` / ``--marker`` option can be provided multiple times. At least three markers are required. 
+Marker locations are required with either ``-m`` / ``--marker`` or ``-g`` / ``--gcp``.  The ``-m`` / ``--marker`` option can be provided multiple times. At least four markers are required. 
 
 
 ### Examples
 
-Supply interior parameters with ``-fl`` / ``--focal-len`` and ``-ss`` / ``--sensor-size``, and marker locations with ``-m`` / ``--marker``:
+Supply marker locations with ``-m`` / ``--marker``:
 
 ```commandline
-plan-rect --image source.jpg --focal-len 50 --sensor-size 31.290 23.491 --marker A 0 0 1002 1221 --marker B 2.659 0 4261 1067 --marker C 2.321 5.198 3440 3706 --marker D -0.313 4.729 1410 3663
+plan-rect --image source.jpg --marker A 0 0 1002 1221 --marker B 2.659 0 4261 1067 --marker C 2.321 5.198 3440 3706 --marker D -0.313 4.729 1410 3663
 ```
 
-Supply interior parameters with ``-ip`` / ``--int-param``  and marker locations with ``-g`` / ``--gcp``:
+Supply marker locations with ``-g`` / ``--gcp``:
 
 ```commandline
-plan-rect --image source.jpg --int-param int_param.yaml --gcp gcps.geojson
+plan-rect --image source.jpg --gcp gcps.geojson
 ```
 
 Set the rectified image pixel size with ``-r`` / ``--res``:
 
 ```commandline
-plan-rect --image source.jpg --res 0.01 --int-param int_param.yaml --gcp gcps.geojson
+plan-rect --image source.jpg --res 0.01 --gcp gcps.geojson
 ```
 
-Export interior parameters and marker locations to Orthority format files in the ``data`` directory, overwriting existing outputs:
+Export marker locations to an Orthority GCP file in the ``data`` directory, overwriting any existing file:
 
 ```commandline
-plan-rect --image source.jpg --export-params --out-dir data --overwrite --focal-len 50 --sensor-size 31.290 23.491 --marker A 0 0 1002 1221 --marker B 2.659 0 4261 1067 --marker C 2.321 5.198 3440 3706 --marker D -0.313 4.729 1410 3663
+plan-rect --image source.jpg --export-params --out-dir data --overwrite --marker A 0 0 1002 1221 --marker B 2.659 0 4261 1067 --marker C 2.321 5.198 3440 3706 --marker D -0.313 4.729 1410 3663
 ```
 
 ## Licence

plan_rect/camera.py

@@ -0,0 +1,48 @@
+# Copyright The Plan-Rect Contributors.
+#
+# This file is part of Plan-Rect.
+#
+# Plan-Rect is free software: you can redistribute it and/or modify it under the
+# terms of the GNU Affero General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later version.
+#
+# Plan-Rect is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
+# PARTICULAR PURPOSE.  See the GNU Affero General Public License for more details.
+#
+# You should have received a copy of the GNU Affero General Public License along with
+# Plan-Rect. If not, see <https://www.gnu.org/licenses/>.
+from __future__ import annotations
+
+import numpy as np
+from orthority.camera import Camera
+
+
+class PerspectiveCamera(Camera):
+    """
+    Perspective camera for projecting between planes.
+
+    :param im_size:
+        Image (width, height) in pixels.
+    :param tform:
+        Perspective transform from world to pixel coordinates, as a 3-by-3 array.
+    """
+
+    def __init__(self, im_size: tuple[int, int], tform: np.ndarray):
+        self._im_size = im_size
+        self._tform = tform
+
+    def world_to_pixel(self, xyz: np.ndarray) -> np.ndarray:
+        xyz_ = np.vstack((xyz[:2], np.ones((1, xyz.shape[1]))))
+        ji_ = self._tform.dot(xyz_)
+        ji = ji_[:2] / ji_[2]
+        return ji
+
+    def pixel_to_world_z(self, ji: np.ndarray, z: float | np.ndarray) -> np.ndarray:
+        # only allow projecting to the world z=0 plane
+        assert np.all(z == 0)
+        ji_ = np.vstack((ji, np.ones((1, ji.shape[1]))))
+        xyz = np.linalg.inv(self._tform).dot(ji_)
+        xyz /= xyz[2]
+        xyz[2] = 0
+        return xyz

plan_rect/cli.py

@@ -20,17 +20,17 @@
 from pathlib import Path
 
 import click
+import cv2
 import fsspec
 import numpy as np
 import rasterio as rio
 from fsspec.core import OpenFile
 from orthority import param_io
-from orthority.camera import create_camera
 from orthority.common import OpenRaster, join_ofile
-from orthority.enums import CameraType, Interp
-from orthority.fit import fit_frame_exterior
+from orthority.enums import Interp
 from rasterio.errors import NotGeoreferencedWarning
 
+from plan_rect.camera import PerspectiveCamera
 from plan_rect.param_io import write_rectification_data
 from plan_rect.rectify import rectify
 from plan_rect.version import __version__
@@ -90,30 +90,6 @@ def _dir_cb(ctx: click.Context, param: click.Parameter, uri_path: str) -> OpenFi
     callback=partial(_file_cb, mode='rb'),
     help='Path / URI of the source image.',
 )
-@click.option(
-    '-fl',
-    '--focal-len',
-    type=click.FLOAT,
-    default=None,
-    help='Camera focal length.',
-)
-@click.option(
-    '-ss',
-    '--sensor-size',
-    type=(float, float),
-    metavar='WIDTH HEIGHT',
-    default=None,
-    help='Camera sensor size in the same units as --focal-len.',
-)
-@click.option(
-    '-ip',
-    '--int-param',
-    'int_param_file',
-    type=click.Path(dir_okay=False),
-    default=None,
-    callback=partial(_file_cb, mode='rt'),
-    help='Path / URI of an Orthority format interior parameter file.',
-)
 @click.option(
     '-m',
     '--marker',
@@ -166,10 +142,9 @@ def _dir_cb(ctx: click.Context, param: click.Parameter, uri_path: str) -> OpenFi
     '-ep',
     '--export-params',
     is_flag=True,
-    type=click.BOOL,
     default=False,
     show_default=True,
-    help='Export interior parameters and markers to Orthority format files and exit.',
+    help='Export markers to an Orthority GCP file and exit.',
 )
 @click.option(
     '-od',
@@ -184,7 +159,6 @@ def _dir_cb(ctx: click.Context, param: click.Parameter, uri_path: str) -> OpenFi
     '-o',
     '--overwrite',
     is_flag=True,
-    type=click.BOOL,
     default=False,
     show_default=True,
     help='Overwrite existing output(s).',
@@ -194,9 +168,6 @@ def _dir_cb(ctx: click.Context, param: click.Parameter, uri_path: str) -> OpenFi
 def cli(
     ctx: click.Context,
     image_file: OpenFile,
-    focal_len: float,
-    sensor_size: tuple[float, float],
-    int_param_file: OpenFile,
     markers: tuple[tuple[str, float, float, float, float]],
     gcp_file: OpenFile,
     export_params: bool,
@@ -206,46 +177,24 @@ def cli(
 ):
     """Rectify an image onto a plane.
 
-    Camera interior parameters are required with either '-fl' / '--focal-len' and
-    '-ss' / '--sensor-size', or '-ip' / '--int-param'.
-
     Marker locations are required with either '-m' / '--marker' or '-g' / '--gcp'.
-    The '-m' / '--marker' option can be provided multiple times.  At least three
+    The '-m' / '--marker' option can be provided multiple times.  At least four
     markers are required.
     """
     # silence not georeferenced warnings
     warnings.simplefilter('ignore', category=NotGeoreferencedWarning)
     # enter rasterio environment
     ctx.with_resource(rio.Env(GDAL_NUM_THREADS='ALL_CPUS', GTIFF_FORCE_RGBA=False))
 
-    # form an interior parameter dictionary
-    if (focal_len is None or not sensor_size) and not int_param_file:
-        raise click.UsageError(
-            "Interior parameters are required with '-fl' / '--focal-len' and '-ss' "
-            "/'--sensor-size', or '-ip' / '--int-param'."
-        )
-    if focal_len and sensor_size:
-        with OpenRaster(image_file, 'r') as src_im:
-            im_size = src_im.shape[::-1]
-        int_param = dict(
-            cam_type=CameraType.pinhole,
-            im_size=im_size,
-            focal_len=focal_len,
-            sensor_size=sensor_size,
-        )
-        int_param_dict = {'PR-Camera': int_param}
-    else:
-        int_param_dict = param_io.read_oty_int_param(int_param_file)
-        int_param = next(iter(int_param_dict.values()))
-        im_size = int_param['im_size']
-
     # form a GCP dictionary
     image_path = Path(image_file.path)
     if not markers and not gcp_file:
         raise click.UsageError(
             "Marker locations are required with either '-m' / '--marker' or '-g' / "
             "'--gcp'."
         )
+    with OpenRaster(image_file, 'r') as src_im:
+        im_size = src_im.shape[::-1]
 
     if markers:
         # convert marker locations to GCPs with pixel coordinates converted from BL to
@@ -272,29 +221,28 @@ def cli(
             )
         gcp_dict = {image_path.name: gcps}
 
-    if len(gcps) < 3:
-        raise click.UsageError('At least three markers are required.')
+    if len(gcps) < 4:
+        raise click.UsageError('At least four markers are required.')
 
     if export_params:
         # export interior parameters and GCPs as orthority format files, and exit
         gcp_file = join_ofile(out_dir, 'gcps.geojson', mode='wt')
         param_io.write_gcps(gcp_file, gcp_dict, overwrite=overwrite)
-        int_param_file = join_ofile(out_dir, 'int_param.yaml', mode='wt')
-        param_io.write_int_param(int_param_file, int_param_dict, overwrite=overwrite)
-        click.echo(f"Orthority format files written to: '{out_dir.path}'.")
+        click.echo(f"Orthority GCP file written to: '{gcp_file.path}'.")
         return
 
-    # fit exterior parameters & create camera
-    ext_param_dict = fit_frame_exterior(int_param_dict, gcp_dict)
-    ext_param = next(iter(ext_param_dict.values()))
-    cam = create_camera(**int_param, xyz=ext_param['xyz'], opk=ext_param['opk'])
+    # fit perspective transform & create camera
+    gcp_ji = np.array([gcp['ji'] for gcp in gcps]).T
+    gcp_xyz = np.array([gcp['xyz'] for gcp in gcps]).T
+    tform, _ = cv2.findHomography(
+        gcp_xyz[:2].T.astype('float32'), gcp_ji.T.astype('float32')
+    )
+    cam = PerspectiveCamera(im_size, tform)
 
     # find & print fitting error
-    gcp_xyz = np.array([gcp['xyz'] for gcp in gcps]).T
-    gcp_ji = np.array([gcp['ji'] for gcp in gcps]).T
     cam_ji = cam.world_to_pixel(gcp_xyz)
     err = np.sqrt(np.sum((gcp_ji - cam_ji) ** 2, axis=0)).mean()
-    click.echo(f'RMS fitting error: {err:.4f} (pixels).')
+    click.echo(f'RMS fitting error: {err:.4e} (pixels).')
 
     # rectify
     rect_array, transform = rectify(image_file, cam, **kwargs)
@@ -315,9 +263,8 @@ def cli(
     # find pixel coordinates of markers in the rectified image (0.5 translation
     # shifts the transform from the GDAL / Rasterio convention to give integer
     # coordinates that refer to pixel centers)
-    cam_xyz = cam.pixel_to_world_z(gcp_ji, z=0)
     inv_transform = ~(rio.Affine(*transform) * rio.Affine.translation(0.5, 0.5))
-    rect_ji = np.array(inv_transform * cam_xyz[:2])
+    rect_ji = np.array(inv_transform * gcp_xyz[:2])
 
     # create a rectified marker list, converting pixel coordinates from TL to BL
     # origin convention
@@ -330,11 +277,7 @@ def cli(
     # write rectification data
     rect_data_file = join_ofile(out_dir, 'pixeldata.txt', mode='wt')
     write_rectification_data(
-        rect_data_file,
-        image_path.name,
-        int_param,
-        rect_markers,
-        overwrite=overwrite,
+        rect_data_file, image_path.name, im_size, rect_markers, overwrite=overwrite
     )
     click.echo(f"Output files written to: '{out_dir.path}'.")
 

plan_rect/param_io.py

@@ -24,7 +24,7 @@
 def write_rectification_data(
     file: str | PathLike | OpenFile,
     src_name: str,
-    int_param: dict[str, Any],
+    im_size: tuple[int, int],
     markers: list[dict[str, Any]],
     overwrite: bool = False,
 ):
@@ -37,8 +37,8 @@ def write_rectification_data(
         mode (``'wt'``).
     :param src_name:
         Source image file name.
-    :param int_param:
-        Camera interior parameters.
+    :param im_size:
+        Image (width, height) in pixels.
     :param markers:
         Markers as a list of dictionaries with ``id``: <marker name>, and ``ji``:
         <pixel coordinate> items.
@@ -47,20 +47,7 @@ def write_rectification_data(
     """
     with Open(file, 'wt', overwrite=overwrite) as f:
         f.write(f'Photo;{src_name}\n')
-        im_size = int_param['im_size']
         f.write(f'Size;{im_size[0]},{im_size[1]};px\n')
 
-        focal_len = int_param['focal_len']
-        focal_str = (
-            f'{focal_len:.4f}'
-            if not isinstance(focal_len, tuple)
-            else f'{focal_len[0]:.4f},{focal_len[1]:.4f}'
-        )
-        f.write(f'Lens;{focal_str};?\n')
-
-        if 'sensor_size' in int_param:
-            sensor_size = int_param['sensor_size']
-            f.write(f'Sensor;{sensor_size[0]:.4f},{sensor_size[1]:.4f};?\n')
-
         for m in markers:
             f.write(f'{m["id"]};{m["ji"][0]:.4f},{m["ji"][1]:.4f},0\n')

plan_rect/rectify.py

@@ -20,7 +20,7 @@
 import numpy as np
 import rasterio as rio
 from fsspec.core import OpenFile
-from orthority.camera import FrameCamera
+from orthority.camera import Camera
 from orthority.common import OpenRaster
 from orthority.enums import Interp
 
@@ -38,7 +38,7 @@ def _area_poly(coords: np.ndarray) -> float:
 
 def rectify(
     src_file: str | PathLike | OpenFile,
-    camera: FrameCamera,
+    camera: Camera,
     resolution: tuple[float, float] | None = None,
     interp: str | Interp = Interp.cubic,
     nodata: int | float | None = None,