diff --git a/contrib/test_data_generation/sofast_fringe/generate_test_data_multi_facet.py b/contrib/test_data_generation/sofast_fringe/generate_test_data_multi_facet.py index 51e7a3e88..d2acc6346 100644 --- a/contrib/test_data_generation/sofast_fringe/generate_test_data_multi_facet.py +++ b/contrib/test_data_generation/sofast_fringe/generate_test_data_multi_facet.py @@ -73,7 +73,7 @@ def generate_dataset( image = np.zeros(mask.shape) * np.nan for idx in range(sofast.num_facets): mask = sofast.data_image_processing_facet[idx].mask_processed - slopes_xy = sofast.data_characterization_facet[idx].slopes_facet_xy + slopes_xy = sofast.data_calculation_facet[idx].slopes_facet_xy slopes = np.sqrt(np.sum(slopes_xy**2, 0)) image[mask] = slopes plt.figure() diff --git a/contrib/test_data_generation/sofast_fringe/generate_test_data_single_facet.py b/contrib/test_data_generation/sofast_fringe/generate_test_data_single_facet.py index f92c41a90..f3a550df8 100644 --- a/contrib/test_data_generation/sofast_fringe/generate_test_data_single_facet.py +++ b/contrib/test_data_generation/sofast_fringe/generate_test_data_single_facet.py @@ -3,7 +3,7 @@ """ from os.path import join, dirname, exists -from typing import Literal +from typing import Literal, Final import matplotlib.pyplot as plt import numpy as np @@ -73,8 +73,9 @@ def generate_dataset( print(f'All data saved to: {file_dataset_out:s}') # Show slope map - mask = sofast.data_image_processing_facet[0].mask_processed - slopes_xy = sofast.data_characterization_facet[0].slopes_facet_xy + idx_facet: Final = 0 + mask = sofast.data_image_processing_facet[idx_facet].mask_processed + slopes_xy = sofast.data_calculation_facet[idx_facet].slopes_facet_xy slopes = np.sqrt(np.sum(slopes_xy**2, 0)) image = np.zeros(mask.shape) * np.nan image[mask] = slopes diff --git a/contrib/test_data_generation/sofast_fringe/generate_test_data_undefined.py b/contrib/test_data_generation/sofast_fringe/generate_test_data_undefined.py index 1282f6f05..b6ba77c26 100644 --- a/contrib/test_data_generation/sofast_fringe/generate_test_data_undefined.py +++ b/contrib/test_data_generation/sofast_fringe/generate_test_data_undefined.py @@ -54,7 +54,7 @@ def generate_dataset( # Show slope map mask = sofast.data_image_processing_facet[0].mask_processed - slopes_xy = sofast.data_characterization_facet[0].slopes_facet_xy + slopes_xy = sofast.data_calculation_facet[0].slopes_facet_xy slopes = np.sqrt(np.sum(slopes_xy**2, 0)) image = np.zeros(mask.shape) * np.nan image[mask] = slopes diff --git a/example/sofast_fringe/example_process_facet_ensemble.py b/example/sofast_fringe/example_process_facet_ensemble.py index 1e1ce4c3c..4cf83c4f3 100644 --- a/example/sofast_fringe/example_process_facet_ensemble.py +++ b/example/sofast_fringe/example_process_facet_ensemble.py @@ -8,6 +8,7 @@ from opencsp.app.sofast.lib.MeasurementSofastFringe import MeasurementSofastFringe from opencsp.app.sofast.lib.ProcessSofastFringe import ProcessSofastFringe from opencsp.app.sofast.lib.SpatialOrientation import SpatialOrientation +from opencsp.app.sofast.lib.SofastConfiguration import SofastConfiguration from opencsp.common.lib.camera.Camera import Camera from opencsp.common.lib.csp.FacetEnsemble import FacetEnsemble from opencsp.common.lib.deflectometry.Surface2DParabolic import Surface2DParabolic @@ -132,21 +133,16 @@ def example_process_facet_ensemble(): data['robust_least_squares'] = bool(data['robust_least_squares']) surfaces.append(Surface2DParabolic(**data)) - # Update search parameters - # sofast.params.mask_hist_thresh = 0.83 - # sofast.params.geometry.perimeter_refine_perpendicular_search_dist = 15.0 - # sofast.params.geometry.facet_corns_refine_frac_keep = 1.0 - # sofast.params.geometry.facet_corns_refine_perpendicular_search_dist = 3.0 - # sofast.params.geometry.facet_corns_refine_step_length = 5.0 - # Process sofast.process_optic_multifacet(facet_data, ensemble_data, surfaces) # 3. Log best-fit parabolic focal lengths # ======================================= - for idx in range(sofast.num_facets): - surf_coefs = sofast.data_characterization_facet[idx].surf_coefs_facet - focal_lengths_xy = [1 / 4 / surf_coefs[2], 1 / 4 / surf_coefs[5]] + sofast_config = SofastConfiguration() + sofast_config.load_sofast_object(sofast) + sofast_stats = sofast_config.get_measurement_stats() + for stat in sofast_stats: + focal_lengths_xy = stat['focal_lengths_parabolic_xy'] lt.info(f'Facet {idx:d} xy focal lengths (meters): {focal_lengths_xy[0]:.3f}, {focal_lengths_xy[1]:.3f}') # 4. Plot slope magnitude @@ -165,8 +161,10 @@ def example_process_facet_ensemble(): axis_control_m = rca.meters() # Plot slope map + res = 0.002 # meter, make the plot with 2mm spatial resolution + clim = 7 # mrad, draw the plot with +/-7mrad scale bars, this mirror has erorrs that extend to about +/-7mrad fig_record = fm.setup_figure(figure_control, axis_control_m, title='') - ensemble.plot_orthorectified_slope(res=0.002, clim=7, axis=fig_record.axis) + ensemble.plot_orthorectified_slope(res=res, clim=clim, axis=fig_record.axis) fig_record.save(dir_save, 'slope_magnitude', 'png') # 5. Plot 3d representation of facet ensemble diff --git a/example/sofast_fringe/example_process_single_facet.py b/example/sofast_fringe/example_process_single_facet.py index 3974ea6d6..a2772487a 100644 --- a/example/sofast_fringe/example_process_single_facet.py +++ b/example/sofast_fringe/example_process_single_facet.py @@ -74,7 +74,7 @@ def example_process_single_facet(): # 3. Log best-fit parabolic focal lengths # ======================================= - surf_coefs = sofast.data_characterization_facet[0].surf_coefs_facet + surf_coefs = sofast.data_calculation_facet[0].surf_coefs_facet focal_lengths_xy = [1 / 4 / surf_coefs[2], 1 / 4 / surf_coefs[5]] lt.info(f'Facet xy focal lengths (meters): ' f'{focal_lengths_xy[0]:.3f}, {focal_lengths_xy[1]:.3f}') diff --git a/example/sofast_fringe/example_process_undefined_shape.py b/example/sofast_fringe/example_process_undefined_shape.py index cae41551e..45efa07f0 100644 --- a/example/sofast_fringe/example_process_undefined_shape.py +++ b/example/sofast_fringe/example_process_undefined_shape.py @@ -72,7 +72,7 @@ def example_process_undefined_shape_facet(): # 3. Log best-fit parabolic focal lengths # ======================================= - surf_coefs = sofast.data_characterization_facet[0].surf_coefs_facet + surf_coefs = sofast.data_calculation_facet[0].surf_coefs_facet focal_lengths_xy = [1 / 4 / surf_coefs[2], 1 / 4 / surf_coefs[5]] lt.info(f'Facet xy focal lengths (meters): ' f'{focal_lengths_xy[0]:.3f}, {focal_lengths_xy[1]:.3f}') diff --git a/opencsp/app/sofast/lib/BlobIndex.py b/opencsp/app/sofast/lib/BlobIndex.py index 0425bed9f..226522c16 100644 --- a/opencsp/app/sofast/lib/BlobIndex.py +++ b/opencsp/app/sofast/lib/BlobIndex.py @@ -1,12 +1,24 @@ +from enum import Enum from typing import Literal import matplotlib.pyplot as plt import numpy as np from opencsp.common.lib.geometry.Vxy import Vxy +from opencsp.common.lib.geometry.LoopXY import LoopXY from opencsp.common.lib.tool import log_tools as lt +class Step(Enum): + """Gives step direction + - Left/right for horizontal searches + - Up/down for vertical searches + """ + + RIGHT_OR_DOWN = 1 + LEFT_OR_UP = -1 + + class BlobIndex: """Class containing blob indexing algorithms to assign indices to blobs in a rough grid pattern. X/Y axes correspond to image axes; +x is to right, +y is down. Class takes in points (in units @@ -16,12 +28,15 @@ class BlobIndex: Attributes ---------- search_thresh : float - + Threshold, in pixels. As algorithm calculates the expected location of the next blob, if a blob + is within "search_thresh" of the expected location, that blob is considered found. + max_num_iters : int + Maximum number of iterations to use when fitting found blobs to grid pattern. search_perp_axis_ratio : float Ratio of point distances: (perpendicular to axis) / (along axis) used to search for points. apply_filter : bool - To filter bad points (experimental) + To filter bad points (experimental, not implemented yet) """ def __init__(self, points: Vxy, x_min: int, x_max: int, y_min: int, y_max: int) -> 'BlobIndex': @@ -46,9 +61,15 @@ def __init__(self, points: Vxy, x_min: int, x_max: int, y_min: int, y_max: int) self._is_assigned = np.zeros(self._num_pts, dtype=bool) self._neighbor_dists = np.zeros((self._num_pts, 4)) * np.nan # left, right, up, down - self.search_thresh = 5.0 # pixels - self.search_perp_axis_ratio = 3.0 - self.apply_filter = False + self.search_thresh: float = 5.0 # pixels + """Threshold, in pixels. As algorithm calculates the expected location of the next blob, if a blob + is within 'search_thresh' of the expected location, that blob is considered found. Default 5.0""" + self.max_num_iters: int = 100 + """Maximum number of iterations to use when fitting found blobs to grid pattern. Default 100""" + self.search_perp_axis_ratio: float = 3.0 + """Ratio of point distances: (perpendicular to axis) / (along axis) used to search for points. Default 3.0""" + self.apply_filter: bool = False + """To filter bad points (experimental, not implemented yet). Default False""" self._offset_x = -x_min # index self._offset_y = -y_min # index @@ -131,17 +152,14 @@ def _point_index_from_xy_index(self, idx_x: int, idx_y: int) -> tuple[bool, int] return (not np.isnan(idx)), int(idx) def _assign(self, idx_pt: int, idx_x: int, idx_y: int) -> None: - """Assigns given blob index an xy index + # Assigns given blob index an xy index + # idx_pt : int + # Index of point (indexing self._points) + # idx_x : int + # X blob index + # idx_y : int + # Y blob index - Parameters - ---------- - idx_pt : int - Index of point (indexing self._points) - idx_x : int - X blob index - idx_y : int - Y blob index - """ # Assign vectors self._idx_x[idx_pt] = idx_x self._idx_y[idx_pt] = idx_y @@ -153,13 +171,8 @@ def _assign(self, idx_pt: int, idx_x: int, idx_y: int) -> None: lt.debug(f'Blob number {idx_pt:d} was assigned ({idx_x:d}, {idx_y:d})') def _unassign(self, idx_pt: int) -> None: - """Unassigns a point index + # Unassigns a point index - Parameters - ---------- - idx_pt : int - Index of point (indexing self._points) - """ # Unassign matrices idx_mat_x = self._idx_x[idx_pt] + self._offset_x idx_mat_y = self._idx_y[idx_pt] + self._offset_y @@ -177,18 +190,12 @@ def _unassign(self, idx_pt: int) -> None: lt.debug(f'Blob number {idx_pt:d} was unassigned') - def _assign_center(self, pt_origin: Vxy) -> None: - """Assigns the center point to (0, 0) - - Parameters - ---------- - pt_origin : Vxy - Location of origin, pixels - """ + def _assign_center(self, pt_origin: Vxy, idx_x: int, idx_y: int) -> None: + # Assigns given origin point to given xy index idx, dist = self._nearest_unassigned_idx_from_xy_point(pt_origin) if dist > self.search_thresh: - lt.warn(f'Assigning point {idx:d} to index (0, 0) resulted in {dist:.2f} pixels error.') - self._assign(idx, 0, 0) + lt.warn(f'Assigning point {idx:d} to index ({idx_x:d}, {idx_y:d}) resulted in {dist:.2f} pixels error.') + self._assign(idx, idx_x, idx_y) def _find_nearest_in_direction( self, idx_pt: int, direction: Literal['right', 'left', 'up', 'down'] @@ -248,15 +255,14 @@ def _num_unassigned(self) -> int: """Returns number of unassigned points (referencing self._points)""" return np.logical_not(self._is_assigned).sum() - def _find_3x3_center_block(self) -> None: - """Finds the center 3x3 block around center point using nearest in direction method - a b c - d e f - g h i - """ - ret, idx_e = self._point_index_from_xy_index(0, 0) + def _find_3x3_center_block(self, idx_x: int, idx_y: int) -> None: + # Finds the center 3x3 block around center point using nearest in direction method + # a b c + # d e f + # g h i + ret, idx_e = self._point_index_from_xy_index(idx_x, idx_y) if not ret: - lt.warn('Could not find 3x3 center block. Could not find point index (0, 0).') + lt.warn(f'Could not find 3x3 center block. Could not find point index xy=({idx_x:d}, {idx_y:d}).') # Right idx_f, x, y = self._find_nearest_in_direction(idx_e, 'right') self._assign(idx_f, x, y) @@ -282,22 +288,18 @@ def _find_3x3_center_block(self) -> None: idx_g, x, y = self._find_nearest_in_direction(idx_d, 'down') self._assign(idx_g, x, y) - def _extend_data(self, direction: Literal['x', 'y'], step: Literal[1, -1]) -> None: - """Extends found blob rows/collumns in given direction - - Steps in the given axis, a or b - - Parameters - ---------- - direction : Literal['x', 'y'] - Axis to search - step : Literal[1, -1] - Direction to search - - 1 = right/down - - -1 = left/up - """ - if step not in [-1, 1]: - raise ValueError(f'Step must be -1 or 1, not {step}') + def _extend_data(self, direction: Literal['x', 'y'], step: Step) -> None: + # Extends found blob rows/collumns in given direction + # Steps in the given axis, a or b + # + # Parameters + # ---------- + # direction : Literal['x', 'y'] + # Axis to search + # step : Step + # Direction to search + if not isinstance(step, Step): + raise ValueError(f'Step must be a Step class, not type {type(step)}') if direction == 'x': idxs_a = self._idx_y @@ -317,8 +319,8 @@ def _extend_data(self, direction: Literal['x', 'y'], step: Literal[1, -1]) -> No is_b = idxs_b[mask] # indices of points on axis # Step through all points on axis for i_b in is_b: - if not i_b + step in is_b: # If adjacent point is not assigned, find it - idx_b_prev = i_b - step # Index used for slope calc + if not i_b + step.value in is_b: # If adjacent point is not assigned, find it + idx_b_prev = i_b - step.value # Index used for slope calc if idx_b_prev in is_b: # If history exists, find points for idx_b_next in range(500): # First iteration, use previously assigned points @@ -344,36 +346,33 @@ def _extend_data(self, direction: Literal['x', 'y'], step: Literal[1, -1]) -> No # Assign point if dist < self.search_thresh: if direction == 'x': - idx_x = int(i_b + step + (idx_b_next * step)) + idx_x = int(i_b + step.value + (idx_b_next * step.value)) idx_y = int(idx_a) else: idx_x = int(idx_a) - idx_y = int(i_b + step + (idx_b_next * step)) + idx_y = int(i_b + step.value + (idx_b_next * step.value)) self._assign(idx_new, idx_x, idx_y) else: break def _exp_pt_from_pt_pair(self, pt_cur: Vxy, pt_prev: Vxy) -> Vxy: - """Calculates the expected point from a given current and previous point pair - - Parameters - ---------- - pt_cur : Vxy - Current point, pixels - pt_prev : Vxy - Previous point, pixels - - Returns - ------- - Vxy - Refined expected point location - """ + # Calculates the expected point from a given current and previous point pair + # Parameters + # ---------- + # pt_cur : Vxy + # Current point, pixels + # pt_prev : Vxy + # Previous point, pixels + # Returns + # ------- + # Vxy + # Refined expected point location del_y = (pt_cur.y - pt_prev.y)[0] # pixels del_x = (pt_cur.x - pt_prev.x)[0] # pixels return pt_cur + Vxy((del_x, del_y)) # pixels def _filter_bad_points(self) -> None: - """Filters erroneous assigned points""" + # Filters erroneous assigned points del_1_x = np.diff(self._points_mat, axis=1) del_1_y = np.diff(self._points_mat, axis=0) del_2_x = np.diff(del_1_x, axis=1) @@ -397,27 +396,29 @@ def _filter_bad_points(self) -> None: for idx in idxs_bad_pts: self._unassign(idx) - def run(self, pt_origin: Vxy) -> None: + def run(self, pt_known: Vxy, x_known: int, y_known: int) -> None: """Runs blob indexing sequence Parameters ---------- - pt_origin : Vxy - Location of origin point with blob index of (0, 0), pixels + pt_known : Vxy + Location of known point with known blob index, pixels + x/y_known : int + XY indies of known points """ # Assign center point - self._assign_center(pt_origin) + self._assign_center(pt_known, x_known, y_known) # Find 3x3 core point block - self._find_3x3_center_block() + self._find_3x3_center_block(x_known, y_known) # Extend rows prev_num_unassigned = self._num_unassigned() - for idx in range(100): - self._extend_data('x', -1) - self._extend_data('x', 1) + for idx in range(self.max_num_iters): + self._extend_data('x', Step(-1)) + self._extend_data('x', Step(1)) if self.apply_filter: self._filter_bad_points() - self._extend_data('y', -1) - self._extend_data('y', 1) + self._extend_data('y', Step(-1)) + self._extend_data('y', Step(1)) if self.apply_filter: self._filter_bad_points() @@ -426,8 +427,7 @@ def run(self, pt_origin: Vxy) -> None: if prev_num_unassigned == cur_num_unassigned: lt.debug(f'All possible points found in {idx + 1:d} iterations.') break - else: - prev_num_unassigned = cur_num_unassigned + prev_num_unassigned = cur_num_unassigned def plot_points_labels(self, labels: bool = False) -> None: """Plots points and labels @@ -497,6 +497,25 @@ def get_data(self) -> tuple[Vxy, Vxy]: points = Vxy((x_pts[mask_assigned], y_pts[mask_assigned])) return points, indices + def get_data_in_region(self, loop: LoopXY) -> tuple[Vxy, Vxy]: + """Returns found points and indices within given region + + Parameters + ---------- + loop : LoopXY + The loop to search within + + Returns + ------- + points : Vxy + Lenght N vector, located points xy locations, pixels + indices_xy : Vxy + Length N vector, located points xy blob indices, int + """ + points, indices = self.get_data() + mask = loop.is_inside(points) + return points[mask], indices[mask] + def get_data_mat(self) -> tuple[np.ndarray, np.ndarray]: """Returns found points and indices in matrix form diff --git a/opencsp/app/sofast/lib/CalibrateSofastFixedDots.py b/opencsp/app/sofast/lib/CalibrateSofastFixedDots.py index c5629808c..395a1beb6 100644 --- a/opencsp/app/sofast/lib/CalibrateSofastFixedDots.py +++ b/opencsp/app/sofast/lib/CalibrateSofastFixedDots.py @@ -61,6 +61,7 @@ def __init__( x_max: int, y_min: int, y_max: int, + xy_pts_known: tuple[tuple[int, int]] = None, ) -> 'CalibrateSofastFixedDots': """Instantiates the calibration class @@ -69,7 +70,7 @@ def __init__( files_images : list[str] File paths to images of Aruco markers and dots origin_pts : Vxy - Points on the images corresponding to index (0, 0) + Points on the images corresponding to index `xy_pts_known` camera : Camera Camera calibration parameters of camera used to take images pts_xyz_corners : Vxyz @@ -80,6 +81,10 @@ def __init__( Expected min/max x index values (follows screen x axis) y_min/y_max : int Expected min/max y index values (follows screen y axis) + xy_pts_known : tuple[tuple[int, int]] + The origin point xy indices. One xy point for each input image. + For example, for two images using (0, 0) - xy_pts_known = ((0, 0), (0, 0)). + If None, defaults to all zeros. """ # Load images self._images: list[ImageMarker] = [] @@ -92,8 +97,12 @@ def __init__( for pt_xyz, pt_id in zip(pts_xyz_corners, pts_ids_corners): im.set_point_id_located(pt_id, pt_xyz.data.squeeze()) + if xy_pts_known is None: + xy_pts_known = ((0, 0),) * len(origin_pts) + # Save data self._origin_pts = origin_pts + self._xy_pts_known = xy_pts_known self._camera = camera self._pts_xyz_corners = pts_xyz_corners self._pts_ids_corners = pts_ids_corners @@ -150,7 +159,9 @@ def _find_dots_in_images(self) -> None: # Index all found points blob_index = BlobIndex(pts, -self._x_max, -self._x_min, self._y_min, self._y_max) blob_index.search_thresh = self.blob_search_threshold - blob_index.run(origin_pt) + blob_index.run( + origin_pt, x_known=self._xy_pts_known[idx_image][0], y_known=self._xy_pts_known[idx_image][1] + ) points, indices = blob_index.get_data_mat() # Save points and indices diff --git a/opencsp/app/sofast/lib/ProcessSofastAbstract.py b/opencsp/app/sofast/lib/ProcessSofastAbstract.py new file mode 100644 index 000000000..a7badf09e --- /dev/null +++ b/opencsp/app/sofast/lib/ProcessSofastAbstract.py @@ -0,0 +1,298 @@ +""" +Abstract module for processing SOFAST data. + +This module provides the `ProcessSofastAbstract` class, which defines common attributes and methods +used for processing SOFAST data. It includes functionality for saving data to HDF5 files, calculating +facet pointing, and generating OpenCSP representations of the optic under test. + +Notes +----- +- The `ProcessSofastAbstract` class is designed to be subclassed and extended with specific + processing logic. + +ChatGPT 4o assisted in generating some docstrings in this module. +""" + +from typing import Literal + +import numpy as np + + +import opencsp.app.sofast.lib.calculation_data_classes as cdc +from opencsp.app.sofast.lib.DefinitionEnsemble import DefinitionEnsemble +from opencsp.app.sofast.lib.DefinitionFacet import DefinitionFacet +from opencsp.common.lib.csp.Facet import Facet +from opencsp.common.lib.csp.FacetEnsemble import FacetEnsemble +from opencsp.common.lib.csp.MirrorPoint import MirrorPoint +from opencsp.common.lib.deflectometry.SlopeSolverData import SlopeSolverData +from opencsp.common.lib.deflectometry.Surface2DAbstract import Surface2DAbstract +from opencsp.common.lib.geometry.RegionXY import RegionXY +from opencsp.common.lib.geometry.TransformXYZ import TransformXYZ +from opencsp.common.lib.geometry.Uxyz import Uxyz +from opencsp.common.lib.geometry.Vxy import Vxy +from opencsp.common.lib.geometry.Vxyz import Vxyz +from opencsp.common.lib.tool.hdf5_tools import HDF5_SaveAbstract +import opencsp.common.lib.tool.log_tools as lt + + +class ProcessSofastAbstract(HDF5_SaveAbstract): + """ + Abstract class for ProcessSofast classes defining common attributes and methods. + + Attributes + ---------- + num_facets : int + Number of facets in the optic. + optic_type : Literal['undefined', 'single', 'multi'] + Type of optic being processed. + data_facet_def : list[DefinitionFacet] + List of facet definitions. + data_ensemble_def : DefinitionEnsemble + Ensemble definition of the optic. + data_surfaces : list[Surface2DAbstract] + List of surface definitions. + data_geometry_general : cdc.CalculationDataGeometryGeneral + General geometry data for the calculation. + data_image_processing_general : cdc.CalculationImageProcessingGeneral + General image processing data for the calculation. + data_geometry_facet : list[cdc.CalculationDataGeometryFacet] + List of geometry data for each facet. + data_image_processing_facet : list[cdc.CalculationImageProcessingFacet] + List of image processing data for each facet. + data_error : cdc.CalculationError + Error data for the calculation. + data_calculation_facet : list[SlopeSolverData] + List of slope solver data for each facet. + data_calculation_ensemble : list[cdc.CalculationFacetEnsemble] + List of calculation data for the ensemble. + params : HDF5_SaveAbstract + Parameters for the calculation. + + Methods + ------- + save_to_hdf(file: str, prefix: str = '') + Saves data to the given HDF file. + get_optic(interp_type: Literal['bilinear', 'clough_tocher', 'nearest'] = 'nearest') -> FacetEnsemble | Facet + Returns the OpenCSP representation of the optic under test. + """ + + def __init__(self): + self.num_facets: int = None + self.optic_type: Literal['undefined', 'single', 'multi'] = None + self.data_facet_def: list[DefinitionFacet] = None + self.data_ensemble_def: DefinitionEnsemble = None + + self.data_surfaces: list[Surface2DAbstract] = None + + self.data_geometry_general: cdc.CalculationDataGeometryGeneral = None + self.data_image_processing_general: cdc.CalculationImageProcessingGeneral = None + self.data_geometry_facet: list[cdc.CalculationDataGeometryFacet] = None + self.data_image_processing_facet: list[cdc.CalculationImageProcessingFacet] = None + self.data_error: cdc.CalculationError = None + + self.data_calculation_facet: list[SlopeSolverData] = None + self.data_calculation_ensemble: list[cdc.CalculationFacetEnsemble] = None + + self.params: HDF5_SaveAbstract = None + + def save_to_hdf(self, file: str, prefix: str = ''): + """Saves data to given file. Data is stored as: PREFIX + Folder/Field_1 + + Parameters + ---------- + file : str + HDF file to save to + prefix : str + Prefix to append to folder path within HDF file (folders must be separated by "/") + """ + # Log + lt.info(f'Saving Sofast data to: {file:s}, in HDF5 folder: "{prefix:s}"') + + # One per measurement + if self.data_error is not None: + self.data_error.save_to_hdf(file, f'{prefix:s}DataSofastCalculation/general/') + self.data_geometry_general.save_to_hdf(file, f'{prefix:s}DataSofastCalculation/general/') + self.data_image_processing_general.save_to_hdf(file, f'{prefix:s}DataSofastCalculation/general/') + + # Sofast parameters + self.params.save_to_hdf(file, f'{prefix:s}DataSofastInput/') + + # Facet definition + if self.data_facet_def is not None: + for idx_facet, facet_data in enumerate(self.data_facet_def): + facet_data.save_to_hdf(file, f'{prefix:s}DataSofastInput/optic_definition/facet_{idx_facet:03d}/') + + # Ensemble definition + if self.data_ensemble_def is not None: + self.data_ensemble_def.save_to_hdf(file, f'{prefix:s}DataSofastInput/optic_definition/') + + # Surface definition + for idx_facet, surface in enumerate(self.data_surfaces): + surface.save_to_hdf(file, f'{prefix:s}DataSofastInput/optic_definition/facet_{idx_facet:03d}/') + + # Calculations, one per facet + for idx_facet in range(self.num_facets): + # Save facet slope data + if self.data_calculation_facet is not None: + self.data_calculation_facet[idx_facet].save_to_hdf( + file, f'{prefix:s}DataSofastCalculation/facet/facet_{idx_facet:03d}/' + ) + # Save facet geometry data + if self.data_geometry_facet is not None: + self.data_geometry_facet[idx_facet].save_to_hdf( + file, f'{prefix:s}DataSofastCalculation/facet/facet_{idx_facet:03d}/' + ) + # Save facet image processing data + if self.data_image_processing_facet is not None: + self.data_image_processing_facet[idx_facet].save_to_hdf( + file, f'{prefix:s}DataSofastCalculation/facet/facet_{idx_facet:03d}/' + ) + # Save ensemle data + if self.data_calculation_ensemble is not None: + self.data_calculation_ensemble[idx_facet].save_to_hdf( + file, f'{prefix:s}DataSofastCalculation/facet/facet_{idx_facet:03d}/' + ) + + def _calculate_facet_pointing(self, reference: Literal['average'] | int = 'average') -> None: + """ + Calculates facet pointing relative to the given reference. + + Parameters + ---------- + reference : 'average' | int + If 'average', the pointing reference is the average of all + facet pointing directions. If, int, that facet index is assumed + to have perfect pointing. + """ + if self.data_calculation_facet is None: + lt.error_and_raise(ValueError, 'Slopes must be solved first by running "solve_slopes".') + if reference != 'average' and not isinstance(reference, int): + lt.error_and_raise(ValueError, 'Given reference must be int or "average".') + if isinstance(reference, int) and reference >= self.num_facets: + lt.error_and_raise( + ValueError, f'Given facet index, {reference:d}, is out of range of 0-{self.num_facets - 1:d}.' + ) + + # Instantiate data list + self.data_calculation_ensemble = [] + + trans_facet_ensemble_list = [] + v_pointing_matrix = np.zeros((3, self.num_facets)) + for idx in range(self.num_facets): + # Get transformation from user-input and slope solving + trans_1 = TransformXYZ.from_R_V( + self.data_ensemble_def.r_facet_ensemble[idx], self.data_ensemble_def.v_facet_locations[idx] + ) + trans_2 = self.data_calculation_facet[idx].trans_alignment + # Calculate inverse of slope solving transform + trans_2 = TransformXYZ.from_V(-trans_2.V) * TransformXYZ.from_R(trans_2.R.inv()) + # Create local to global transformation + trans_facet_ensemble_list.append(trans_2 * trans_1) + + # Calculate pointing vector in ensemble coordinates + v_pointing = Vxyz((0, 0, 1)).rotate(trans_facet_ensemble_list[idx].R) + v_pointing_matrix[:, idx] = v_pointing.data.squeeze() + + # Calculate reference pointing direction + if isinstance(reference, int): + v_pointing_ref = Vxyz(v_pointing_matrix[:, reference]) + elif reference == 'average': + v_pointing_ref = Vxyz(v_pointing_matrix.mean(1)) + # Calculate rotation to align pointing vectors + r_align_pointing = v_pointing_ref.align_to(Vxyz((0, 0, 1))) + trans_align_pointing = TransformXYZ.from_R(r_align_pointing) + + # Apply alignment rotation to total transformation + trans_facet_ensemble_list = [trans_align_pointing * t for t in trans_facet_ensemble_list] + + # Calculate global slope and surface points + for idx in range(self.num_facets): + # Get slope data + slopes = self.data_calculation_facet[idx].slopes_facet_xy # facet coordinats + + # Calculate surface normals in local (facet) coordinates + number_data_points = slopes.shape[1] + u_surf_norms = np.ones((3, number_data_points)) + u_surf_norms[:2] = -slopes + u_surf_norms = Uxyz(u_surf_norms).as_Vxyz() + + # Apply rotation to normal vectors + u_surf_norms_global = u_surf_norms.rotate(trans_facet_ensemble_list[idx].R) + # Convert normal vectors to global (ensemble) slopes + slopes_ensemble_xy = -u_surf_norms_global.projXY().data / u_surf_norms_global.z + + # Convert surface points to global (ensemble) coordinates + v_surf_points_ensemble = trans_facet_ensemble_list[idx].apply( + self.data_calculation_facet[idx].v_surf_points_facet + ) + + # Calculate pointing vectors in ensemble coordinates + v_facet_pointing_ensemble = Vxyz((0, 0, 1)).rotate(trans_facet_ensemble_list[idx].R) + + data = cdc.CalculationFacetEnsemble( + trans_facet_ensemble_list[idx], slopes_ensemble_xy, v_surf_points_ensemble, v_facet_pointing_ensemble + ) + self.data_calculation_ensemble.append(data) + + def get_optic( + self, interp_type: Literal['bilinear', 'clough_tocher', 'nearest'] = 'nearest' + ) -> FacetEnsemble | Facet: + """Returns the OpenCSP representation of the optic under test. Returns either + a Facet or FacetEnsemble object depending on the optic type. Each mirror is + represented by a MirrorPoint object. Each mirror origin is co-located with its + parent facet origin. + + Parameters + ---------- + interp_type : {'bilinear', 'clough_tocher', 'nearest'}, optional + Mirror interpolation type, by default 'nearest' + + Returns + ------- + FacetEnsemble | Facet + Optic object + """ + if self.data_calculation_facet is None: + lt.error_and_raise(ValueError, 'Sofast data must be processed before optic is available.') + + facets = [] + trans_list = [] + for idx_mirror in range(self.num_facets): + # Get surface points + pts: Vxyz = self.data_calculation_facet[idx_mirror].v_surf_points_facet + # Get normals from slopes + slopes: np.ndarray = self.data_calculation_facet[idx_mirror].slopes_facet_xy + number_data_points = slopes.shape[1] + norm_data = np.ones((3, number_data_points)) + norm_data[:2] = -slopes + norm_vecs = Uxyz(norm_data) + # Get mirror shape + if self.optic_type == 'undefined': + # Find bounding box + x1 = pts.x.min() + x2 = pts.x.max() + y1 = pts.y.min() + y2 = pts.y.max() + vertices = Vxy(([x1, x1, x2, x2], [y1, y2, y2, y1])) + else: + # Get optic region from optic definition + vertices = self.data_facet_def[idx_mirror].v_facet_corners.projXY() + shape = RegionXY.from_vertices(vertices) + # Create mirror + mirror = MirrorPoint(pts, norm_vecs, shape, interp_type) + # Create facet + facet = Facet(mirror) + # Locate facet within ensemble + if self.optic_type == 'multi': + trans: TransformXYZ = self.data_calculation_ensemble[idx_mirror].trans_facet_ensemble + trans_list.append(trans) + # Save facets + facets.append(facet) + + # Return optics + if self.optic_type == 'multi': + ensemble = FacetEnsemble(facets) + ensemble.set_facet_transform_list(trans_list) + return ensemble + else: + return facets[0] diff --git a/opencsp/app/sofast/lib/ProcessSofastFixed.py b/opencsp/app/sofast/lib/ProcessSofastFixed.py index 1a8986e03..95ae9d58a 100644 --- a/opencsp/app/sofast/lib/ProcessSofastFixed.py +++ b/opencsp/app/sofast/lib/ProcessSofastFixed.py @@ -1,41 +1,32 @@ """Class that handles the processing of fixed pattern deflectometry data. """ -from typing import Literal - import cv2 as cv import numpy as np from numpy import ndarray from opencsp.app.sofast.lib.BlobIndex import BlobIndex -import opencsp.app.sofast.lib.calculation_data_classes as cdc +from opencsp.app.sofast.lib.DefinitionEnsemble import DefinitionEnsemble from opencsp.app.sofast.lib.DefinitionFacet import DefinitionFacet from opencsp.app.sofast.lib.DotLocationsFixedPattern import DotLocationsFixedPattern import opencsp.app.sofast.lib.image_processing as ip from opencsp.app.sofast.lib.MeasurementSofastFixed import MeasurementSofastFixed from opencsp.app.sofast.lib.ParamsSofastFixed import ParamsSofastFixed +from opencsp.app.sofast.lib.ProcessSofastAbstract import ProcessSofastAbstract import opencsp.app.sofast.lib.process_optics_geometry as pr from opencsp.app.sofast.lib.SpatialOrientation import SpatialOrientation from opencsp.common.lib.camera.Camera import Camera -from opencsp.common.lib.csp.MirrorPoint import MirrorPoint from opencsp.common.lib.deflectometry.SlopeSolver import SlopeSolver -from opencsp.common.lib.deflectometry.SlopeSolverData import SlopeSolverData from opencsp.common.lib.deflectometry.Surface2DAbstract import Surface2DAbstract -from opencsp.common.lib.geometry.RegionXY import RegionXY -from opencsp.common.lib.geometry.Uxyz import Uxyz -from opencsp.common.lib.tool.hdf5_tools import HDF5_SaveAbstract +from opencsp.common.lib.geometry.Vxy import Vxy import opencsp.common.lib.tool.log_tools as lt -class ProcessSofastFixed(HDF5_SaveAbstract): +class ProcessSofastFixed(ProcessSofastAbstract): """Fixed Pattern Deflectrometry data processing class""" def __init__( - self, - orientation: SpatialOrientation, - camera: Camera, - fixed_pattern_dot_locs: DotLocationsFixedPattern, - data_facet: DefinitionFacet, + self, orientation: SpatialOrientation, camera: Camera, fixed_pattern_dot_locs: DotLocationsFixedPattern ) -> 'ProcessSofastFixed': """Instantiates class @@ -47,17 +38,14 @@ def __init__( Camera object fixed_pattern_dot_locs : DotLocationsFixedPattern Image projection dictionary - data_facet : DefinitionFacet - DefinitionFacet object """ + super().__init__() + self.orientation = orientation self.camera = camera - self.fixed_pattern_dot_locs = fixed_pattern_dot_locs - self.data_facet = data_facet - self.params = ParamsSofastFixed() - - # Measurement data self.measurement: MeasurementSofastFixed + self.fixed_pattern_dot_locs = fixed_pattern_dot_locs + self.params: ParamsSofastFixed = ParamsSofastFixed() # Define blob detector self.blob_detector: cv.SimpleBlobDetector_Params = cv.SimpleBlobDetector_Params() @@ -69,39 +57,34 @@ def __init__( self.blob_detector.filterByConvexity = False self.blob_detector.filterByInertia = False - # Calculations - self.data_slope_solver: SlopeSolverData - self.data_geometry_general: cdc.CalculationDataGeometryGeneral - self.data_image_proccessing_general: cdc.CalculationImageProcessingGeneral - self.data_geometry_facet: list[cdc.CalculationDataGeometryFacet] - self.data_image_processing_facet: list[cdc.CalculationImageProcessingFacet] - self.data_error: cdc.CalculationError - self.blob_index: BlobIndex - self.slope_solver: SlopeSolver - - # Input parameters - self.data_surface: Surface2DAbstract = None - - def find_blobs(self) -> BlobIndex: - """Finds blobs in image""" + # Instantiate other data containers + self.slope_solvers: list[SlopeSolver] = None + self.blob_index: BlobIndex = None + + def find_blobs(self, pts_known: Vxy, xys_known: tuple[tuple[int, int]]) -> BlobIndex: + """Finds blobs in image + + Parameters + ---------- + pts_known : Vxy + Length N, xy pixel location of known point(s) with known xy dot index locations + xys_known : tuple[tuple[int]] + Length N integer xy dot indices + + NOTE: N=number of facets + """ pts_blob = ip.detect_blobs(self.measurement.image, self.blob_detector) # Index blobs blob_index = BlobIndex(pts_blob, *self.fixed_pattern_dot_locs.dot_extent) blob_index.search_thresh = self.params.blob_search_thresh blob_index.search_perp_axis_ratio = self.params.search_perp_axis_ratio - blob_index.run(self.measurement.origin) + for pt_known, xy_known in zip(pts_known, xys_known): + blob_index.run(pt_known, xy_known[0], xy_known[1]) return blob_index - def calculate_mask(self) -> ndarray: - """Calculate mask image - - Parameters - ---------- - image : ndarray - Measurement image - """ + def _calculate_mask(self) -> ndarray: # Calculate mask im_dark = self.measurement.image * 0 images = np.concatenate((im_dark[..., None], self.measurement.image[..., None]), axis=2) @@ -113,45 +96,37 @@ def calculate_mask(self) -> ndarray: ] mask = ip.calc_mask_raw(images, *params) - if self.params.mask.keep_largest_area: + if (self.optic_type == 'multi') and self.params.mask.keep_largest_area: + lt.warn( + '"keep_largest_area" mask processing option cannot be used ' + 'for multifacet ensembles. This will be turned off.' + ) + self.params.mask.keep_largest_area = False + elif self.params.mask.keep_largest_area: mask = ip.keep_largest_mask_area(mask) return mask - def generate_geometry(self, blob_index: BlobIndex, mask_raw: np.ndarray) -> dict: - """Generates blob dataset from sofast dataset. + def load_measurement_data(self, measurement: MeasurementSofastFixed) -> None: + """Saves measurement data in class Parameters ---------- - blob_index : BlobIndex - BlobIndex object with all blobs assigned - mask_raw : ndarray - Mask array - - Returns - ------- - dict - Key word argument dictionary for SlopeSolver (does not include "surface_data") - - Sets Attributes - --------------- - self.data_geometry_general - self.data_image_proccessing_general - self.data_geometry_facet - self.data_image_processing_facet - self.data_error + measurement: MeasurementSofastFixed + Fixed pattern measurement object """ - pts_image, pts_index_xy = blob_index.get_data() + self.measurement = measurement - # Process optic geometry + def _process_optic_singlefacet_geometry(self, blob_index: BlobIndex, mask_raw: np.ndarray) -> dict: + # Process optic geometry (find mask corners, etc.) ( self.data_geometry_general, - self.data_image_proccessing_general, - self.data_geometry_facet, - self.data_image_processing_facet, + self.data_image_processing_general, + self.data_geometry_facet, # list + self.data_image_processing_facet, # list self.data_error, ) = pr.process_singlefacet_geometry( - self.data_facet, + self.data_facet_def[0], mask_raw, self.measurement.v_measure_point_facet, self.measurement.dist_optic_screen, @@ -161,12 +136,15 @@ def generate_geometry(self, blob_index: BlobIndex, mask_raw: np.ndarray) -> dict self.params.debug_geometry, ) + # Get image points and blob indices + pts_image, pts_index_xy = blob_index.get_data() + # Define optic orientation w.r.t. camera rot_optic_cam = self.data_geometry_general.r_optic_cam_refine_1 v_cam_optic_cam = self.data_geometry_general.v_cam_optic_cam_refine_2 u_cam_measure_point_facet = self.data_geometry_facet[0].u_cam_measure_point_facet - # Get screen-camera pose + # Get screen/camera poses rot_cam_optic = rot_optic_cam.inv() rot_optic_screen = self.orientation.r_cam_screen * rot_optic_cam rot_screen_optic = rot_optic_screen.inv() @@ -183,94 +161,201 @@ def generate_geometry(self, blob_index: BlobIndex, mask_raw: np.ndarray) -> dict u_pixel_pointing_cam = self.camera.vector_from_pixel(pts_image) u_pixel_pointing_facet = u_pixel_pointing_cam.rotate(rot_cam_optic).as_Vxyz() - self.params.debug_slope_solver.optic_data = self.data_facet + # Update debug data + self.params.debug_slope_solver.optic_data = self.data_facet_def[0] + # Construct surface kwargs return { 'v_optic_cam_optic': v_optic_cam_optic, 'u_active_pixel_pointing_optic': u_pixel_pointing_facet, 'u_measure_pixel_pointing_optic': u_cam_measure_point_facet, 'v_screen_points_facet': v_screen_points_facet, 'v_optic_screen_optic': v_optic_screen_optic, - 'v_align_point_optic': self.data_facet.v_facet_centroid, + 'v_align_point_optic': self.data_facet_def[0].v_facet_centroid, 'dist_optic_screen': self.measurement.dist_optic_screen, 'debug': self.params.debug_slope_solver, + 'surface': self.data_surfaces[0], } - def load_measurement_data(self, measurement: MeasurementSofastFixed) -> None: - """Saves measurement data in class - - Parameters - ---------- - measurement: MeasurementSofastFixed - Fixed pattern measurement object - """ - self.measurement = measurement + def _process_optic_multifacet_geometry(self, blob_index: BlobIndex, mask_raw: np.ndarray) -> list[dict]: + # Process optic geometry (find mask corners, etc.) + ( + self.data_geometry_general, + self.data_image_processing_general, + self.data_geometry_facet, # list + self.data_image_processing_facet, # list + self.data_error, + ) = pr.process_multifacet_geometry( + self.data_facet_def, + self.data_ensemble_def, + mask_raw, + self.measurement.v_measure_point_facet, + self.orientation, + self.camera, + self.measurement.dist_optic_screen, + self.params.geometry, + self.params.debug_geometry, + ) - def process_single_facet_optic(self, surface: Surface2DAbstract) -> None: - """Processes single facet optic. Saves data to self.data_slope_solver + kwargs_list = [] + for idx_facet in range(self.num_facets): + # Get pixel region of current facet + loop = self.data_image_processing_facet[idx_facet].loop_facet_image_refine + + # Get image points and blob indices + pts_image, pts_index_xy = blob_index.get_data_in_region(loop) + + # Define optic orientation w.r.t. camera + rot_facet_ensemble = self.data_ensemble_def.r_facet_ensemble[idx_facet] + rot_ensemble_cam = self.data_geometry_general.r_optic_cam_refine_2 + rot_facet_cam = rot_ensemble_cam * rot_facet_ensemble + + v_cam_ensemble_cam = self.data_geometry_general.v_cam_optic_cam_refine_3 + v_ensemble_facet_ensemble = self.data_ensemble_def.v_facet_locations[idx_facet] + v_ensemble_facet_cam = v_ensemble_facet_ensemble.rotate(rot_ensemble_cam) + v_cam_facet_cam = v_cam_ensemble_cam + v_ensemble_facet_cam + + u_cam_measure_point_facet = self.data_geometry_facet[idx_facet].u_cam_measure_point_facet + + # Get screen/camera poses + rot_cam_facet = rot_facet_cam.inv() + rot_facet_screen = self.orientation.r_cam_screen * rot_facet_cam + rot_screen_facet = rot_facet_screen.inv() + + v_facet_cam_facet = -v_cam_facet_cam.rotate(rot_cam_facet) + v_cam_screen_facet = self.orientation.v_cam_screen_cam.rotate(rot_cam_facet) + v_facet_screen_facet = v_facet_cam_facet + v_cam_screen_facet + + # Calculate xyz screen points + v_screen_points_screen = self.fixed_pattern_dot_locs.xy_indices_to_screen_coordinates(pts_index_xy) + v_screen_points_facet = v_facet_screen_facet + v_screen_points_screen.rotate(rot_screen_facet) + + # Calculate active pixel pointing + u_pixel_pointing_cam = self.camera.vector_from_pixel(pts_image) + u_pixel_pointing_facet = u_pixel_pointing_cam.rotate(rot_cam_facet).as_Vxyz() + + # Update debug data + self.params.debug_slope_solver.optic_data = self.data_facet_def[idx_facet] + + # Construct list of surface kwargs + kwargs_list.append( + { + 'v_optic_cam_optic': v_facet_cam_facet, + 'u_active_pixel_pointing_optic': u_pixel_pointing_facet, + 'u_measure_pixel_pointing_optic': u_cam_measure_point_facet, + 'v_screen_points_facet': v_screen_points_facet, + 'v_optic_screen_optic': v_facet_screen_facet, + 'v_align_point_optic': self.data_geometry_facet[idx_facet].v_align_point_facet, + 'dist_optic_screen': self.data_geometry_facet[idx_facet].measure_point_screen_distance, + 'debug': self.params.debug_slope_solver, + 'surface': self.data_surfaces[idx_facet], + } + ) + return kwargs_list + + def process_single_facet_optic( + self, data_facet_def: DefinitionFacet, surface: Surface2DAbstract, pt_known: Vxy, xy_known: tuple[int, int] + ) -> None: + """Processes single facet optic. Saves data to self.data_calculation_facet Parameters ---------- + data_facet_def : DefinitionFacet objec + Facet definition surface : Surface2DAbstract Surface 2d class + pt_known : Vxy + Length 1, xy pixel location of known point(s) with known xy dot index locations + xy_known : tuple[int, int] + Integer xy dot indices """ + + # Check inputs + if len(pt_known) != 1: + lt.error_and_raise( + ValueError, f'Only 1 pt_known can be given for single facet processing but {len(pt_known):d} were given' + ) + + self.optic_type = 'single' + self.num_facets = 1 + self.data_facet_def = [data_facet_def.copy()] + self.data_surfaces = [surface] + # Find blobs - self.blob_index = self.find_blobs() + self.blob_index = self.find_blobs(pt_known, (xy_known,)) # Calculate mask - mask_raw = self.calculate_mask() - mask_raw = ip.keep_largest_mask_area(mask_raw) + mask_raw = self._calculate_mask() # Generate geometry and slope solver inputs - kwargs = self.generate_geometry(self.blob_index, mask_raw) - - # Add surface fitting parameters - kwargs.update({'surface': surface}) + kwargs = self._process_optic_singlefacet_geometry(self.blob_index, mask_raw) # Calculate slope - self.slope_solver = SlopeSolver(**kwargs) - self.slope_solver.fit_surface() - self.slope_solver.solve_slopes() - self.data_slope_solver = self.slope_solver.get_data() - - self.data_surface = surface + slope_solver = SlopeSolver(**kwargs) + slope_solver.fit_surface() + slope_solver.solve_slopes() + self.slope_solvers = [slope_solver] + self.data_calculation_facet = [slope_solver.get_data()] - def save_to_hdf(self, file: str, prefix: str = ''): - """Saves data to given HDF5 file. Data is stored in CalculationsFixedPattern/... + def process_multi_facet_optic( + self, + data_facet_def: list[DefinitionFacet], + surfaces: list[Surface2DAbstract], + data_ensemble_def: DefinitionEnsemble, + pts_known: Vxy, + xys_known: tuple[tuple[int, int]], + ) -> None: + """Processes multi facet optic. Saves data to self.data_calculation_facet Parameters ---------- - file : str - HDF file to save to - prefix : str, optional - Prefix to append to folder path within HDF file (folders must be separated by "/"). - Default is empty string ''. + data_facet_def : list[DefinitionFacet] + List of facet data objects. + data_ensemble_def : DefinitionEnsemble + Ensemble data object. + surfaces : list[Surface2dAbstract] + List of surface type definitions + pts_known : Vxy + Length N, xy pixel location of known point(s) with known xy dot index locations + xys_known : tuple[tuple[int, int]] + List of N integer xy dot indices corresponding to pts_known + + NOTE: N=number of facets """ - # Sofast input - self.params.save_to_hdf(file, f'{prefix:s}DataSofastInput/') - self.data_surface.save_to_hdf(file, f'{prefix:s}DataSofastInput/optic_definition/facet_000/') - self.data_facet.save_to_hdf(file, f'{prefix:s}DataSofastInput/optic_definition/facet_000/') - - # General - self.data_error.save_to_hdf(file, f'{prefix:s}DataSofastCalculation/general/') - self.data_geometry_general.save_to_hdf(file, f'{prefix:s}DataSofastCalculation/general/') - self.data_image_proccessing_general.save_to_hdf(file, f'{prefix:s}DataSofastCalculation/general/') - - # Calculations - self.data_slope_solver.save_to_hdf(file, f'{prefix:s}DataSofastCalculation/facet/facet_000/') - self.data_geometry_facet[0].save_to_hdf(file, f'{prefix:s}DataSofastCalculation/facet/facet_000/') - self.data_image_processing_facet[0].save_to_hdf(file, f'{prefix:s}DataSofastCalculation/facet/facet_000/') - - lt.info(f'SofastFixed data saved to: {file:s} with prefix: {prefix:s}') - - def get_mirror( - self, interpolation_type: Literal['given', 'bilinear', 'clough_tocher', 'nearest'] = 'nearest' - ) -> MirrorPoint: - """Returns mirror object with slope data""" - v_surf_pts = self.data_slope_solver.v_surf_points_facet - v_normals_data = np.ones((3, len(v_surf_pts))) - v_normals_data[:2, :] = self.data_slope_solver.slopes_facet_xy - v_normals_data[:2, :] *= -1 - v_normals = Uxyz(v_normals_data) - shape = RegionXY.from_vertices(self.data_facet.v_facet_corners.projXY()) - return MirrorPoint(v_surf_pts, v_normals, shape, interpolation_type) + + # Check inputs + if len(data_facet_def) != len(surfaces) != len(pts_known) != len(xys_known): + lt.error_and_raise( + ValueError, + 'Length of data_facet_def does not equal length of data_surfaces' + + f'data_facet_def={len(data_facet_def)}, surface_data={len(surfaces)}, ' + + f'pts_known={len(pts_known)}, xys_known={len(xys_known)}', + ) + + self.optic_type = 'multi' + self.num_facets = len(data_facet_def) + self.data_facet_def = [d.copy() for d in data_facet_def] + self.data_ensemble_def = data_ensemble_def.copy() + self.data_surfaces = surfaces + + # Find blobs + self.blob_index = self.find_blobs(pts_known, xys_known) + + # Calculate mask + mask_raw = self._calculate_mask() + + # Generate geometry and slope solver inputs + kwargs_list = self._process_optic_multifacet_geometry(self.blob_index, mask_raw) + + # Calculate slope + self.slope_solvers = [] + self.data_calculation_facet = [] + for kwargs in kwargs_list: + slope_solver = SlopeSolver(**kwargs) + slope_solver.fit_surface() + slope_solver.solve_slopes() + self.slope_solvers.append(slope_solver) + self.data_calculation_facet.append(slope_solver.get_data()) + + # Calculate facet pointing + self._calculate_facet_pointing() diff --git a/opencsp/app/sofast/lib/ProcessSofastFringe.py b/opencsp/app/sofast/lib/ProcessSofastFringe.py index 6b9b6f34b..163d0e331 100644 --- a/opencsp/app/sofast/lib/ProcessSofastFringe.py +++ b/opencsp/app/sofast/lib/ProcessSofastFringe.py @@ -2,36 +2,25 @@ to calculate surface slopes. """ -from typing import Literal - import numpy as np -import opencsp.app.sofast.lib.calculation_data_classes as cdc from opencsp.app.sofast.lib.DefinitionEnsemble import DefinitionEnsemble from opencsp.app.sofast.lib.DefinitionFacet import DefinitionFacet from opencsp.app.sofast.lib.DisplayShape import DisplayShape as Display import opencsp.app.sofast.lib.image_processing as ip from opencsp.app.sofast.lib.MeasurementSofastFringe import MeasurementSofastFringe from opencsp.app.sofast.lib.ParamsSofastFringe import ParamsSofastFringe +from opencsp.app.sofast.lib.ProcessSofastAbstract import ProcessSofastAbstract import opencsp.app.sofast.lib.process_optics_geometry as po from opencsp.app.sofast.lib.SpatialOrientation import SpatialOrientation from opencsp.common.lib.camera.Camera import Camera -from opencsp.common.lib.csp.Facet import Facet -from opencsp.common.lib.csp.FacetEnsemble import FacetEnsemble -from opencsp.common.lib.csp.MirrorPoint import MirrorPoint from opencsp.common.lib.deflectometry.SlopeSolver import SlopeSolver -from opencsp.common.lib.deflectometry.SlopeSolverData import SlopeSolverData from opencsp.common.lib.deflectometry.Surface2DAbstract import Surface2DAbstract -from opencsp.common.lib.geometry.RegionXY import RegionXY -from opencsp.common.lib.geometry.TransformXYZ import TransformXYZ -from opencsp.common.lib.geometry.Uxyz import Uxyz from opencsp.common.lib.geometry.Vxy import Vxy -from opencsp.common.lib.geometry.Vxyz import Vxyz -from opencsp.common.lib.tool.hdf5_tools import HDF5_SaveAbstract import opencsp.common.lib.tool.log_tools as lt -class ProcessSofastFringe(HDF5_SaveAbstract): +class ProcessSofastFringe(ProcessSofastAbstract): """Class that processes measurement data captured by a SOFAST system. Computes optic surface slope and saves data to HDF5 format. @@ -63,8 +52,8 @@ class ProcessSofastFringe(HDF5_SaveAbstract): - data_error - errors between optic/sceario definitions and internal calculations - data_image_processing_general - general optic image processing calculations - data_image_processing_facet - facet specific image processing calculations - - data_characterization_facet - facet specific slope calculations in facet reference frame - - data_characterization_ensemble - facet specific slope/pointing calculations in ensemble reference frame + - data_calculation_facet - facet specific slope calculations in facet reference frame + - data_calculation_ensemble - facet specific slope/pointing calculations in ensemble reference frame External Data Storage --------------------- @@ -175,31 +164,13 @@ def __init__( display : Display Display object used to capture data. """ - # Store data + super().__init__() + + self.orientation = orientation + self.camera = camera self.measurement = measurement self.display = display - self.camera = camera - self.orientation = orientation - - # Define default calculation parameters - self.params = ParamsSofastFringe() - - # Instantiate data containers - self.num_facets: int = 0 - self.optic_type: Literal['undefined', 'single', 'multi'] = None - self.data_facet_def: list[DefinitionFacet] = None - self.data_ensemble_def: DefinitionEnsemble = None - - self.data_surfaces: list[Surface2DAbstract] = None - - self.data_geometry_general: cdc.CalculationDataGeometryGeneral = None - self.data_image_processing_general: cdc.CalculationImageProcessingGeneral = None - self.data_geometry_facet: list[cdc.CalculationDataGeometryFacet] = None - self.data_image_processing_facet: list[cdc.CalculationImageProcessingFacet] = None - self.data_error: cdc.CalculationError = None - - self.data_characterization_facet: list[SlopeSolverData] = None - self.data_characterization_ensemble: list[cdc.CalculationFacetEnsemble] = None + self.params: ParamsSofastFringe = ParamsSofastFringe() def help(self) -> None: """Prints Sofast doc string""" @@ -258,9 +229,8 @@ def process_optic_multifacet( List of facet data objects. ensemble_data : DefinitionEnsemble Ensemble data object. - surface_data : dict - See Sofast documentation or Sofast.help() for more details. - + surface_data : list[Surface2dAbstract] + List of surface type definitions. """ # Check inputs if len(facet_data) != len(surfaces): @@ -515,7 +485,7 @@ def _solve_slopes(self, surfaces: list[Surface2DAbstract]) -> None: lt.error_and_raise(ValueError, 'Not all facets geometrically processed; cannot solve slopes.') # Loop through all input facets and solve slopes - self.data_characterization_facet = [] + self.data_calculation_facet = [] for facet_idx in range(self.num_facets): # Check debug status if self.params.debug_slope_solver.debug_active: @@ -544,202 +514,7 @@ def _solve_slopes(self, surfaces: list[Surface2DAbstract]) -> None: slope_solver.solve_slopes() # Save slope data - self.data_characterization_facet.append(slope_solver.get_data()) + self.data_calculation_facet.append(slope_solver.get_data()) # Save input surface parameters data self.data_surfaces = surfaces - - def _calculate_facet_pointing(self, reference: Literal['average'] | int = 'average') -> None: - """ - Calculates facet pointing relative to the given reference. - - Parameters - ---------- - reference : 'average' | int - If 'average', the pointing reference is the average of all - facet pointing directions. If, int, that facet index is assumed - to have perfect pointing. - """ - if self.data_characterization_facet is None: - lt.error_and_raise(ValueError, 'Slopes must be solved first by running "solve_slopes".') - if reference != 'average' and not isinstance(reference, int): - lt.error_and_raise(ValueError, 'Given reference must be int or "average".') - if isinstance(reference, int) and reference >= self.num_facets: - lt.error_and_raise( - ValueError, f'Given facet index, {reference:d}, is out of range of 0-{self.num_facets - 1:d}.' - ) - - # Instantiate data list - self.data_characterization_ensemble = [] - - trans_facet_ensemble_list = [] - v_pointing_matrix = np.zeros((3, self.num_facets)) - for idx in range(self.num_facets): - # Get transformation from user-input and slope solving - trans_1 = TransformXYZ.from_R_V( - self.data_ensemble_def.r_facet_ensemble[idx], self.data_ensemble_def.v_facet_locations[idx] - ) - trans_2 = self.data_characterization_facet[idx].trans_alignment - # Calculate inverse of slope solving transform - trans_2 = TransformXYZ.from_V(-trans_2.V) * TransformXYZ.from_R(trans_2.R.inv()) - # Create local to global transformation - trans_facet_ensemble_list.append(trans_2 * trans_1) - - # Calculate pointing vector in ensemble coordinates - v_pointing = Vxyz((0, 0, 1)).rotate(trans_facet_ensemble_list[idx].R) - v_pointing_matrix[:, idx] = v_pointing.data.squeeze() - - # Calculate reference pointing direction - if isinstance(reference, int): - v_pointing_ref = Vxyz(v_pointing_matrix[:, reference]) - elif reference == 'average': - v_pointing_ref = Vxyz(v_pointing_matrix.mean(1)) - # Calculate rotation to align pointing vectors - r_align_pointing = v_pointing_ref.align_to(Vxyz((0, 0, 1))) - trans_align_pointing = TransformXYZ.from_R(r_align_pointing) - - # Apply alignment rotation to total transformation - trans_facet_ensemble_list = [trans_align_pointing * t for t in trans_facet_ensemble_list] - - # Calculate global slope and surface points - for idx in range(self.num_facets): - # Get slope data - slopes = self.data_characterization_facet[idx].slopes_facet_xy # facet coordinats - - # Calculate surface normals in local (facet) coordinates - u_surf_norms = np.ones((3, slopes.shape[1])) - u_surf_norms[:2] = -slopes - u_surf_norms = Uxyz(u_surf_norms).as_Vxyz() - - # Apply rotation to normal vectors - u_surf_norms_global = u_surf_norms.rotate(trans_facet_ensemble_list[idx].R) - # Convert normal vectors to global (ensemble) slopes - slopes_ensemble_xy = -u_surf_norms_global.data[:2] / u_surf_norms_global.data[2:] - - # Convert surface points to global (ensemble) coordinates - v_surf_points_ensemble = trans_facet_ensemble_list[idx].apply( - self.data_characterization_facet[idx].v_surf_points_facet - ) - - # Calculate pointing vectors in ensemble coordinates - v_facet_pointing_ensemble = Vxyz((0, 0, 1)).rotate(trans_facet_ensemble_list[idx].R) - - data = cdc.CalculationFacetEnsemble( - trans_facet_ensemble_list[idx], slopes_ensemble_xy, v_surf_points_ensemble, v_facet_pointing_ensemble - ) - self.data_characterization_ensemble.append(data) - - def get_optic( - self, interp_type: Literal['bilinear', 'clough_tocher', 'nearest'] = 'nearest' - ) -> FacetEnsemble | Facet: - """Returns the OpenCSP representation of the optic under test. For - single facet data collects, returns a Facet object, and for multi-facet - collects, returns a FacetEnsemble object. Each mirror is represented - by a MirrorPoint object. Each mirror origin is co-located with its - parent facet origin. - - Parameters - ---------- - interp_type : {'bilinear', 'clough_tocher', 'nearest'}, optional - Mirror interpolation type, by default 'nearest' - - Returns - ------- - FacetEnsemble if ProcessSofastFringe.optic_type = 'multi', otherwise Facet - """ - facets = [] - trans_list = [] - for idx_mirror, data in enumerate(self.data_characterization_facet): - # Get surface points - pts: Vxyz = data.v_surf_points_facet - # Get normals from slopes - slopes: np.ndarray = data.slopes_facet_xy - norm_data = np.ones((3, slopes.shape[1])) - norm_data[:2] = -slopes - norm_vecs = Uxyz(norm_data) - # Get mirror shape - if self.optic_type == 'undefined': - # Find bounding box - x1 = pts.x.min() - x2 = pts.x.max() - y1 = pts.y.min() - y2 = pts.y.max() - vertices = Vxy(([x1, x1, x2, x2], [y1, y2, y2, y1])) - else: - # Get optic region from optic definition - vertices = self.data_facet_def[idx_mirror].v_facet_corners.projXY() - shape = RegionXY.from_vertices(vertices) - # Create mirror - mirror = MirrorPoint(pts, norm_vecs, shape, interp_type) - # Create facet - facet = Facet(mirror) - # Locate facet within ensemble - if self.optic_type == 'multi': - trans: TransformXYZ = self.data_characterization_ensemble[idx_mirror].trans_facet_ensemble - trans_list.append(trans) - # Save facets - facets.append(facet) - - # Return optics - if self.optic_type == 'multi': - ensemble = FacetEnsemble(facets) - ensemble.set_facet_transform_list(trans_list) - return ensemble - else: - return facets[0] - - def save_to_hdf(self, file: str, prefix: str = ''): - """Saves data to given file. Data is stored as: PREFIX + Folder/Field_1 - - Parameters - ---------- - file : str - HDF file to save to - prefix : str - Prefix to append to folder path within HDF file (folders must be separated by "/") - """ - # Log - lt.info(f'Saving SofastFringe data to: {file:s}, in HDF5 folder: "{prefix:s}"') - - # One per measurement - if self.data_error is not None: - self.data_error.save_to_hdf(file, f'{prefix:s}DataSofastCalculation/general/') - self.data_geometry_general.save_to_hdf(file, f'{prefix:s}DataSofastCalculation/general/') - self.data_image_processing_general.save_to_hdf(file, f'{prefix:s}DataSofastCalculation/general/') - - # Sofast parameters - self.params.save_to_hdf(file, f'{prefix:s}DataSofastInput/') - - # Facet definition - if self.data_facet_def is not None: - for idx_facet, facet_data in enumerate(self.data_facet_def): - facet_data.save_to_hdf(file, f'{prefix:s}DataSofastInput/optic_definition/facet_{idx_facet:03d}/') - - # Ensemble definition - if self.data_ensemble_def is not None: - self.data_ensemble_def.save_to_hdf(file, f'{prefix:s}DataSofastInput/optic_definition/') - - # Surface definition - for idx_facet, surface in enumerate(self.data_surfaces): - surface.save_to_hdf(file, f'{prefix:s}DataSofastInput/optic_definition/facet_{idx_facet:03d}/') - - # Calculations, one per facet - for idx_facet in range(self.num_facets): - # Save facet slope data - self.data_characterization_facet[idx_facet].save_to_hdf( - file, f'{prefix:s}DataSofastCalculation/facet/facet_{idx_facet:03d}/' - ) - # Save facet geometry data - self.data_geometry_facet[idx_facet].save_to_hdf( - file, f'{prefix:s}DataSofastCalculation/facet/facet_{idx_facet:03d}/' - ) - # Save facet image processing data - self.data_image_processing_facet[idx_facet].save_to_hdf( - file, f'{prefix:s}DataSofastCalculation/facet/facet_{idx_facet:03d}/' - ) - - if self.data_characterization_ensemble: - # Save ensemle data - self.data_characterization_ensemble[idx_facet].save_to_hdf( - file, f'{prefix:s}DataSofastCalculation/facet/facet_{idx_facet:03d}/' - ) diff --git a/opencsp/app/sofast/lib/SofastConfiguration.py b/opencsp/app/sofast/lib/SofastConfiguration.py index f10798b6b..f12d7e703 100644 --- a/opencsp/app/sofast/lib/SofastConfiguration.py +++ b/opencsp/app/sofast/lib/SofastConfiguration.py @@ -62,7 +62,7 @@ def get_measurement_stats(self) -> list[dict]: for idx_facet in range(num_facets): if self._is_fringe: # Get data - data_calc = self.data_sofast_object.data_characterization_facet[idx_facet] + data_calc = self.data_sofast_object.data_calculation_facet[idx_facet] data_im_proc = self.data_sofast_object.data_image_processing_facet[idx_facet] data_surf = self.data_sofast_object.data_surfaces[idx_facet] @@ -77,10 +77,10 @@ def get_measurement_stats(self) -> list[dict]: num_samps = len(data_calc.v_surf_points_facet) else: # Get data - data_surf = self.data_sofast_object.slope_solver.surface - data_calc = self.data_sofast_object.data_slope_solver + data_surf = self.data_sofast_object.slope_solvers[idx_facet].surface + data_calc = self.data_sofast_object.data_calculation_facet[idx_facet] # Sample resolution - surf_points = self.data_sofast_object.data_slope_solver.v_surf_points_facet + surf_points = self.data_sofast_object.data_calculation_facet[idx_facet].v_surf_points_facet pts_index_xy = self.data_sofast_object.blob_index.get_data()[1] point_indices_mat = self.data_sofast_object.blob_index.get_data_mat()[1] offset_x = self.data_sofast_object.blob_index._offset_x diff --git a/opencsp/app/sofast/lib/image_processing.py b/opencsp/app/sofast/lib/image_processing.py index 942c64290..b5ad0158c 100644 --- a/opencsp/app/sofast/lib/image_processing.py +++ b/opencsp/app/sofast/lib/image_processing.py @@ -55,6 +55,11 @@ def calc_mask_raw( # Calculate histogram of delta image hist, edges = np.histogram(delta.flatten(), bins=N_BINS_IMAGE, density=True) + # Make sure first and last values of histogram are zero + hist = np.concatenate([[0], hist, [0]]) + bin_step = edges[1] - edges[0] + edges = np.concatenate([[edges[0] - bin_step], edges, [edges[-1] + bin_step]]) + # Find two peaks in histogram (light and dark regions) for dist in np.arange(N_PEAK_STEP, N_BINS_IMAGE, N_PEAK_STEP): peaks = find_peaks(x=hist, height=HIST_PEAK_THRESH, distance=dist)[0] diff --git a/opencsp/app/sofast/lib/process_optics_geometry.py b/opencsp/app/sofast/lib/process_optics_geometry.py index 76fdc27e6..fd5fc085a 100644 --- a/opencsp/app/sofast/lib/process_optics_geometry.py +++ b/opencsp/app/sofast/lib/process_optics_geometry.py @@ -359,7 +359,7 @@ def process_undefined_geometry( def process_multifacet_geometry( - facet_data: DefinitionFacet, + facet_data: list[DefinitionFacet], ensemble_data: DefinitionEnsemble, mask_raw: ndarray, v_meas_pt_ensemble: Vxyz, @@ -451,7 +451,8 @@ def process_multifacet_geometry( # Calculate ensemble corners in ensemble coordinates v_ensemble_corns_ensemble = [] - for r_facet_ensemble_cur, (idx_facet, idx_corn) in zip(r_facet_ensemble, ensemble_corns_indices): + for idx_facet, idx_corn in ensemble_corns_indices: + r_facet_ensemble_cur = r_facet_ensemble[idx_facet] v_ensemble_corns_ensemble.append( ( v_facet_locs_ensemble[idx_facet] diff --git a/opencsp/app/sofast/test/data/input/SofastConfiguration/stats_fixed.json b/opencsp/app/sofast/test/data/input/SofastConfiguration/stats_fixed.json index 325723a85..3175f0d09 100644 --- a/opencsp/app/sofast/test/data/input/SofastConfiguration/stats_fixed.json +++ b/opencsp/app/sofast/test/data/input/SofastConfiguration/stats_fixed.json @@ -3,6 +3,6 @@ "delta_x_sample_points_average": 0.019620769949538475, "delta_y_sample_points_average": 0.019562008548198912, "number_samples": 3782, - "focal_lengths_parabolic_xy": [-136.85679647500112, -146.7666251432111] + "focal_lengths_parabolic_xy": [-136.8543448388678, -146.76468150138166] } ] diff --git a/opencsp/app/sofast/test/test_ProcessSofastFixed.py b/opencsp/app/sofast/test/test_ProcessSofastFixed.py index 5166159b3..395bb782d 100644 --- a/opencsp/app/sofast/test/test_ProcessSofastFixed.py +++ b/opencsp/app/sofast/test/test_ProcessSofastFixed.py @@ -1,4 +1,4 @@ -"""Example script that runs fixed pattern deflectometry analysis on saved data +"""Test script that runs fixed pattern deflectometry analysis on saved data """ from os.path import join @@ -7,6 +7,7 @@ import matplotlib.pyplot as plt import numpy as np +from opencsp.app.sofast.lib.DefinitionEnsemble import DefinitionEnsemble from opencsp.app.sofast.lib.DefinitionFacet import DefinitionFacet from opencsp.app.sofast.lib.DotLocationsFixedPattern import DotLocationsFixedPattern from opencsp.app.sofast.lib.MeasurementSofastFixed import MeasurementSofastFixed @@ -15,6 +16,7 @@ from opencsp.app.sofast.lib.SpatialOrientation import SpatialOrientation from opencsp.common.lib.camera.Camera import Camera from opencsp.common.lib.deflectometry.Surface2DParabolic import Surface2DParabolic +from opencsp.common.lib.geometry.Vxy import Vxy from opencsp.common.lib.opencsp_path.opencsp_root_path import opencsp_code_dir import opencsp.common.lib.render.figure_management as fm import opencsp.common.lib.render_control.RenderControlAxis as rca @@ -28,25 +30,77 @@ class TestProcessSofastFixed(unittest.TestCase): @classmethod def setUpClass(cls): """Loads data and runs ProcessSofastFixed""" - dir_sofast_fixed = join(opencsp_code_dir(), 'test/data/sofast_fixed') - dir_sofast_common = join(opencsp_code_dir(), 'test/data/sofast_common') - - # Definitions - file_camera = join(dir_sofast_common, "camera_sofast.h5") - file_facet = join(dir_sofast_common, "Facet_NSTTF.json") - file_ori = join(dir_sofast_common, 'spatial_orientation.h5') - file_dot_locs = join(dir_sofast_fixed, 'data_measurement/fixed_pattern_dot_locations.h5') - file_meas = join(dir_sofast_fixed, 'data_measurement/measurement_facet.h5') - file_exp = join(dir_sofast_fixed, 'data_expected/calculation_facet.h5') - - cls.save_dir = join(opencsp_code_dir(), 'app/sofast/test/data/output/process_sofast_fixed_single_facet') + cls.dir_sofast_fixed = join(opencsp_code_dir(), 'test/data/sofast_fixed') + cls.dir_sofast_common = join(opencsp_code_dir(), 'test/data/sofast_common') + cls.save_dir = join(opencsp_code_dir(), 'app/sofast/test/data/output/process_sofast_fixed') ft.create_directories_if_necessary(cls.save_dir) lt.logger(join(cls.save_dir, 'sofast_fixed_process_log.txt'), lt.log.ERROR) - surface = Surface2DParabolic(initial_focal_lengths_xy=(150.0, 150), robust_least_squares=False, downsample=1) + cls.sofast_single_facet: ProcessSofastFixed = None + cls.sofast_facet_ensemble: ProcessSofastFixed = None + cls.exp_slopes_xy_single_facet: np.ndarray = None + cls.exp_slopes_xy_facet_ensemble: list[np.ndarray] = None + cls.exp_facet_pointing_trans: list[np.ndarray] = None + + cls._process_facet_ensemble(cls) + cls._process_single_facet(cls) + + def _process_facet_ensemble(self): + # Definitions + file_camera = join(self.dir_sofast_common, "camera_sofast.h5") + file_facet = join(self.dir_sofast_common, "Facet_split_NSTTF.json") + file_ensemble = join(self.dir_sofast_common, 'Ensemble_split_NSTTF_facet.json') + file_ori = join(self.dir_sofast_common, 'spatial_orientation.h5') + file_dot_locs = join(self.dir_sofast_fixed, 'data_measurement/fixed_pattern_dot_locations.h5') + file_meas = join(self.dir_sofast_fixed, 'data_measurement/measurement_facet_ensemble.h5') + file_exp = join(self.dir_sofast_fixed, 'data_expected/calculation_facet_ensemble.h5') + + # Load data + camera = Camera.load_from_hdf(file_camera) + ensemble_data = DefinitionEnsemble.load_from_json(file_ensemble) + facet_data = [DefinitionFacet.load_from_json(file_facet)] * ensemble_data.num_facets + dot_locs = DotLocationsFixedPattern.load_from_hdf(file_dot_locs) + orientation = SpatialOrientation.load_from_hdf(file_ori) + measurement = MeasurementSofastFixed.load_from_hdf(file_meas) + surfaces = [ + Surface2DParabolic(initial_focal_lengths_xy=(150.0, 150), robust_least_squares=False, downsample=1) + ] * ensemble_data.num_facets + + # Load expected data + self.exp_slopes_xy_facet_ensemble = [] + self.exp_facet_pointing_trans = [] + for idx_facet in range(ensemble_data.num_facets): + datasets = [ + f'DataSofastCalculation/facet/facet_{idx_facet:03d}/SlopeSolverData/slopes_facet_xy', + f'DataSofastCalculation/facet/facet_{idx_facet:03d}/CalculationEnsemble/trans_facet_ensemble', + ] + data = h5.load_hdf5_datasets(datasets, file_exp) + self.exp_slopes_xy_facet_ensemble.append(data['slopes_facet_xy']) + self.exp_facet_pointing_trans.append(data['trans_facet_ensemble']) + + # Instantiate class + params = ParamsSofastFixed.load_from_hdf(file_exp, 'DataSofastInput/') + self.sofast_facet_ensemble = ProcessSofastFixed(orientation, camera, dot_locs) + self.sofast_facet_ensemble.params = params + self.sofast_facet_ensemble.load_measurement_data(measurement) + + # Process + pts_known = Vxy(((853, 1031), (680, 683))) + xys_known = ((0, 0), (15, 0)) + self.sofast_facet_ensemble.process_multi_facet_optic(facet_data, surfaces, ensemble_data, pts_known, xys_known) + + def _process_single_facet(self): + # Definitions + file_camera = join(self.dir_sofast_common, "camera_sofast.h5") + file_facet = join(self.dir_sofast_common, "Facet_NSTTF.json") + file_ori = join(self.dir_sofast_common, 'spatial_orientation.h5') + file_dot_locs = join(self.dir_sofast_fixed, 'data_measurement/fixed_pattern_dot_locations.h5') + file_meas = join(self.dir_sofast_fixed, 'data_measurement/measurement_facet.h5') + file_exp = join(self.dir_sofast_fixed, 'data_expected/calculation_facet.h5') # Load data + surface = Surface2DParabolic(initial_focal_lengths_xy=(150.0, 150), robust_least_squares=False, downsample=1) camera = Camera.load_from_hdf(file_camera) facet_data = DefinitionFacet.load_from_json(file_facet) dot_locs = DotLocationsFixedPattern.load_from_hdf(file_dot_locs) @@ -56,37 +110,82 @@ def setUpClass(cls): # Load expected data datasets = ['DataSofastCalculation/facet/facet_000/SlopeSolverData/slopes_facet_xy'] data = h5.load_hdf5_datasets(datasets, file_exp) - cls.exp_slopes_xy = data['slopes_facet_xy'] + self.exp_slopes_xy_single_facet = data['slopes_facet_xy'] # Instantiate class params = ParamsSofastFixed.load_from_hdf(file_exp, 'DataSofastInput/') - cls.process_sofast_fixed = ProcessSofastFixed(orientation, camera, dot_locs, facet_data) - cls.process_sofast_fixed.params = params - cls.process_sofast_fixed.load_measurement_data(measurement) + self.sofast_single_facet = ProcessSofastFixed(orientation, camera, dot_locs) + self.sofast_single_facet.params = params + self.sofast_single_facet.load_measurement_data(measurement) # Process - cls.process_sofast_fixed.process_single_facet_optic(surface) + pt_known = measurement.origin + xy_known = (0, 0) + self.sofast_single_facet.process_single_facet_optic(facet_data, surface, pt_known, xy_known) - def test_save_as_hdf(self): + def test_save_facet_ensemble_as_hdf(self): """Tests saving to HDF file""" - self.process_sofast_fixed.save_to_hdf(join(self.save_dir, 'data_calculation_sofast_fixed.h5')) + self.sofast_facet_ensemble.save_to_hdf(join(self.save_dir, 'data_calculation_facet_ensemble.h5')) - def test_save_figure(self): - """Tests generating and saving a figure""" + def test_save_single_facet_as_hdf(self): + """Tests saving to HDF file""" + self.sofast_single_facet.save_to_hdf(join(self.save_dir, 'data_calculation_single_facet.h5')) + + def test_save_facet_ensemble_slope_figure(self): + """Tests genreating and saving a figure of a facet ensemble""" figure_control = rcfg.RenderControlFigure(tile_array=(1, 1), tile_square=True) axis_control_m = rca.meters() fig_mng = fm.setup_figure(figure_control, axis_control_m, title='') - mirror = self.process_sofast_fixed.get_mirror('bilinear') - mirror.plot_orthorectified_slope(res=0.002, clim=7, axis=fig_mng.axis) - fig_mng.save(self.save_dir, 'sofast_fixed_orthorectified_slope_magnitude', 'png') + ensemble = self.sofast_facet_ensemble.get_optic('bilinear') + res_meter = 0.002 # meter, 2mm resolution + clim_mrad = 7 # +/- 7 mrad color bar limits + ensemble.plot_orthorectified_slope(res=res_meter, clim=clim_mrad, axis=fig_mng.axis) + fig_mng.save(self.save_dir, 'slope_magnitude_facet_ensemble', 'png') + + def test_save_single_facet_slope_figure(self): + """Tests generating and saving a figure of a single facet""" + figure_control = rcfg.RenderControlFigure(tile_array=(1, 1), tile_square=True) + axis_control_m = rca.meters() - def test_slopes_xy(self): + fig_mng = fm.setup_figure(figure_control, axis_control_m, title='') + facet = self.sofast_single_facet.get_optic('bilinear') + res_meter = 0.002 # meter, 2mm resolution + clim_mrad = 7 # +/- 7 mrad color bar limits + facet.plot_orthorectified_slope(res=res_meter, clim=clim_mrad, axis=fig_mng.axis) + fig_mng.save(self.save_dir, 'slope_magnitude_single_facet', 'png') + + def test_slopes_xy_facet_ensemble(self): + """Tests slope data""" + for idx_facet in range(self.sofast_facet_ensemble.num_facets): + with self.subTest(i=idx_facet): + np.testing.assert_allclose( + self.sofast_facet_ensemble.data_calculation_facet[idx_facet].slopes_facet_xy, + self.exp_slopes_xy_facet_ensemble[idx_facet], + rtol=0, # relative tolerance (i.e. +/- a fixed fraction of the expected value) + atol=1e-6, # absolute tolerance (i.e. +/- a fixed value) + ) + + def test_slopes_xy_single_facet(self): """Tests slope data""" np.testing.assert_allclose( - self.process_sofast_fixed.data_slope_solver.slopes_facet_xy, self.exp_slopes_xy, rtol=0, atol=1e-6 + self.sofast_single_facet.data_calculation_facet[0].slopes_facet_xy, + self.exp_slopes_xy_single_facet, + rtol=0, # relative tolerance (i.e. +/- a fixed fraction of the expected value) + atol=1e-6, # absolute tolerance (i.e. +/- a fixed value) ) + def test_facet_pointing_ensemble(self): + """Tests facet pointing""" + for idx_facet in range(self.sofast_facet_ensemble.num_facets): + with self.subTest(i=idx_facet): + np.testing.assert_allclose( + self.sofast_facet_ensemble.data_calculation_ensemble[idx_facet].trans_facet_ensemble.matrix, + self.exp_facet_pointing_trans[idx_facet], + rtol=0, # relative tolerance (i.e. +/- a fixed fraction of the expected value) + atol=1e-6, # absolute tolerance (i.e. +/- a fixed value) + ) + def tearDown(self) -> None: # Make sure we release all matplotlib resources. plt.close('all') diff --git a/opencsp/app/sofast/test/test_SofastConfiguration.py b/opencsp/app/sofast/test/test_SofastConfiguration.py index f0b0fd300..e4cec6035 100644 --- a/opencsp/app/sofast/test/test_SofastConfiguration.py +++ b/opencsp/app/sofast/test/test_SofastConfiguration.py @@ -61,12 +61,18 @@ def _get_process_sofast_fixed(): measurement = MeasurementSofastFixed.load_from_hdf(file_meas) # Instantiate SofastFixed class and load measurement data - sofast = ProcessSofastFixed(orientation, camera, fixed_pattern_dot_locs, facet_data) + sofast = ProcessSofastFixed(orientation, camera, fixed_pattern_dot_locs) sofast.load_measurement_data(measurement) # Process + # - Since we are using an NSTTF mirror with about a 150m focal length, + # we will seed the initial focal lengths with 150 in x and y + # - Since we are testing SofastFixed (which has sparse data compared + # to SofastFringe), we will not downsample the data (set downsample to 1) surface = Surface2DParabolic(initial_focal_lengths_xy=(150.0, 150), robust_least_squares=False, downsample=1) - sofast.process_single_facet_optic(surface) + pt_known = measurement.origin + xy_known = (0, 0) + sofast.process_single_facet_optic(facet_data, surface, pt_known, xy_known) return sofast @@ -93,6 +99,10 @@ def _get_process_sofast_fringe(): facet_data = DefinitionFacet.load_from_json(file_facet) # Define surface definition (parabolic surface) + # - Since we are using an NSTTF mirror with about a 300m focal length, + # we will seed the initial focal lengths with 300 in x and y + # - Since we are testing SofastFringe (which has dense data sampling), + # we will downsample the data by a factor of 10 surface = Surface2DParabolic(initial_focal_lengths_xy=(300.0, 300.0), robust_least_squares=True, downsample=10) # Calibrate fringes diff --git a/opencsp/app/sofast/test/test_integration_multi_facet.py b/opencsp/app/sofast/test/test_integration_multi_facet.py index 0dc1c1b6c..fcdb71f35 100644 --- a/opencsp/app/sofast/test/test_integration_multi_facet.py +++ b/opencsp/app/sofast/test/test_integration_multi_facet.py @@ -132,8 +132,8 @@ def setUpClass(cls, base_dir: str | None = None): cls.file_dataset = file_dataset for idx in range(sofast.num_facets): - cls.data_test['slopes_facet_xy'].append(sofast.data_characterization_facet[idx].slopes_facet_xy) - cls.data_test['surf_coefs_facet'].append(sofast.data_characterization_facet[idx].surf_coefs_facet) + cls.data_test['slopes_facet_xy'].append(sofast.data_calculation_facet[idx].slopes_facet_xy) + cls.data_test['surf_coefs_facet'].append(sofast.data_calculation_facet[idx].surf_coefs_facet) def test_slope(self): for idx in range(self.num_facets): diff --git a/opencsp/app/sofast/test/test_integration_single_facet.py b/opencsp/app/sofast/test/test_integration_single_facet.py index 8e8c9284b..227f5796e 100644 --- a/opencsp/app/sofast/test/test_integration_single_facet.py +++ b/opencsp/app/sofast/test/test_integration_single_facet.py @@ -80,9 +80,9 @@ def setUpClass(cls, base_dir: str | None = None): sofast.process_optic_singlefacet(facet_data, surface) # Store test data - cls.slopes.append(sofast.data_characterization_facet[0].slopes_facet_xy) - cls.surf_coefs.append(sofast.data_characterization_facet[0].surf_coefs_facet) - cls.v_surf_points_facet.append(sofast.data_characterization_facet[0].v_surf_points_facet.data) + cls.slopes.append(sofast.data_calculation_facet[0].slopes_facet_xy) + cls.surf_coefs.append(sofast.data_calculation_facet[0].surf_coefs_facet) + cls.v_surf_points_facet.append(sofast.data_calculation_facet[0].v_surf_points_facet.data) def test_slopes(self): datasets = ['DataSofastCalculation/facet/facet_000/SlopeSolverData/slopes_facet_xy'] diff --git a/opencsp/app/sofast/test/test_integration_undefined.py b/opencsp/app/sofast/test/test_integration_undefined.py index 15ce0aa21..845bd118e 100644 --- a/opencsp/app/sofast/test/test_integration_undefined.py +++ b/opencsp/app/sofast/test/test_integration_undefined.py @@ -91,8 +91,8 @@ def test_undefined(self): sofast.process_optic_undefined(surface) # Test - slopes = sofast.data_characterization_facet[0].slopes_facet_xy - slope_coefs = sofast.data_characterization_facet[0].slope_coefs_facet + slopes = sofast.data_calculation_facet[0].slopes_facet_xy + slope_coefs = sofast.data_calculation_facet[0].slope_coefs_facet np.testing.assert_allclose(data['slopes_facet_xy'], slopes, atol=1e-7, rtol=0) np.testing.assert_allclose(data['slope_coefs_facet'], slope_coefs, atol=1e-8, rtol=0) diff --git a/opencsp/common/lib/csp/Facet.py b/opencsp/common/lib/csp/Facet.py index bd9753aac..437bb0984 100644 --- a/opencsp/common/lib/csp/Facet.py +++ b/opencsp/common/lib/csp/Facet.py @@ -1,25 +1,19 @@ """Facet class inherited by all facet classes""" -from typing import Callable -from warnings import warn - import numpy as np -from scipy.spatial.transform import Rotation from opencsp.common.lib.csp.MirrorAbstract import MirrorAbstract from opencsp.common.lib.csp.RayTraceable import RayTraceable from opencsp.common.lib.csp.VisualizeOrthorectifiedSlopeAbstract import VisualizeOrthorectifiedSlopeAbstract -from opencsp.common.lib.geometry.LoopXY import LoopXY from opencsp.common.lib.geometry.Pxyz import Pxyz -from opencsp.common.lib.geometry.RegionXY import RegionXY from opencsp.common.lib.geometry.RegionXY import Resolution from opencsp.common.lib.geometry.Vxy import Vxy from opencsp.common.lib.geometry.Vxyz import Vxyz from opencsp.common.lib.geometry.TransformXYZ import TransformXYZ from opencsp.common.lib.render.View3d import View3d from opencsp.common.lib.render_control.RenderControlFacet import RenderControlFacet -from opencsp.common.lib.render_control.RenderControlMirror import RenderControlMirror from opencsp.common.lib.csp.OpticOrientationAbstract import OpticOrientationAbstract +import opencsp.common.lib.tool.log_tools as lt UP = Vxyz([0, 0, 1]) @@ -45,6 +39,7 @@ def __init__(self, mirror: MirrorAbstract, name: str = None) -> 'Facet': @property def transform_mirror_to_self(self) -> TransformXYZ: + """Returns the mirror to facet 3d transform""" return self.mirror._self_to_parent_transform # override from VisualizeOrthorectifiedSlopeAbstract @@ -61,6 +56,10 @@ def axis_aligned_bounding_box(self) -> tuple[float, float, float, float]: # Get XYZ locations of all points making up mirror region points_xy = Vxy.merge([loop.vertices for loop in self.mirror.region.loops]) # mirror frame points_z = self.mirror.surface_displacement_at(points_xy) # mirror frame + # If the corners aren't in range of the mirror's interpolation function, set to 0 + if np.any(np.isnan(points_z)): + lt.warn('Could not find corner z values for facet; filling with zeros.') + points_z = np.nan_to_num(points_z, nan=0) points_xyz = Vxyz((points_xy.x, points_xy.y, points_z)) # mirror frame # Transform "mirror base" to "facet child" coordinates @@ -206,69 +205,3 @@ def draw(self, view: View3d, facet_style: RenderControlFacet = None, transform: if facet_style.draw_mirror_curvature: self.mirror.draw(view, facet_style.mirror_styles, transform * self.mirror._self_to_parent_transform) - - # pass # end function - - ### POINTING FUNCTION METHODS - # TODO TJL: Pointing Function methods are not tested with the updated base classes. - # There will need to be an addition to `Facet` that allows users to specify the ways - # a facet mounts the mirror it contains. Defining some function might - # be the way to do this, but that is a task for the future. - - def define_pointing_function_UNVERIFIED(self, func: Callable[..., TransformXYZ]) -> None: - """Sets the canting function to use. I.e., defines the - "set_pointing" function. - - Parameters - ---------- - func : Callable - Function that returns a "child to base" TransformXYZ object. - """ - self.pointing_function = func - - def set_pointing_UNVERIFIED(self, *args) -> None: - """Sets current facet canting (i.e. sets - self.ori.transform_child_to_base using the given arguments. - """ - # warn("Depricated, do not use OpticOrientation instance, use OpticOrietionAbstract.") - if self.pointing_function is None: - raise ValueError('self.pointing_function is not defined. Use self.define_pointing_function.') - # self.ori.transform_child_to_base = self.pointing_function(*args) - self.mirror._self_to_parent_transform = self.pointing_function(*args) - - @classmethod - def generate_az_el_UNVERIFIED(cls, mirror: MirrorAbstract) -> 'Facet': - """Generates Facet object defined by a simple azimuth then elevation - canting strategy. The "pointing_function" accessed by self.set_pointing - has the following inputs - - az - float - azimuth angle (rotation about z axis) in radians - - el - float - elevation angle (rotation about x axis) in radians - """ - - def pointing_function(az: float, el: float) -> TransformXYZ: - r = Rotation.from_euler('zx', [az, el], degrees=False) - return TransformXYZ.from_R(r) - - # Create facet - facet = cls(mirror) - facet.define_pointing_function(pointing_function) - - return facet - - @classmethod - def generate_rotation_defined_UNVERIFIED(cls, mirror: MirrorAbstract) -> 'Facet': - """Generates FacetCantable object defined by a given scipy Rotation object. - The "pointing_function" accessed by self.set_pointing has the following input - - rotation - scipy.spatial.transform.Rotation - rotation object - """ - - def pointing_function(rotation: Rotation) -> TransformXYZ: - return TransformXYZ.from_R(rotation) - - # Create facet - facet = cls(mirror) - facet.define_pointing_function(pointing_function) - - return facet - - ### END POINTING FUNCTION METHODS diff --git a/opencsp/common/lib/csp/FacetEnsemble.py b/opencsp/common/lib/csp/FacetEnsemble.py index d06d9d768..e28f320c5 100644 --- a/opencsp/common/lib/csp/FacetEnsemble.py +++ b/opencsp/common/lib/csp/FacetEnsemble.py @@ -14,6 +14,7 @@ from opencsp.common.lib.geometry.Vxyz import Vxyz from opencsp.common.lib.render.View3d import View3d from opencsp.common.lib.render_control.RenderControlFacetEnsemble import RenderControlFacetEnsemble +import opencsp.common.lib.tool.log_tools as lt class FacetEnsemble(RayTraceable, VisualizeOrthorectifiedSlopeAbstract, OpticOrientationAbstract): @@ -56,10 +57,15 @@ def axis_aligned_bounding_box(self) -> tuple[float, float, float, float]: """ # Get XYZ locations of all points making up mirror region xyz = [] # facet frame - for facet in self.facets: + for facet_idx, facet in enumerate(self.facets): # Get all mirror region vertices points_xy = Pxy.merge([loop.vertices for loop in facet.mirror.region.loops]) # mirror frame points_z = facet.mirror.surface_displacement_at(points_xy) # mirror frame + # If the corners aren't in range of the mirror's interpolation function, set to 0 + if np.any(np.isnan(points_z)): + lt.warn(f'Could not find corner z values for facet number {facet_idx:3d}; filling with zeros.') + points_z = np.nan_to_num(points_z, nan=0) + # Translate xyz points to their locations in the ensemble points_xyz = Pxyz((points_xy.x, points_xy.y, points_z)) # mirror frame points_xyz = facet.mirror.get_transform_relative_to(self).apply(points_xyz) # facet frame xyz.append(points_xyz) # facet frame diff --git a/opencsp/common/lib/csp/OpticOrientationAbstract.py b/opencsp/common/lib/csp/OpticOrientationAbstract.py index 2c6156804..766aa7935 100644 --- a/opencsp/common/lib/csp/OpticOrientationAbstract.py +++ b/opencsp/common/lib/csp/OpticOrientationAbstract.py @@ -1,11 +1,11 @@ -from abc import abstractmethod +from abc import abstractmethod, ABC import copy from warnings import warn from opencsp.common.lib.geometry.TransformXYZ import TransformXYZ -class OpticOrientationAbstract: +class OpticOrientationAbstract(ABC): """ Classes that extend OpticOrientationAbstract are objects that can be in different orientations, and can contain child objects (that also extend OpticOrientationAbstract) diff --git a/opencsp/test/data/sofast_common/Ensemble_split_NSTTF_facet.json b/opencsp/test/data/sofast_common/Ensemble_split_NSTTF_facet.json new file mode 100644 index 000000000..af19eb8c2 --- /dev/null +++ b/opencsp/test/data/sofast_common/Ensemble_split_NSTTF_facet.json @@ -0,0 +1,18 @@ +{ + "v_facet_locations": { + "x": [-0.320, 0.320], + "y": [ 0.000, 0.000], + "z": [ 0.000, 0.000] + }, + "r_facet_ensemble": { + "0": [0.0, 0.0, 0.0], + "1": [0.0, 0.0, 0.0] + }, + "ensemble_perimeter": { + "facet_indices": [1, 0, 0, 1], + "corner_indices": [0, 1, 2, 3] + }, + "v_centroid_ensemble": { + "x": [0.0], "y": [0.0], "z": [0.0] + } +} \ No newline at end of file diff --git a/opencsp/test/data/sofast_common/Facet_split_NSTTF.json b/opencsp/test/data/sofast_common/Facet_split_NSTTF.json new file mode 100644 index 000000000..fdf232e9e --- /dev/null +++ b/opencsp/test/data/sofast_common/Facet_split_NSTTF.json @@ -0,0 +1,10 @@ +{ + "v_facet_corners": { + "x": [ 0.299, -0.299, -0.299, 0.299], + "y": [-0.610, -0.610, 0.610, 0.610], + "z": [ 0.0, 0.0, 0.0, 0.0 ] + }, + "v_centroid_facet": { + "x": [0.0], "y": [0.0], "z": [0.0] + } +} \ No newline at end of file diff --git a/opencsp/test/data/sofast_fixed/data_expected/calculation_facet.h5 b/opencsp/test/data/sofast_fixed/data_expected/calculation_facet.h5 index 94b5f651c..40a253f06 100644 Binary files a/opencsp/test/data/sofast_fixed/data_expected/calculation_facet.h5 and b/opencsp/test/data/sofast_fixed/data_expected/calculation_facet.h5 differ diff --git a/opencsp/test/data/sofast_fixed/data_expected/calculation_facet_ensemble.h5 b/opencsp/test/data/sofast_fixed/data_expected/calculation_facet_ensemble.h5 new file mode 100644 index 000000000..193dfaeba Binary files /dev/null and b/opencsp/test/data/sofast_fixed/data_expected/calculation_facet_ensemble.h5 differ diff --git a/opencsp/test/data/sofast_fixed/data_measurement/measurement_facet_ensemble.h5 b/opencsp/test/data/sofast_fixed/data_measurement/measurement_facet_ensemble.h5 new file mode 100644 index 000000000..46bc8cbda Binary files /dev/null and b/opencsp/test/data/sofast_fixed/data_measurement/measurement_facet_ensemble.h5 differ