diff --git a/opencsp/app/sofast/lib/BlobIndex.py b/opencsp/app/sofast/lib/BlobIndex.py
index 226522c16..d465873ff 100644
--- a/opencsp/app/sofast/lib/BlobIndex.py
+++ b/opencsp/app/sofast/lib/BlobIndex.py
@@ -6,6 +6,7 @@
 
 from opencsp.common.lib.geometry.Vxy import Vxy
 from opencsp.common.lib.geometry.LoopXY import LoopXY
+from opencsp.common.lib.render.lib.AbstractPlotHandler import AbstractPlotHandler
 from opencsp.common.lib.tool import log_tools as lt
 
 
@@ -19,7 +20,30 @@ class Step(Enum):
     LEFT_OR_UP = -1
 
 
-class BlobIndex:
+class DebugBlobIndex:
+    """Debug class for BlobIndex
+
+    Parameters
+    ----------
+        debug_active : bool
+            Flag to turn on debugging, default false
+        figures : list[Figure]
+            Container to hold debug figures
+        bg_image : ndarray
+            2d ndarray, by default None. Background image (i.e. the image with the blobs present)
+            to plot in debugging plots. If None, background image is not plotted.
+        name : str
+            Name of current instance to prepend to plot titles
+    """
+
+    def __init__(self):
+        self.debug_active: bool = False
+        self.figures: list = []
+        self.bg_image: np.ndarray = None
+        self.name: str = ''
+
+
+class BlobIndex(AbstractPlotHandler):
     """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
     of pixels) that have been previously found with a blob detector and attempts to assign all found
@@ -51,6 +75,8 @@ def __init__(self, points: Vxy, x_min: int, x_max: int, y_min: int, y_max: int)
         x_min/x_max/y_min/y_max : int
             Expected min/max of blob indices in x/y directions
         """
+        super().__init__()
+
         self._points = points
 
         self._num_pts = len(points)
@@ -70,14 +96,18 @@ def __init__(self, points: Vxy, x_min: int, x_max: int, y_min: int, y_max: int)
         """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.debug: DebugBlobIndex = DebugBlobIndex()
+        """BlobIndex debug object"""
+        self.shape_yx_data_mat: tuple[int, int] = (y_max - y_min + 1, x_max - x_min + 1)
+        """The yx shape of the internal data matrix that holds the point pixel locations and xy indices"""
 
         self._offset_x = -x_min  # index
         self._offset_y = -y_min  # index
         idx_x_vec = np.arange(x_min, x_max + 1)  # index
         idx_y_vec = np.arange(y_min, y_max + 1)  # index
         self._idx_x_mat, self._idx_y_mat = np.meshgrid(idx_x_vec, idx_y_vec)  # index
-        self._points_mat = np.zeros((y_max - y_min + 1, x_max - x_min + 1, 2)) * np.nan  # pixels
-        self._point_indices_mat = np.zeros((y_max - y_min + 1, x_max - x_min + 1)) * np.nan  # index
+        self._points_mat = np.zeros(self.shape_yx_data_mat + (2,)) * np.nan  # pixels
+        self._point_indices_mat = np.zeros(self.shape_yx_data_mat) * np.nan  # index
 
     def _get_assigned_point_indices(self) -> np.ndarray[int]:
         """Returns found point indices"""
@@ -238,7 +268,6 @@ def _find_nearest_in_direction(
             mask = np.logical_and(unassigned_deltas.y > 0, unassigned_deltas.y > (2 * np.abs(unassigned_deltas.x)))
             idx_x_out = idx_x
             idx_y_out = idx_y + 1
-        # Down
         elif direction == 'down':
             mask = np.logical_and(unassigned_deltas.y < 0, -unassigned_deltas.y > (2 * np.abs(unassigned_deltas.x)))
             idx_x_out = idx_x
@@ -406,10 +435,42 @@ def run(self, pt_known: Vxy, x_known: int, y_known: int) -> None:
         x/y_known : int
             XY indies of known points
         """
+        # Plot all input blobs
+        if self.debug.debug_active:
+            fig = plt.figure()
+            self._register_plot(fig)
+            if self.debug.bg_image is not None:
+                plt.imshow(self.debug.bg_image)
+            self.plot_all_points()
+            plt.title('1: All found dots in image' + self.debug.name)
+            self.debug.figures.append(fig)
+
         # Assign center point
         self._assign_center(pt_known, x_known, y_known)
+
+        # Plot assigned known center point
+        if self.debug.debug_active:
+            fig = plt.figure()
+            self._register_plot(fig)
+            if self.debug.bg_image is not None:
+                plt.imshow(self.debug.bg_image)
+            self.plot_assigned_points_labels(labels=True)
+            plt.title('2: Locations of known, center dot' + self.debug.name)
+            self.debug.figures.append(fig)
+
         # Find 3x3 core point block
         self._find_3x3_center_block(x_known, y_known)
+
+        # Plot 3x3 core block
+        if self.debug.debug_active:
+            fig = plt.figure()
+            self._register_plot(fig)
+            if self.debug.bg_image is not None:
+                plt.imshow(self.debug.bg_image)
+            self.plot_assigned_points_labels(labels=True)
+            plt.title('3: Locations of 3x3 center block' + self.debug.name)
+            self.debug.figures.append(fig)
+
         # Extend rows
         prev_num_unassigned = self._num_unassigned()
         for idx in range(self.max_num_iters):
@@ -429,8 +490,18 @@ def run(self, pt_known: Vxy, x_known: int, y_known: int) -> None:
                 break
             prev_num_unassigned = cur_num_unassigned
 
-    def plot_points_labels(self, labels: bool = False) -> None:
-        """Plots points and labels
+        # Plot all found points
+        if self.debug.debug_active:
+            fig = plt.figure()
+            self._register_plot(fig)
+            if self.debug.bg_image is not None:
+                plt.imshow(self.debug.bg_image)
+            self.plot_points_connections()
+            plt.title('4: Row assignments' + self.debug.name)
+            self.debug.figures.append(fig)
+
+    def plot_assigned_points_labels(self, labels: bool = False) -> None:
+        """Plots all assigned points with [optionally] labels
 
         Parameters
         ----------
@@ -444,6 +515,10 @@ def plot_points_labels(self, labels: bool = False) -> None:
             ):
                 plt.text(*pt.data, f'({x:.0f}, {y:.0f})')
 
+    def plot_all_points(self) -> None:
+        """Plots all points"""
+        plt.scatter(*self._points.data, color='red')
+
     def plot_points_connections(self, labels: bool = False) -> None:
         """Plots points and connections for rows/collumns
 
@@ -497,6 +572,22 @@ def get_data(self) -> tuple[Vxy, Vxy]:
         points = Vxy((x_pts[mask_assigned], y_pts[mask_assigned]))
         return points, indices
 
+    def pts_index_to_mat_index(self, pts_index: Vxy) -> tuple[np.ndarray, np.ndarray]:
+        """Returns corresponding matrix indices (see self.get_data_mat) given point
+        indices (assigned x/y indies of points.)
+
+        Parameters
+        ----------
+        pts_index : Vxy
+            Assigned point xy indices, length N.
+
+        Returns
+        -------
+        x, y
+            Length N 1d arrays of corresponding data matrix indices (see self.get_data_mat)
+        """
+        return pts_index.x - self._offset_x, pts_index.y - self._offset_y
+
     def get_data_in_region(self, loop: LoopXY) -> tuple[Vxy, Vxy]:
         """Returns found points and indices within given region
 
diff --git a/opencsp/app/sofast/lib/CalibrateSofastFixedDots.py b/opencsp/app/sofast/lib/CalibrateSofastFixedDots.py
index 395a1beb6..d8802cd50 100644
--- a/opencsp/app/sofast/lib/CalibrateSofastFixedDots.py
+++ b/opencsp/app/sofast/lib/CalibrateSofastFixedDots.py
@@ -140,7 +140,7 @@ def __init__(
         self._dot_image_points_indices: Vxy
         self._dot_image_points_indices_x: ndarray
         self._dot_image_points_indices_y: ndarray
-        self._dot_points_xyz_mat = np.ndarray((x_max - x_min + 1, y_max - y_min + 1, 3)) * np.nan
+        self._dot_points_xyz_mat = np.ndarray((y_max - y_min + 1, x_max - x_min + 1, 3)) * np.nan
         self._num_dots: int
         self._rots_cams: list[Rotation] = []
         self._vecs_cams: list[Vxyz] = []
diff --git a/opencsp/app/sofast/lib/DefinitionEnsemble.py b/opencsp/app/sofast/lib/DefinitionEnsemble.py
index bc22c9f19..af2ffd271 100644
--- a/opencsp/app/sofast/lib/DefinitionEnsemble.py
+++ b/opencsp/app/sofast/lib/DefinitionEnsemble.py
@@ -4,9 +4,11 @@
 from copy import deepcopy
 import json
 
+import matplotlib.pyplot as plt
 import numpy as np
 from scipy.spatial.transform import Rotation
 
+from opencsp.app.sofast.lib.DefinitionFacet import DefinitionFacet
 from opencsp.common.lib.geometry.Vxyz import Vxyz
 from opencsp.common.lib.tool import hdf5_tools
 
@@ -164,6 +166,14 @@ def load_from_hdf(cls, file: str, prefix: str = '') -> 'DefinitionEnsemble':
         v_centroid_ensemble = Vxyz(data['v_centroid_ensemble'])
         return cls(v_facet_locations, r_facet_ensemble, ensemble_perimeter, v_centroid_ensemble)
 
+    def plot_facet_corners_xy_proj(self, facets: list[DefinitionFacet]) -> None:
+        """Plots the xy projection of all facet corners given accompanying facet definitions"""
+        for idx_facet, (facet, R, T) in enumerate(zip(facets, self.r_facet_ensemble, self.v_facet_locations)):
+            facet: DefinitionFacet
+            corners_cur_facet = facet.v_facet_corners
+            corners_cur_ensemble: Vxyz = corners_cur_facet.rotate(R) + T
+            plt.scatter(corners_cur_ensemble.x, corners_cur_ensemble.y, label=f'Facet {idx_facet:d}')
+
 
 def _Vxyz_to_dict(V: Vxyz) -> dict:
     d = {'x': V.x.tolist(), 'y': V.y.tolist(), 'z': V.z.tolist()}
diff --git a/opencsp/app/sofast/lib/ParamsOpticGeometry.py b/opencsp/app/sofast/lib/ParamsOpticGeometry.py
index 8e48ec041..953ab00fa 100644
--- a/opencsp/app/sofast/lib/ParamsOpticGeometry.py
+++ b/opencsp/app/sofast/lib/ParamsOpticGeometry.py
@@ -8,10 +8,19 @@ class ParamsOpticGeometry(hdf5_tools.HDF5_IO_Abstract):
     """Parameter dataclass for processing optic geometry"""
 
     perimeter_refine_axial_search_dist: float = 50.0
+    """The length of the search box (along the search direction) to use when finding optic
+    perimeter. Units pixels. Default 50.0"""
     perimeter_refine_perpendicular_search_dist: float = 50.0
+    """The half-width of the search box (perpendicular to the search direction) to use when finding
+    optic perimeter. Units pixels. Default 50.0"""
     facet_corns_refine_step_length: float = 10.0
+    """The length of the search box (along the search direction) to use when refining facet corner
+    locations (when processing a facet ensemble). Units pixels. Default 10.0"""
     facet_corns_refine_perpendicular_search_dist: float = 10.0
+    """The half-width of the search box (perpendicular to the search direction) to use when
+    refining facet corner locations (when processing a facet ensemble). Units pixels. Default 10.0"""
     facet_corns_refine_frac_keep: float = 0.5
+    """The fraction of pixels to consider within search box when finding optic edges. Default 0.5"""
 
     def save_to_hdf(self, file: str, prefix: str = ''):
         """Saves data to given HDF5 file. Data is stored in PREFIX + ParamsOpticGeometry/...
diff --git a/opencsp/app/sofast/lib/ProcessSofastFixed.py b/opencsp/app/sofast/lib/ProcessSofastFixed.py
index 95ae9d58a..061b957b7 100644
--- a/opencsp/app/sofast/lib/ProcessSofastFixed.py
+++ b/opencsp/app/sofast/lib/ProcessSofastFixed.py
@@ -4,8 +4,10 @@
 import cv2 as cv
 import numpy as np
 from numpy import ndarray
+from scipy.spatial.transform import Rotation
 
-from opencsp.app.sofast.lib.BlobIndex import BlobIndex
+from opencsp.app.sofast.lib.BlobIndex import BlobIndex, DebugBlobIndex
+from opencsp.app.sofast.lib.calculation_data_classes import CalculationBlobAssignment
 from opencsp.app.sofast.lib.DefinitionEnsemble import DefinitionEnsemble
 from opencsp.app.sofast.lib.DefinitionFacet import DefinitionFacet
 from opencsp.app.sofast.lib.DotLocationsFixedPattern import DotLocationsFixedPattern
@@ -18,7 +20,9 @@
 from opencsp.common.lib.camera.Camera import Camera
 from opencsp.common.lib.deflectometry.SlopeSolver import SlopeSolver
 from opencsp.common.lib.deflectometry.Surface2DAbstract import Surface2DAbstract
+from opencsp.common.lib.geometry.LoopXY import LoopXY
 from opencsp.common.lib.geometry.Vxy import Vxy
+from opencsp.common.lib.geometry.Vxyz import Vxyz
 import opencsp.common.lib.tool.log_tools as lt
 
 
@@ -59,28 +63,33 @@ def __init__(
 
         # 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
-        """
+        self.blob_index: list[BlobIndex] = None
+        self.debug_blob_index: DebugBlobIndex = DebugBlobIndex()
+        self.data_calculation_blob_assignment: list[CalculationBlobAssignment] = None
+
+    def _find_blobs(self, pt_known: Vxy, xy_known: tuple[int, int], loop: LoopXY = None) -> BlobIndex:
+        # Finds and assigns indices to blobs in image
+        # Parameters:
+        # pt_known, Vxy, xy pixel location of known point in image, units of pixels.
+        # xy_known, tuple[int, int], length 2 tuple. Integer xy dot indices
+        # loop, LoopXY | None, only consider points inside this loop. If None, consider all points.
+        # Returns:
+        # BlobIndex class with found/assigned blobs
+
+        # Find all blobs in image
         pts_blob = ip.detect_blobs(self.measurement.image, self.blob_detector)
 
+        # Filter blobs if loop is given
+        if loop is not None:
+            mask = loop.is_inside(pts_blob)
+            pts_blob = pts_blob[mask]
+
         # 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
-        for pt_known, xy_known in zip(pts_known, xys_known):
-            blob_index.run(pt_known, xy_known[0], xy_known[1])
+        blob_index.debug = self.debug_blob_index
+        blob_index.run(pt_known, xy_known[0], xy_known[1])
 
         return blob_index
 
@@ -117,32 +126,14 @@ def load_measurement_data(self, measurement: MeasurementSofastFixed) -> None:
         """
         self.measurement = measurement
 
-    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_processing_general,
-            self.data_geometry_facet,  # list
-            self.data_image_processing_facet,  # list
-            self.data_error,
-        ) = pr.process_singlefacet_geometry(
-            self.data_facet_def[0],
-            mask_raw,
-            self.measurement.v_measure_point_facet,
-            self.measurement.dist_optic_screen,
-            self.orientation,
-            self.camera,
-            self.params.geometry,
-            self.params.debug_geometry,
-        )
+    def _process_optic_common_geometry(
+        self, rot_optic_cam: Rotation, v_cam_optic_cam: Vxyz, idx_facet: int
+    ) -> tuple[dict, CalculationBlobAssignment]:
+        # Get blob index data for current facet
+        pts_image, pts_index_xy = self.blob_index[idx_facet].get_data()
 
-        # 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 measure point unit vector
+        u_cam_measure_point_facet = self.data_geometry_facet[idx_facet].u_cam_measure_point_facet
 
         # Get screen/camera poses
         rot_cam_optic = rot_optic_cam.inv()
@@ -157,54 +148,74 @@ def _process_optic_singlefacet_geometry(self, blob_index: BlobIndex, mask_raw: n
         v_screen_points_screen = self.fixed_pattern_dot_locs.xy_indices_to_screen_coordinates(pts_index_xy)
         v_screen_points_facet = v_optic_screen_optic + v_screen_points_screen.rotate(rot_screen_optic)
 
+        # Check for nans returning from screen point calculation
+        nan_mask = np.isnan(v_screen_points_screen.data).sum(0).astype(bool)
+        active_point_mask: np.ndarray = np.logical_not(nan_mask)
+
+        # Remove nans if any present
+        if np.any(nan_mask):
+            lt.warn(
+                'ProcessSofastFixed._process_optic_common_geometry(): '
+                f'{nan_mask.sum():d} / {nan_mask.size:d} screen points are undefined '
+                f'for facet {idx_facet:d}. These data points will be removed.'
+            )
+            # Remove nan data points from screen points
+            pts_image = pts_image[active_point_mask]
+            pts_index_xy = pts_index_xy[active_point_mask]
+            v_screen_points_facet = v_screen_points_facet[active_point_mask]
+
+        # Make 2d mask of active points (w.r.t. BlobIndex internal 2d arrays)
+        x_idx_mat, y_idx_mat = self.blob_index[idx_facet].pts_index_to_mat_index(pts_index_xy)
+        active_point_mask_mat = np.zeros(self.blob_index[idx_facet].shape_yx_data_mat, dtype=bool)
+        active_point_mask_mat[y_idx_mat, x_idx_mat] = True
+
+        # Save Sofast Fixed calculation parameters (specific to Sofast Fixed calculations)
+        data_calculation_blob_assignment = CalculationBlobAssignment(pts_image, pts_index_xy, active_point_mask_mat)
+
         # 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_optic).as_Vxyz()
 
         # Update debug data
-        self.params.debug_slope_solver.optic_data = self.data_facet_def[0]
+        self.params.debug_slope_solver.optic_data = self.data_facet_def[idx_facet]
 
         # Construct surface kwargs
-        return {
+        kwargs = {
             '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_def[0].v_facet_centroid,
-            'dist_optic_screen': self.measurement.dist_optic_screen,
+            '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[0],
+            'surface': self.data_surfaces[idx_facet],
         }
 
-    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,
+        return kwargs, data_calculation_blob_assignment
+
+    def _process_optic_singlefacet_geometry(self) -> dict:
+        # 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
+
+        # Calculate optic geometry for single facet (facet index = 0)
+        kwargs, data_calculation_blob_assignment = self._process_optic_common_geometry(
+            rot_optic_cam, v_cam_optic_cam, idx_facet=0
         )
 
-        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
+        # Save blob assignment data
+        self.data_calculation_blob_assignment = [data_calculation_blob_assignment]
 
-            # Get image points and blob indices
-            pts_image, pts_index_xy = blob_index.get_data_in_region(loop)
+        return kwargs
 
+    def _process_optic_multifacet_geometry(self) -> list[dict]:
+        # Clear blob assignment calculations data container
+        self.data_calculation_blob_assignment = []
+
+        # Loop through all facets and calcualte geometry data
+        kwargs_list = []
+        for idx_facet in range(self.num_facets):
             # 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
@@ -215,42 +226,13 @@ def _process_optic_multifacet_geometry(self, blob_index: BlobIndex, mask_raw: np
             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],
-                }
+            # Calculate optic geometry for all facets
+            kwargs, data_calculation_blob_assignment = self._process_optic_common_geometry(
+                rot_facet_cam, v_cam_facet_cam, idx_facet
             )
+            kwargs_list.append(kwargs)
+            self.data_calculation_blob_assignment.append(data_calculation_blob_assignment)
+
         return kwargs_list
 
     def process_single_facet_optic(
@@ -281,14 +263,32 @@ def process_single_facet_optic(
         self.data_facet_def = [data_facet_def.copy()]
         self.data_surfaces = [surface]
 
-        # Find blobs
-        self.blob_index = self.find_blobs(pt_known, (xy_known,))
-
         # Calculate mask
         mask_raw = self._calculate_mask()
 
+        # 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_singlefacet_geometry(
+            self.data_facet_def[0],  # 0 because there is only one facet
+            mask_raw,
+            self.measurement.v_measure_point_facet,
+            self.measurement.dist_optic_screen,
+            self.orientation,
+            self.camera,
+            self.params.geometry,
+            self.params.debug_geometry,
+        )
+
+        # Find blobs
+        self.blob_index = [self._find_blobs(pt_known, xy_known)]
+
         # Generate geometry and slope solver inputs
-        kwargs = self._process_optic_singlefacet_geometry(self.blob_index, mask_raw)
+        kwargs = self._process_optic_singlefacet_geometry()
 
         # Calculate slope
         slope_solver = SlopeSolver(**kwargs)
@@ -327,7 +327,7 @@ def process_multi_facet_optic(
         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'
+                'ProcessSofastFixed: 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)}',
             )
@@ -338,14 +338,38 @@ def process_multi_facet_optic(
         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()
 
+        # 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,
+        )
+
+        # Find blobs
+        self.blob_index: list[BlobIndex] = []
+
+        for idx_facet, geom in enumerate(self.data_image_processing_facet):
+            loop = geom.loop_facet_image_refine
+            self.debug_blob_index.name = f' - Facet {idx_facet:d}'
+            self.blob_index.append(self._find_blobs(pts_known[idx_facet], xys_known[idx_facet], loop))
+
         # Generate geometry and slope solver inputs
-        kwargs_list = self._process_optic_multifacet_geometry(self.blob_index, mask_raw)
+        kwargs_list = self._process_optic_multifacet_geometry()
 
         # Calculate slope
         self.slope_solvers = []
@@ -359,3 +383,24 @@ def process_multi_facet_optic(
 
         # Calculate facet pointing
         self._calculate_facet_pointing()
+
+    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 "/")
+        """
+        # Save default data
+        super().save_to_hdf(file, prefix)
+
+        # Calculations, one per facet
+        for idx_facet in range(self.num_facets):
+            # Save facet blob index calcluations (specific to Sofast Fixed)
+            if self.data_calculation_blob_assignment is not None:
+                self.data_calculation_blob_assignment[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 f12d7e703..832cc4f18 100644
--- a/opencsp/app/sofast/lib/SofastConfiguration.py
+++ b/opencsp/app/sofast/lib/SofastConfiguration.py
@@ -51,22 +51,17 @@ def get_measurement_stats(self) -> list[dict]:
         - focal_lengths_parabolic_xy
         """
         self._check_sofast_object_loaded()
-
-        if self._is_fringe:
-            num_facets = self.data_sofast_object.num_facets
-        elif isinstance(self.data_sofast_object, ProcessSofastFixed):
-            num_facets = 1
+        num_facets = self.data_sofast_object.num_facets
 
         stats = []
-
         for idx_facet in range(num_facets):
             if self._is_fringe:
-                # Get data
+                # Get surface data
                 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]
 
-                # Sample resolution
+                # Assemble surface points in 2d arrays
                 mask = data_im_proc.mask_processed
                 im_x = np.zeros(mask.shape) * np.nan
                 im_y = np.zeros(mask.shape) * np.nan
@@ -76,19 +71,18 @@ def get_measurement_stats(self) -> list[dict]:
                 # Number of points
                 num_samps = len(data_calc.v_surf_points_facet)
             else:
-                # Get data
+                # Get surface data
                 data_surf = self.data_sofast_object.slope_solvers[idx_facet].surface
                 data_calc = self.data_sofast_object.data_calculation_facet[idx_facet]
-                # Sample resolution
+
+                # Assemble surface points in 2d arrays
                 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
-                offset_y = self.data_sofast_object.blob_index._offset_y
-                im_x = np.zeros(point_indices_mat.shape[:2]) * np.nan
-                im_y = np.zeros(point_indices_mat.shape[:2]) * np.nan
-                im_y[pts_index_xy.y - offset_y, pts_index_xy.x - offset_x] = surf_points.y
-                im_x[pts_index_xy.y - offset_y, pts_index_xy.x - offset_x] = surf_points.x
+                mask = self.data_sofast_object.data_calculation_blob_assignment[idx_facet].active_point_mask
+                im_x = mask.astype(float) * np.nan
+                im_y = mask.astype(float) * np.nan
+                im_x[mask] = surf_points.x
+                im_y[mask] = surf_points.y
+
                 # Number of points
                 num_samps = len(surf_points)
 
@@ -183,7 +177,8 @@ def visualize_setup(
             p_screen_outline = display.interp_func(Vxy(([0, 0.95, 0.95, 0, 0], [0, 0, 0.95, 0.95, 0])))
             p_screen_cent = display.interp_func(Vxy((0.5, 0.5)))
         elif self._is_fixed:
-            p_screen_outline = Vxyz([np.nan, np.nan, 0])
+            locs = self.data_sofast_object.fixed_pattern_dot_locs.xyz_dot_loc
+            p_screen_outline = Vxyz((locs[..., 0], locs[..., 1], locs[..., 2]))
             p_screen_cent = self.data_sofast_object.fixed_pattern_dot_locs.xy_indices_to_screen_coordinates(
                 Vxy([0, 0], dtype=int)
             )
@@ -239,8 +234,12 @@ def visualize_setup(
         ax.plot([x, x], [y, y], [z, z + lz2], color='black')
         ax.text(x, y, z + lz1, 'z')
 
-        # Add screen outline
-        ax.plot(*p_screen_outline.data)
+        if self._is_fixed:
+            # Add screen points
+            ax.scatter(*p_screen_outline.data, marker='.', alpha=0.5, color='blue', label='Screen Points')
+        else:
+            # Add screen outline
+            ax.plot(*p_screen_outline.data, color='red', label='Screen Outline')
 
         # Add camera position origin
         ax.scatter(*v_screen_cam_screen.data, color='black')
diff --git a/opencsp/app/sofast/lib/calculation_data_classes.py b/opencsp/app/sofast/lib/calculation_data_classes.py
index 2fb457b03..0679d8f48 100644
--- a/opencsp/app/sofast/lib/calculation_data_classes.py
+++ b/opencsp/app/sofast/lib/calculation_data_classes.py
@@ -232,6 +232,36 @@ def save_to_hdf(self, file: str, prefix: str = ''):
         _save_data_in_file(data, datasets, file)
 
 
+@dataclass
+class CalculationBlobAssignment(hdf5_tools.HDF5_SaveAbstract):
+    """Data class for holding calculated parameters from Sofast Fixed blob assignment"""
+
+    pts_image: Vxy = None
+    """Positions in the measured image that correspond to blobs (units of image pixels from upper-left corner)"""
+    pts_index_xy: Vxy = None
+    """XY indices relative to user-defined origin point (0, 0) corresponding to positions in the image (pts_image)"""
+    active_point_mask: ndarray[bool] = None
+    """2d ndarray, mask of active points"""
+
+    def save_to_hdf(self, file: str, prefix: str = '') -> None:
+        """Saves data to given HDF5 file. Data is stored in PREFIX + CalculationBlobAssignment/...
+
+        Parameters
+        ----------
+        file : str
+            HDF file to save to
+        prefix : str
+            Prefix to append to folder path within HDF file (folders must be separated by "/")
+        """
+        data = [self.pts_image.data, self.pts_index_xy.data, self.active_point_mask]
+        datasets = [
+            prefix + 'CalculationBlobAssignment/pts_image',
+            prefix + 'CalculationBlobAssignment/pts_index_xy.data',
+            prefix + 'CalculationBlobAssignment/active_point_mask',
+        ]
+        _save_data_in_file(data, datasets, file)
+
+
 @dataclass
 class CalculationFacetEnsemble(hdf5_tools.HDF5_SaveAbstract):
     """Data class used in deflectometry calculations. Holds calculations
diff --git a/opencsp/app/sofast/test/data/input/SofastConfiguration/fixed_setup_visualize.png b/opencsp/app/sofast/test/data/input/SofastConfiguration/fixed_setup_visualize.png
index 393d12a47..18c939be5 100755
Binary files a/opencsp/app/sofast/test/data/input/SofastConfiguration/fixed_setup_visualize.png and b/opencsp/app/sofast/test/data/input/SofastConfiguration/fixed_setup_visualize.png differ
diff --git a/opencsp/app/sofast/test/data/input/SofastConfiguration/fringe_setup_visualize.png b/opencsp/app/sofast/test/data/input/SofastConfiguration/fringe_setup_visualize.png
index bfab85798..da9349646 100755
Binary files a/opencsp/app/sofast/test/data/input/SofastConfiguration/fringe_setup_visualize.png and b/opencsp/app/sofast/test/data/input/SofastConfiguration/fringe_setup_visualize.png differ