Python code for performing canadian grid calculation

canadiangrid.py

import math
from dataclasses import dataclass
from typing import List

__author__ = "Brian Schlining"
__copyright__ = "Copyright 2019, Monterey Bay Aquarium Research Institute"


@dataclass
class Camera:
    """
    Calculates view and image parameters for a camera mounted above a view plane. Here's a few
    important terms:
      1. View plane - The plane that we're trying to image (e.g. the seafloor)
      2. Image plane - The plane of the image (abstract concept)
      3. Principal point - The center of the image on our image plane. (e.g. if you have an image
         400x300 pixels then pixel at coordinates 200, 150 is the principal point). A line can
         be extended from this point to intersect with our view plane.

    References: [1] W. Wakefield and A. Genin, “The use of a Canadian (perspective) grid in deep-sea
      photography,” Deep Sea Research Part A. Oceanographic …, vol. 34, no. 3, pp. 469–478, 1987.

    height -- The height of the camera above the plane
    alpha -- The vertical view angle of the camera in radians
    beta -- The horizontal view angle of the camera in radians
    theta -- The tilt of the camera in radians (0 degree is parallel to the plane, 90 is pointed
                 directly at plane)
    """

    height: float
    alpha: float
    beta: float
    theta: float

    def planeDistance(self) -> float:
        """The distance from the point directly under the camera on the view plane to the principal point
        on the view plane."""
        return self.height / math.tan(self.theta)

    def lensDistance(self) -> float:
        """The distance from the camera to the principal point on the view plane"""
        return self.height / math.sin(self.theta)

    def nearViewEdgeDistance(self) -> float:
        """The distance from the point directly under the camera to the nearest edge of the image on the
        view plane"""
        return self.height / math.tan(self.theta + self.alpha / 2.0)

    def farViewEdgeDistance(self) -> float:
        """The distance from the point directly under the camera to the farthest edge of the image on the
        view plane"""
        return self.height / math.tan(self.theta - self.alpha / 2.0)

    def viewHeight(self) -> float:
        """The distance between the near and far edges of our image projected onto the view plane"""
        return self.farViewEdgeDistance() - self.nearViewEdgeDistance()

    def viewWidth(self) -> float:
        """The width of the image on the view plane at the principal point"""
        return self.lensDistance() * math.tan(self.beta / 2.0) * 2.0

    @staticmethod
    def fromDegrees(height: float, alpha: float, beta: float, theta: float):
        return Camera(height, math.radians(alpha), math.radians(beta), math.radians(theta))


@dataclass
class Pixel:
    """
    Given a pixel within an image and a Camera object, Pixel
    calculates the actual position of the pixel relative to the camera's location. This class uses
    the image coordinate system which has the origin at the top-left corner of the image and +X is
    right and +Y is down.

    camera -- The camera whose capturing images
    width -- The width of the image in pixels
    height -- The height of the image in pixels
    x -- The x coordinate of the pixel of interest (0 index)
    y -- The y coordinate of the pixel of interest (0 index)
    """
    camera: Camera
    width: int
    height: int
    x: int
    y: int

    def alpha(self) -> float:
        """The vertical angle between the principal point and the pixel"""
        bp = self.height / 2.0
        ip = bp - self.y
        return math.atan(ip * math.tan(self.camera.alpha / 2.0) / bp)

    def beta(self) -> float:
        """The horizontal angle between the principal point and the pixel"""
        gp = self.width / 2.0
        ip = gp - self.x
        # flip sign (+ is right)
        return math.atan(ip * math.tan(self.camera.beta / 2.0) / gp) * -1

    def xDistance(self) -> float:
        """The side-to-side distance from the central meridian (a line extending from the point directly
        beneath the camera on the view plane through the principal point)"""
        oiDistance = self.camera.height / \
            math.sin(self.camera.theta - self.alpha())
        return oiDistance * math.tan(self.beta())

    def yDistance(self) -> float:
        """The forward distance from the the point on the view plane directly under the camera to the
        row in our image containing the pixel."""
        return self.camera.height / math.tan(self.camera.theta - self.alpha())

    @staticmethod
    def imageCorners(camera: Camera, width: int, height: int):
        return [Pixel(camera, width, height, 0, 0),
                Pixel(camera, width, height, width - 1, 0),
                Pixel(camera, width, height, width - 1, height - 1),
                Pixel(camera, width, height, 0, height - 1)]


def area(pixels: List[Pixel]) -> float:
    """Calculates the areas of a polygon, as defined by a list of pixels representing the
    vertices of the polygon"""
    n = len(pixels)
    p1 = 0
    p2 = 0
    for i1 in range(0, n):
        i2 = i1 + 1
        if i2 >= n:
            i2 = 0
        pixel1 = pixels[i1]
        pixel2 = pixels[i2]
        p1 = p1 + (pixel1.xDistance() * pixel2.yDistance())
        p2 = p2 + (pixel2.xDistance() * pixel1.yDistance())
    return abs((p1 - p2) / 2.0)


if __name__ == "__main__":
    width = 1782
    height = 1239
    c = Camera.fromDegrees(260, 35, 49, 31)
    pixels = [Pixel(c, width, height, 0, 38),
              Pixel(c, width, height, width - 1, 38),
              Pixel(c, width, height, width - 1, 93),
              Pixel(c, width, height, 0, 93),
              Pixel(c, width, height, 0, 38)]

    a = area(pixels)
    print(a)