Merge pull request #118 from braden6521/visualize_slope_map

Visualize slope map

opencsp/common/lib/csp/VisualizeOrthorectifiedSlopeAbstract.py

@@ -8,6 +8,8 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
+from opencsp.common.lib.csp.visualize_orthorectified_image import add_quivers, plot_orthorectified_image
+
 
 class VisualizeOrthorectifiedSlopeAbstract:
     """Abstract class inherited by all objects which can have orthorectified slope
@@ -23,6 +25,31 @@ def orthorectified_slope_array(self, x_vec: np.ndarray, y_vec: np.ndarray) -> np
     def axis_aligned_bounding_box(self) -> tuple[float, float, float, float]:
         pass
 
+    def get_orthorectified_slope_array(self, res) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+        """Returns slope array given resolution
+
+        Parameters
+        ----------
+        res : float, optional
+            The xy resolution of the plot, meters, by default 0.1
+
+        Returns
+        -------
+        slope_array : ndarray
+            Shape (2, n, m) array of x/y slopes
+        x_vec : ndarray
+            X interpolation vector
+        y_vec : ndarray
+            Y interpolation vector
+        """
+        # Create interpolation axes
+        left, right, bottom, top = self.axis_aligned_bounding_box
+        x_vec = np.arange(left, right, res)  # meters
+        y_vec = np.arange(bottom, top, res)  # meters
+
+        # Calculate slope image
+        return self.orthorectified_slope_array(x_vec, y_vec), x_vec, y_vec
+
     def plot_orthorectified_slope_error(
         self,
         reference: 'VisualizeOrthorectifiedSlopeAbstract',
@@ -115,11 +142,11 @@ def plot_orthorectified_slope_error(
 
         # Plot image on axis
         extent = (left - res / 2, right + res / 2, bottom - res / 2, top + res / 2)
-        self._plot_orthorectified_image(image, axis, cmap, extent, clims, 'mrad')
+        plot_orthorectified_image(image, axis, cmap, extent, clims, 'mrad')
 
         # Add quiver arrows
         if quiver_density is not None:
-            self._add_quivers(x_image, y_image, x_vec, y_vec, quiver_density, axis, quiver_scale, quiver_color)
+            add_quivers(x_image, y_image, x_vec, y_vec, quiver_density, axis, quiver_scale, quiver_color)
 
         # Label axes
         axis.set_title(title)
@@ -164,13 +191,8 @@ def plot_orthorectified_slope(
         if axis is None:
             axis = plt.gca()
 
-        # Create interpolation axes
-        left, right, bottom, top = self.axis_aligned_bounding_box
-        x_vec = np.arange(left, right, res)  # meters
-        y_vec = np.arange(bottom, top, res)  # meters
-
         # Calculate slope image
-        slopes = self.orthorectified_slope_array(x_vec, y_vec)
+        slopes, x_vec, y_vec = self.get_orthorectified_slope_array(res)
 
         # Calculate slope image
         if type_ == 'x':
@@ -202,12 +224,13 @@ def plot_orthorectified_slope(
             title = 'Slope Magnitude'
 
         # Plot image on axes
+        left, right, bottom, top = self.axis_aligned_bounding_box
         extent = (left - res / 2, right + res / 2, bottom - res / 2, top + res / 2)
-        self._plot_orthorectified_image(image, axis, 'jet', extent, clims, 'mrad')
+        plot_orthorectified_image(image, axis, 'jet', extent, clims, 'mrad')
 
         # Add quiver arrows
         if quiver_density is not None:
-            self._add_quivers(x_image, y_image, x_vec, y_vec, quiver_density, axis, quiver_scale, quiver_color)
+            add_quivers(x_image, y_image, x_vec, y_vec, quiver_density, axis, quiver_scale, quiver_color)
 
         # Label axes
         axis.set_title(title)
@@ -281,70 +304,7 @@ def plot_orthorectified_curvature(
             extent = (left, right, bottom, top)
 
         # Plot image on axes
-        self._plot_orthorectified_image(image, axis, 'seismic', extent, clims, 'mrad/meter')
+        plot_orthorectified_image(image, axis, 'seismic', extent, clims, 'mrad/meter')
 
         # Label axes
         axis.set_title(title)
-
-    def _add_quivers(
-        self,
-        im_x: np.ndarray,
-        im_y: np.ndarray,
-        x_vec: np.ndarray,
-        y_vec: np.ndarray,
-        quiver_density: float,
-        axis: plt.Axes | None = None,
-        scale: float | None = None,
-        color: str = 'white',
-    ) -> None:
-        """
-        Adds quiver arrows to data plot.
-
-        Parameters
-        ----------
-        im_x/im_y : ndarray
-            Images to sample x/y quiver directions from.
-        x_vec/y_vec : ndarray
-            X and Y data grid axes, meters.
-        quiver_density : float
-            Spacing of quiver arrows in meters.
-        axis : [plt.Axes | None], optional
-            Axes to plot on. The default is None. If None, uses plt.gca().
-        scale : [float | None], optional
-            Matplotlib "scale" for adding quiver arrows. The default is None.
-            If None, uses the default scale.
-        color : str
-            Color of the quiver arrows.
-        """
-        if axis is None:
-            axis = plt.gca()
-
-        # Calculate quiver points
-        res_x = np.mean(np.abs(np.diff(x_vec)))
-        res_y = np.mean(np.abs(np.diff(y_vec)))
-        Nx = int(quiver_density / res_x)
-        Ny = int(quiver_density / res_y)
-        x1 = int(Nx / 2)
-        y1 = int(Ny / 2)
-
-        x_locs, y_locs = np.meshgrid(x_vec[x1::Nx], y_vec[y1::Ny])
-        u_dirs = -im_x[y1::Ny, x1::Nx]
-        v_dirs = -im_y[y1::Ny, x1::Nx]
-
-        # Add quiver arrows to axes
-        axis.quiver(x_locs, y_locs, u_dirs, v_dirs, color=color, scale=scale, scale_units='x')
-
-    def _plot_orthorectified_image(
-        self,
-        image: np.ndarray,
-        axis: plt.Axes,
-        cmap: str,
-        extent: tuple[float, float, float, float],
-        clims: tuple[float, float],
-        cmap_title: str,
-    ):
-        """Plots orthorectified image on axes"""
-        plt_im = axis.imshow(image, cmap, origin='lower', extent=extent)
-        plt_im.set_clim(clims)
-        plt_cmap = plt.colorbar(plt_im, ax=axis)
-        plt_cmap.ax.set_ylabel(cmap_title, rotation=270, labelpad=15)

opencsp/common/lib/csp/test/test_VisualizeOrthorectifiedImage.py

@@ -0,0 +1,145 @@
+from os.path import join, dirname
+import unittest
+
+import matplotlib.pyplot as plt
+import matplotlib.testing.compare as mplt
+import numpy as np
+
+from opencsp.common.lib.csp.MirrorPoint import MirrorPoint
+from opencsp.common.lib.geometry.RegionXY import RegionXY
+from opencsp.common.lib.geometry.Pxyz import Pxyz
+from opencsp.common.lib.geometry.Vxy import Vxy
+from opencsp.common.lib.geometry.Uxyz import Uxyz
+import opencsp.common.lib.tool.file_tools as ft
+
+
+class TestVisualizeOrthorectifiedSlopeAbstract(unittest.TestCase):
+    """Tests orthorectified plots of
+    - Slope
+    - Slope error
+    - Curvature
+
+    NOTE: To update the unit test data, run the test and copy the .PNG files from
+    the data/output folder to the data/input folder. Run the test again to confirm passing.
+    """
+
+    @classmethod
+    def setUpClass(cls) -> None:
+        # Set up directories
+        cls.dir_output = join(dirname(__file__), 'data/output/VisualizeOrthorectifiedSlopeAbstract')
+        ft.create_directories_if_necessary(cls.dir_output)
+        cls.dir_input = join(dirname(__file__), 'data/input/VisualizeOrthorectifiedSlopeAbstract')
+        ft.create_directories_if_necessary(cls.dir_input)
+
+        # Define optic shape
+        shape = RegionXY.from_vertices(Vxy(([0.6, -0.6, -0.6, 0.6], [-0.6, -0.6, 0.6, 0.6])))
+
+        # Calculate surface xyz points
+        xv = yv = np.arange(-0.6, 0.7, 0.1)
+        X, Y = np.meshgrid(xv, yv)
+        Z = np.zeros(X.shape)
+        surface_points = Pxyz((X, Y, Z))
+
+        # Calculate normal vectors
+        nvecs = np.ones((3, len(surface_points)))
+        nvecs[0] = np.sin(2 * np.pi * X).flatten() * 0.05
+        nvecs[1] = np.sin(2 * np.pi * Y).flatten() * 0.05
+        normal_vectors = Uxyz(nvecs)
+
+        # Create mirror object
+        cls.test_mirror_bilinear = MirrorPoint(surface_points, normal_vectors, shape, 'bilinear')
+        cls.test_mirror_nearest = MirrorPoint(surface_points, normal_vectors, shape, 'nearest')
+
+        # Create reference optic
+        nvecs_flat = np.ones((3, len(surface_points)))
+        nvecs_flat[0:2] = 0
+        normal_vectors_flat = Uxyz(nvecs_flat)
+        cls.reference = MirrorPoint(surface_points, normal_vectors_flat, shape, 'nearest')
+
+    def _get_axes(self) -> tuple[plt.Figure, plt.Axes]:
+        fig = plt.figure(figsize=(6, 6))
+        return fig, fig.gca()
+
+    def test_plot_slope_error_magnitude_linear(self):
+        # Setup files
+        file_out = join(self.dir_output, 'slope_error_linear.png')
+        file_in = join(self.dir_input, 'slope_error_linear.png')
+        # Create image
+        fig, ax = self._get_axes()
+        self.test_mirror_bilinear.plot_orthorectified_slope_error(
+            self.reference, 0.005, 'magnitude', 100, ax, 0.1, 1000, 'black'
+        )
+        fig.savefig(file_out, dpi=300)
+        plt.close(fig)
+        # Test
+        self._compare_actual_expected_images(file_out, file_in)
+
+    def test_plot_slope_error_magnitude_nearest(self):
+        # Setup files
+        file_out = join(self.dir_output, 'slope_error_nearest.png')
+        file_in = join(self.dir_input, 'slope_error_nearest.png')
+        # Create image
+        fig, ax = self._get_axes()
+        self.test_mirror_nearest.plot_orthorectified_slope_error(
+            self.reference, 0.005, 'magnitude', 100, ax, 0.1, 1000, 'black'
+        )
+        fig.savefig(file_out, dpi=300)
+        plt.close(fig)
+        # Test
+        self._compare_actual_expected_images(file_out, file_in)
+
+    def test_plot_slope_magnitude_linear(self):
+        # Setup files
+        file_out = join(self.dir_output, 'slope_linear.png')
+        file_in = join(self.dir_input, 'slope_linear.png')
+        # Create image
+        fig, ax = self._get_axes()
+        self.test_mirror_bilinear.plot_orthorectified_slope(0.005, 'magnitude', 100, ax, 0.1, 1000, 'black')
+        fig.savefig(file_out, dpi=300)
+        plt.close(fig)
+        # Test
+        self._compare_actual_expected_images(file_out, file_in)
+
+    def test_plot_slope_magnitude_nearest(self):
+        # Setup files
+        file_out = join(self.dir_output, 'slope_nearest.png')
+        file_in = join(self.dir_input, 'slope_nearest.png')
+        # Create image
+        fig, ax = self._get_axes()
+        self.test_mirror_nearest.plot_orthorectified_slope(0.005, 'magnitude', 100, ax, 0.1, 1000, 'black')
+        fig.savefig(file_out, dpi=300)
+        plt.close(fig)
+        # Test
+        self._compare_actual_expected_images(file_out, file_in)
+
+    def test_plot_curvature_magnitude_linear(self):
+        # Setup files
+        file_out = join(self.dir_output, 'curvature_linear.png')
+        file_in = join(self.dir_input, 'curvature_linear.png')
+        # Create image
+        fig, ax = self._get_axes()
+        self.test_mirror_bilinear.plot_orthorectified_curvature(0.005, 'combined', 100, ax)
+        fig.savefig(file_out, dpi=300)
+        plt.close(fig)
+        # Test
+        self._compare_actual_expected_images(file_out, file_in)
+
+    def test_plot_curvature_magnitude_nearest(self):
+        # Setup files
+        file_out = join(self.dir_output, 'curvature_nearest.png')
+        file_in = join(self.dir_input, 'curvature_nearest.png')
+        # Create image
+        fig, ax = self._get_axes()
+        self.test_mirror_nearest.plot_orthorectified_curvature(0.005, 'combined', 100, ax)
+        fig.savefig(file_out, dpi=300)
+        plt.close(fig)
+        # Test
+        self._compare_actual_expected_images(file_out, file_in)
+
+    def _compare_actual_expected_images(self, actual_location: str, expected_location: str, tolerance=0.2) -> bool:
+        """Tests if image files match."""
+        self.assertIsNone(mplt.compare_images(expected_location, actual_location, tolerance))
+
+
+if __name__ == '__main__':
+    unittest.main()

opencsp/common/lib/csp/visualize_orthorectified_image.py

@@ -0,0 +1,83 @@
+"""Library that handles plotting orthorectified images (slope images, etc.)"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+
+def add_quivers(
+    im_x: np.ndarray,
+    im_y: np.ndarray,
+    x_vec: np.ndarray,
+    y_vec: np.ndarray,
+    quiver_density: float,
+    axis: plt.Axes | None = None,
+    scale: float | None = None,
+    color: str = 'white',
+) -> None:
+    """
+    Adds quiver arrows to data plot.
+
+    Parameters
+    ----------
+    im_x/im_y : ndarray
+        Images to sample x/y quiver directions from.
+    x_vec/y_vec : ndarray
+        X and Y data grid axes, meters.
+    quiver_density : float
+        Spacing of quiver arrows in meters.
+    axis : [plt.Axes | None], optional
+        Axes to plot on. The default is None. If None, uses plt.gca().
+    scale : [float | None], optional
+        Matplotlib "scale" for adding quiver arrows. The default is None.
+        If None, uses the default scale.
+    color : str
+        Color of the quiver arrows.
+    """
+    if axis is None:
+        axis = plt.gca()
+
+    # Calculate quiver points
+    res_x = np.mean(np.abs(np.diff(x_vec)))
+    res_y = np.mean(np.abs(np.diff(y_vec)))
+    Nx = int(quiver_density / res_x)
+    Ny = int(quiver_density / res_y)
+    x1 = int(Nx / 2)
+    y1 = int(Ny / 2)
+
+    x_locs, y_locs = np.meshgrid(x_vec[x1::Nx], y_vec[y1::Ny])
+    u_dirs = -im_x[y1::Ny, x1::Nx]
+    v_dirs = -im_y[y1::Ny, x1::Nx]
+
+    # Add quiver arrows to axes
+    axis.quiver(x_locs, y_locs, u_dirs, v_dirs, color=color, scale=scale, scale_units='x')
+
+
+def plot_orthorectified_image(
+    image: np.ndarray,
+    axis: plt.Axes,
+    cmap: str,
+    extent: tuple[float, float, float, float],
+    clims: tuple[float, float],
+    cmap_title: str,
+):
+    """Plots orthorectified image on axes
+
+    Parameters
+    ----------
+    image : np.ndarray
+        2d image to plot
+    axis : plt.Axes
+        Matplotlib axis to plot on
+    cmap : str
+        Color map
+    extent : tuple[float, float, float, float]
+        Left, right, bottom, top
+    clims : tuple[float, float]
+        Color bar limits [low, high]
+    cmap_title : str
+        Title of colorbar
+    """
+    plt_im = axis.imshow(image, cmap, origin='lower', extent=extent)
+    plt_im.set_clim(clims)
+    plt_cmap = plt.colorbar(plt_im, ax=axis)
+    plt_cmap.ax.set_ylabel(cmap_title, rotation=270, labelpad=15)