adding more image processors for spot analysis

contrib/common/lib/cv/annotations/MomentsAnnotation.py

@@ -0,0 +1,222 @@
+import copy
+from functools import cache, cached_property
+
+import numpy as np
+import scipy.spatial.transform
+
+from opencsp.common.lib.cv.annotations.AbstractAnnotations import AbstractAnnotations
+import opencsp.common.lib.geometry.LineXY as l2
+import opencsp.common.lib.geometry.Pxy as p2
+import contrib.common.lib.geometry.RectXY as r2
+import opencsp.common.lib.geometry.RegionXY as reg
+import opencsp.common.lib.render.Color as color
+import opencsp.common.lib.render_control.RenderControlFigureRecord as rcfr
+import opencsp.common.lib.render_control.RenderControlPointSeq as rcps
+import opencsp.common.lib.tool.image_tools as it
+import opencsp.common.lib.tool.log_tools as lt
+
+
+class MomentsAnnotation(AbstractAnnotations):
+    def __init__(
+        self,
+        moments,
+        centroid_style: rcps.RenderControlPointSeq = None,
+        rotation_style: rcps.RenderControlPointSeq = None,
+    ):
+        """
+        centroid_style: RenderControlPointSeq, optional
+            Style used for render the centroid point. By default
+            RenderControlPointSeq.defualt(color=magenta).
+        rotation_style: RenderControlPointSeq, optional
+            Style used for render the rotation line. By default centroid_style.
+        """
+        if centroid_style is None:
+            centroid_style = rcps.default(color=color.magenta())
+        if rotation_style is None:
+            rotation_style = copy.deepcopy(centroid_style)
+        super().__init__(centroid_style)
+
+        self.moments = moments
+        self.rotation_style = rotation_style
+
+    def get_bounding_box(self, index=0) -> reg.RegionXY:
+        """
+        Returns a rectangle with 0 area and the same point for the
+        upper-left/lower-right corners. Moments don't have a bounding box.
+        """
+        return r2.RectXY(self.centroid, self.centroid)
+
+    @property
+    def origin(self) -> p2.Pxy:
+        return self.centroid
+
+    @cached_property
+    def rotation_angle_2d(self) -> float:
+        # from https://en.wikipedia.org/wiki/Image_moment
+        u00, u20, u02, u11 = (
+            self.central_moment(0, 0),
+            self.central_moment(2, 0),
+            self.central_moment(0, 2),
+            self.central_moment(1, 1),
+        )
+        up20, up02, up11 = u20 / u00, u02 / u00, u11 / u00
+
+        # phi = 0.5 * np.arctan((2*up11) / (up20 - up02))
+        phi = 0.5 * np.arctan2((2 * up11), (up20 - up02))
+
+        # convert to our standard angle coordinates, with 0 being on the
+        # positive x-axis and the angle increasing counter-clockwise
+        if phi < 0 and phi >= -np.pi / 2:
+            phi = -phi
+        elif phi <= np.pi / 2:
+            phi = (np.pi / 2 - phi) + np.pi / 2
+        else:
+            lt.error_and_raise(
+                RuntimeError,
+                "Error in MomentsAnnotation.rotation_angle_2d: "
+                + f"expected moments angle from atan2'/:: to be between -π/2 and π/2, but got {phi}.",
+            )
+
+        return phi
+
+    @property
+    def rotation(self) -> scipy.spatial.transform.Rotation:
+        """
+        The orientation of the primary axis of a spot in the 2d image, in
+        radians, about the z-axis.
+
+        This is a decent approximation of actual orientation for images with a
+        single spot that is shaped like an enlongated circle.
+        """
+        phi = self.rotation_angle_2d
+        return scipy.spatial.transform.Rotation.from_euler('z', [phi])
+
+    @property
+    def size(self) -> list[float]:
+        """Returns 0, always."""
+        return [0]
+
+    @property
+    def cX(self) -> float:
+        """Centroid X value"""
+        return self.moments["m10"] / self.moments["m00"]
+
+    @property
+    def cY(self) -> float:
+        """Centroid Y value"""
+        return self.moments["m01"] / self.moments["m00"]
+
+    @cached_property
+    def centroid(self) -> p2.Pxy:
+        """
+        The centroid of the image. The centroid is similar in concept to the
+        center of mass.
+        """
+        return p2.Pxy([[self.cX], [self.cY]])
+
+    @cached_property
+    def eccentricity_untested(self) -> float:
+        """
+        How elongated the image is.
+
+        TODO The results from this method are untested. Someone who knows what
+        this value means should contribute unit tests.
+        """
+        # from https://en.wikipedia.org/wiki/Image_moment
+        u00, u20, u02, u11 = (
+            self.central_moment(0, 0),
+            self.central_moment(2, 0),
+            self.central_moment(0, 2),
+            self.central_moment(1, 1),
+        )
+        up20, up02, up11 = u20 / u00, u02 / u00, u11 / u00
+        eigenval_1 = ((up20 + up02) / 2) + (np.sqrt(4 * up11 * up11 + (up20 - up02) ** 2) / 2)
+        eigenval_2 = ((up20 + up02) / 2) - (np.sqrt(4 * up11 * up11 + (up20 - up02) ** 2) / 2)
+        eccentricity = np.sqrt(1 - (eigenval_2 / eigenval_1))
+        return eccentricity
+
+    @cache
+    def central_moment(self, p: int, q: int) -> float:
+        """
+        Returns the central moment for the given order p and order q.
+        """
+        # from https://en.wikipedia.org/wiki/Image_moment
+        if p == 0 and q == 0:
+            return self.moments["m00"]
+        elif p == 0 and q == 1:
+            return 0
+        elif p == 1 and q == 0:
+            return 0
+        elif p == 1 and q == 1:
+            return self.moments["m11"] - self.cY * self.moments["m10"]
+        elif p == 2 and q == 0:
+            return self.moments["m20"] - self.cX * self.moments["m10"]
+        elif p == 0 and q == 2:
+            return self.moments["m02"] - self.cY * self.moments["m01"]
+        elif p == 2 and q == 1:
+            m21, m11, m20, m01 = self.moments["m21"], self.moments["m11"], self.moments["m20"], self.moments["m01"]
+            return m21 - 2 * self.cX * m11 - self.cY * m20 + 2 * self.cX * self.cX * m01
+        elif p == 1 and q == 2:
+            m12, m11, m02, m10 = self.moments["m12"], self.moments["m11"], self.moments["m02"], self.moments["m10"]
+            return m12 - 2 * self.cY * m11 - self.cX * m02 + 2 * self.cY * self.cY * m10
+        elif p == 3 and q == 0:
+            m30, m20, m10 = self.moments["m30"], self.moments["m20"], self.moments["m10"]
+            return m30 - 3 * self.cX * m20 + 2 * self.cX * self.cX * m10
+        elif p == 0 and q == 3:
+            m03, m02, m01 = self.moments["m03"], self.moments["m02"], self.moments["m01"]
+            return m03 - 3 * self.cY * m02 + 2 * self.cY * self.cY * m01
+        else:
+            lt.error_and_raise(
+                ValueError,
+                "Error in MomentsAnnotation.central_moment(): "
+                + f"formula for central moment with order (p={p}, q={q}) has not been implemented.",
+            )
+
+    # other uses of moments...
+
+    def render_to_figure(
+        self, fig_record: rcfr.RenderControlFigureRecord, image: np.ndarray = None, include_label: bool = False
+    ):
+        # draw the centroid marker
+        label = None if not include_label else "centroid"
+        fig_record.view.draw_pq(([self.cX], [self.cY]), self.style, label=label)
+
+        # start by assuming that the plot is >= 30 pixels
+        height, width, rotation_arrow_dist = None, None, 30
+        if image is not None:
+            height, width = it.dims_and_nchannels(image)[0]
+            rotation_arrow_dist = min(width, height) * 0.2
+        rot_mat_2d = self.rotation.as_matrix()[:, :2, :2].squeeze()
+        rot_mat_2d_rev = rot_mat_2d.transpose()
+        rotation_endpoints_1 = p2.Pxy([rotation_arrow_dist, 0]).rotate(rot_mat_2d) + self.centroid
+        rotation_endpoints_2 = p2.Pxy([rotation_arrow_dist, 0]).rotate(rot_mat_2d_rev) + self.centroid
+        rotation_endpoints_list = [(pnt.x[0], pnt.y[0]) for pnt in [rotation_endpoints_1, rotation_endpoints_2]]
+
+        # Use the bounds of the image in determining where to put the end point
+        # of the rotation arrow.
+        if image is not None:
+            image_tl = p2.Pxy([0, 0])
+            image_br = p2.Pxy([width, height])
+            if width > 80 and height > 80:
+                if self.cX > 20 and self.cX < width - 20 and self.cY > 20 and self.cY < height - 20:
+                    # Find where the rotation line intersects close to the edge of the
+                    # image. We use a small buffer from the image edge in order to
+                    # prevent Matplotlib from drawing outside the bounds of the image.
+                    image_tl = p2.Pxy([20, 20])
+                    image_br = p2.Pxy([width - 20, height - 20])
+            image_bounds = r2.RectXY(image_tl, image_br)
+            angle = self.rotation_angle_2d
+            rotation_line = l2.LineXY.from_location_angle(self.centroid, -angle)  # images use an inverted y-axis
+            intersections = p2.Pxy(image_bounds.loops[0].intersect_line(rotation_line))
+
+            if len(intersections) > 0:
+                rotation_endpoints_list = [(intersections.x[i], intersections.y[i]) for i in range(len(intersections))]
+
+        # draw the rotation as an arrow
+        style = copy.deepcopy(self.rotation_style)
+        style.marker = "arrow"
+        style.linestyle = "-"
+        fig_record.view.draw_pq_list(rotation_endpoints_list, style=style)
+
+    def __str__(self):
+        return f"MomentsAnnotation"

contrib/common/lib/cv/annotations/RectangleAnnotations.py

@@ -0,0 +1,140 @@
+import matplotlib.axes
+import matplotlib.collections
+import matplotlib.patches
+import numpy as np
+import scipy.spatial.transform
+
+from opencsp.common.lib.cv.annotations.AbstractAnnotations import AbstractAnnotations
+import opencsp.common.lib.geometry.LoopXY as l2
+import opencsp.common.lib.geometry.Pxy as p2
+import opencsp.common.lib.geometry.RegionXY as reg
+import opencsp.common.lib.render.Color as color
+import opencsp.common.lib.render.figure_management as fm
+import opencsp.common.lib.render.view_spec as vs
+import opencsp.common.lib.render.View3d as v3d
+import opencsp.common.lib.render_control.RenderControlAxis as rca
+import opencsp.common.lib.render_control.RenderControlPointSeq as rcps
+import opencsp.common.lib.tool.log_tools as lt
+
+
+class RectangleAnnotations(AbstractAnnotations):
+    """
+    A collection of pixel bounding boxes where points of interest are located in an image.
+    """
+
+    def __init__(
+        self,
+        style: rcps.RenderControlPointSeq = None,
+        upperleft_lowerright_corners: tuple[p2.Pxy, p2.Pxy] = None,
+        pixels_to_meters: float = None,
+    ):
+        """
+        Parameters
+        ----------
+        style : RenderControlBcs, optional
+            The rendering style, by default {magenta, no corner markers}
+        upperleft_lowerright_corners : Pxy
+            The upper-left and lower-right corners of the bounding box for this rectangle, in pixels
+        pixels_to_meters : float, optional
+            A simple conversion method for how many meters a pixel represents,
+            for use in scale(). By default None.
+        """
+        if style is None:
+            style = rcps.default(marker=None, color=color.magenta())
+        super().__init__(style)
+
+        self.points = upperleft_lowerright_corners
+        self.pixels_to_meters = pixels_to_meters
+
+    def get_bounding_box(self, index=0) -> reg.RegionXY:
+        x1 = self.points[0].x[index]
+        y1 = self.points[0].y[index]
+        x2 = self.points[1].x[index]
+        y2 = self.points[1].y[index]
+        w = x2 - x1
+        h = y2 - y1
+
+        return reg.RegionXY(l2.LoopXY.from_rectangle(x1, y1, w, h))
+
+    @property
+    def origin(self) -> p2.Pxy:
+        return self.points[0]
+
+    @property
+    def rotation(self) -> scipy.spatial.transform.Rotation:
+        raise NotImplementedError("Orientation is not yet implemented for RectangleAnnotations")
+
+    @property
+    def size(self) -> list[float]:
+        width_height = self.points[1] - self.points[0]
+        max_size: np.ndarray = np.max(width_height.data, axis=0)
+        return max_size.tolist()
+
+    @property
+    def scale(self) -> list[float]:
+        if self.pixels_to_meters is None:
+            lt.error_and_raise(
+                RuntimeError,
+                "Error in RectangeAnnotations.scale(): "
+                + "no pixels_to_meters conversion ratio is set, so scale can't be estimated",
+            )
+        return [self.size * self.pixels_to_meters]
+
+    def render_to_figure(self, fig: rcfr.RenderControlFigureRecord, image: np.ndarray = None, include_label=False):
+        label = self.get_label(include_label)
+
+        # get the corner vertices for each bounding box
+        draw_loops: list[list[tuple[int, int]]] = []
+        for index in range(len(self.points[0])):
+            bbox = self.get_bounding_box(index)
+            for loop in bbox.loops:
+                loop_verts = list(zip(loop.vertices.x, loop.vertices.y))
+                loop_verts = [(int(x), int(y)) for x, y in loop_verts]
+                draw_loops.append(loop_verts)
+
+        # draw the bounding boxes
+        for i, draw_loop in enumerate(draw_loops):
+            fig.view.draw_pq_list(draw_loop, close=True, style=self.style, label=label)
+            label = None
+
+
+if __name__ == "__main__":
+    from PIL import Image
+
+    import opencsp.common.lib.opencsp_path.opencsp_root_path as orp
+    from opencsp.common.lib.opencsp_path import opencsp_settings
+    import opencsp.common.lib.tool.file_tools as ft
+    import opencsp.common.lib.render.Color as color
+    import opencsp.common.lib.render.figure_management as fm
+    import opencsp.common.lib.render.view_spec as vs
+    import opencsp.common.lib.render_control.RenderControlAxis as rca
+    import opencsp.common.lib.render_control.RenderControlFigure as rcfg
+    import opencsp.common.lib.render_control.RenderControlFigureRecord as rcfr
+    import opencsp.common.lib.render_control.RenderControlPointSeq as rcps
+
+    src_dir = ft.join(
+        opencsp_settings["opencsp_root_path"]["collaborative_dir"],
+        "NSTTF_Optics/Experiments/2024-05-22_SolarNoonTowerTest_4/2_Data/BCS_data/20240522_115941 TowerTest4_NoSun/Raw Images",
+    )
+    src_file = "20240522_115941.65 TowerTest4_NoSun Raw.JPG"
+
+    # Prepare the rectangles
+    ul1, lr1 = p2.Pxy((447, 898)), p2.Pxy((519, 953))
+    ul2, lr2 = p2.Pxy((1158, 877)), p2.Pxy((1241, 935))
+    ul_lr_corners = (ul1.concatenate(ul2), lr1.concatenate(lr2))
+    rectangles = RectangleAnnotations(upperleft_lowerright_corners=ul_lr_corners)
+
+    # Load and update the image
+    image = np.array(Image.open(ft.join(src_dir, src_file)))
+    image2 = rectangles.render_to_image(image)
+
+    # Draw the image using figure_management
+    axis_control = rca.image(grid=False)
+    figure_control = rcfg.RenderControlFigure()
+    view_spec_2d = vs.view_spec_im()
+    fig_record = fm.setup_figure(
+        figure_control, axis_control, view_spec_2d, title="Flashlight ROIs", code_tag=f"{__file__}", equal=False
+    )
+    fig_record.view.imshow(image2)
+    fig_record.view.show(block=True)
+    fig_record.close()