fit_dmd_affine.py

"""
Determine affine transformation mapping DMD space (object space) to camera space (image space) by projecting
an array of spots onto a relatively flat fluorescent background (e.g. a thin layer of fluorescent dye)

The affine transformation is represented by a matrix :math:`T` which relates object
space coordinates :math:`r_o = (x_o, y_o)` to image space coordinates :math:`r_i = (x_i, y_i)`,

.. math::

  \\begin{pmatrix}
  x_i\\\\
  y_i\\\\
  1
  \\end{pmatrix}
  = T
  \\begin{pmatrix}
  xo\\\\
  yo\\\\
  1
  \\end{pmatrix}

Given a function defined on object space, :math:`g(x_o, y_o)`, we can define a corresponding function on image space
:math:`\\tilde{g}(r_i) = g(T^{-1} r_i)`
"""

from typing import Optional, Union
from collections.abc import Sequence
import json
from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.figure import Figure
from matplotlib.patches import Rectangle
from matplotlib.colors import PowerNorm
from localize_psf.fit import gauss2d
from localize_psf.rois import get_centered_rois
from localize_psf.affine import xform2params, xform_points, xform_mat, fit_xform_points


def get_affine_fit_pattern(dmd_size: Sequence[int, int],
                           radii: Sequence[float] = (1., 1.5, 2.),
                           corner_size: int = 4,
                           point_spacing: int = 61,
                           mark_sep: int = 15):
    """
    Create DMD patterns of a sparse 2D grid of points all with the same radius. This is useful for determining the
    affine transformation between the DMD and the camera

    :param dmd_size: [nx, ny]
    :param radii: list of radii of spots for affine patterns.
     If more than one, more than one pattern will be generated.
    :param corner_size: size of blcosk indicating corners
    :param point_spacing: spacing between points
    :param mark_sep: separation between inversion/flip markers near center
    :return patterns, radii, centers: centers in format [[cx, cy]]
    """
    if isinstance(radii, (float, int)):
        radii = [radii]

    nx, ny = dmd_size

    # set spacing between points. Does not necessarily need to divide Nx and Ny
    xc = (point_spacing - 1) / 2
    yc = (point_spacing - 1) / 2

    cxs = np.arange(xc, nx, point_spacing)
    cys = np.arange(yc, ny, point_spacing)

    cxcx, cycy = np.meshgrid(cxs, cys)
    centers = np.concatenate((cxcx[:, :, None], cycy[:, :, None]), axis=2)

    patterns = []
    for r in radii:
        one_pt = np.zeros((point_spacing, point_spacing))
        xx, yy = np.meshgrid(range(one_pt.shape[1]), range(one_pt.shape[0]))
        rr = np.sqrt(np.square(xx - xc) + np.square(yy - yc))
        one_pt[rr < r] = 1

        mask = np.tile(one_pt, [int(np.ceil(ny / one_pt.shape[0])), int(np.ceil(nx / one_pt.shape[1]))])
        mask = mask[:ny, :nx]

        # add corners
        mask[:corner_size, :corner_size] = 1
        mask[:corner_size, -corner_size:] = 1
        mask[-corner_size:, :corner_size] = 1
        mask[-corner_size:, -corner_size:] = 1

        # add various markers to fix orientation

        # two edges
        mask[:1, :] = 1
        mask[:, :1] = 1

        # marks near center
        cx = nx // 2
        cy = ny // 2

        # block displaced along x-axis
        xstart1 = cx - mark_sep
        xend1 = xstart1 + corner_size
        ystart1 = cy - corner_size//2
        yend1 = ystart1 + corner_size
        mask[ystart1:yend1, xstart1:xend1] = 1

        # second block along x-axis
        xstart4 = cx - 2 * mark_sep
        xend4 = xstart4 + corner_size
        ystart4 = ystart1
        yend4 = yend1
        mask[ystart4:yend4, xstart4:xend4] = 1

        # central block
        xstart2 = cx - corner_size//2
        xend2 = xstart2 + corner_size
        ystart2 = cy - mark_sep
        yend2 = ystart2 + corner_size
        mask[ystart2:yend2, xstart2:xend2] = 1

        # block displaced along y-axis
        xstart3 = cx - corner_size//2
        xend3 = xstart3 + corner_size
        ystart3 = cy - corner_size//2
        yend3 = ystart3 + corner_size
        mask[ystart3:yend3, xstart3:xend3] = 1

        patterns.append(mask)

    patterns = np.asarray(patterns).astype(bool)

    return patterns, radii, centers


def fit_pattern_peaks(img: np.ndarray,
                      centers,
                      centers_init: Sequence[Sequence[float]],
                      indices_init: Sequence[Sequence[int]],
                      roi_size: int,
                      chi_squared_relative_max: float = np.inf,
                      max_position_err: float = 0.1,
                      img_sd: Optional[np.ndarray] = None,
                      debug: bool = False) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """
    Fit peaks of fluorescence image corresponding to affine calibration pattern.

    Affine calibration pattern can be obtained with dmd.get_affine_fit_pattern() function

    :param img: Fluorescence image of calibration pattern
    :param centers: spot positions on DMD
    :param centers_init: [[cy1, cx1], [cy2, cx2], ...] guess position of centers in the fluorescence image.
      Must supply at least three centers. Given one initial center, the other two should be shifted by one index
      along each direction of the DMD
    :param indices_init: indices corresponding to centers init in the calibration pattern
    :param roi_size: ROI size in pixels used in fitting
    :param chi_squared_relative_max: fits with chi squared values larger than this factor * the chi squared of the
      initial guess points will be ignored
    :param max_position_err: points where fits have larger relative position error than this value will be ignored
    :param img_sd:
    :param bool debug:
    :return:
    """

    if img_sd is None:
        img_sd = np.ones(img.shape)

    indices_init = np.array(indices_init, copy=True)

    # indices for dmd_centers from mask
    inds_a = np.arange(centers.shape[0])
    inds_b = np.arange(centers.shape[1])
    ibib, iaia = np.meshgrid(inds_b, inds_a)

    # to store fitting results
    fps = np.zeros((centers.shape[0], centers.shape[1], 7))
    chisqs = np.zeros((centers.shape[0], centers.shape[1]))

    # #############################
    # fit initial centers
    # #############################
    for ii in range(len(centers_init)):
        roi = get_centered_rois(centers_init[ii], [roi_size, roi_size])[0]

        cell = img[roi[0]:roi[1], roi[2]:roi[3]]
        cell_sd = img_sd[roi[0]:roi[1], roi[2]:roi[3]]
        xx, yy = np.meshgrid(range(roi[2], roi[3]), range(roi[0], roi[1]))

        gauss_model = gauss2d()
        result = gauss_model.fit(cell, (yy, xx), init_params=None, sd=cell_sd)
        def fit_fn(x, y): return gauss_model.model((y, x), result["fit_params"])
        pfit = result['fit_params']
        chi_sq = result['chi_squared']

        fps[indices_init[ii, 0], indices_init[ii, 1], :] = pfit
        chisqs[indices_init[ii, 0], indices_init[ii, 1]] = chi_sq

        if debug:
            dx = xx[0, 1] - xx[0, 0]
            dy = yy[1, 0] - yy[0, 0]
            extent = (xx[0, 0] - 0.5 * dx, xx[0, -1] + 0.5 * dx,
                      yy[-1, 0] + 0.5 * dy, yy[0, 0] - 0.5 * dy)

            vmin = pfit[5]
            vmax = pfit[5] + pfit[0] * 1.2

            fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2)
            ax1.imshow(cell, vmin=vmin, vmax=vmax, extent=extent)
            ax1.scatter(pfit[1], pfit[2], c='r', marker='x')
            ax2.imshow(fit_fn(xx, yy), vmin=vmin, vmax=vmax, extent=extent)
            ax3.imshow(cell - fit_fn(xx, yy), extent=(xx[0, 0], xx[0, -1], yy[-1, 0], yy[0, 0]))
            ax4.imshow(img, vmin=vmin, vmax=vmax)
            rec = Rectangle((roi[2], roi[0]), roi[3] - roi[2], roi[1] - roi[0], color='white', fill=0)
            ax4.add_artist(rec)
            fig.suptitle('(ia, ib) = (%d, %d), chi sq=%0.2f' % (indices_init[ii][0], indices_init[ii][1], chi_sq))

    # #############################
    # guess spacing between spots, i.e. vec_a, vec_b, from initial fits
    # #############################
    # iia = np.array(indices_init)
    ias1, ias2 = np.meshgrid(indices_init[:, 0], indices_init[:, 0])
    adiffs = ias1 - ias2
    ibs1, ibs2 = np.meshgrid(indices_init[:, 1], indices_init[:, 1])
    bdiffs = ibs1 - ibs2
    xdiffs = fps[ias1, ibs1, 1] - fps[ias2, ibs2, 1]
    ydiffs = fps[ias1, ibs1, 2] - fps[ias2, ibs2, 2]

    to_use_a = np.logical_and(adiffs != 0, bdiffs == 0)
    vec_a_ref = np.array([np.mean(xdiffs[to_use_a] / adiffs[to_use_a]), np.mean(ydiffs[to_use_a] / adiffs[to_use_a])])

    to_use_b = np.logical_and(adiffs == 0, bdiffs != 0)
    vec_b_ref = np.array([np.mean(xdiffs[to_use_b] / bdiffs[to_use_b]), np.mean(ydiffs[to_use_b] / bdiffs[to_use_b])])

    # set maximum chi-square value that we will consider a 'good fit'
    # the first three fits must succeed
    chi_sq_max = np.nanmax(chisqs[chisqs != 0]) * chi_squared_relative_max

    # loop over points
    # estimate position by guessing vec_a and vec_b from nearest pairs
    while np.sum(chisqs != 0) < chisqs.size:
        with np.errstate(invalid="ignore"):
            successful_fits = np.logical_and(chisqs > 0, chisqs < chi_sq_max)
            completed_fits = (chisqs != 0)

        # find nearest point to already fitted points, only including successful fits
        ia_fitted = iaia[successful_fits]
        ib_fitted = ibib[successful_fits]
        # use broadcasting to minimize distance sum
        dists = np.asarray(
            np.sum(np.square(iaia[:, :, None] - ia_fitted) + np.square(ibib[:, :, None] - ib_fitted), axis=2),
            dtype=float)
        # exclude points already considered
        dists[completed_fits] = np.nan
        # find minimum
        ind = np.nanargmin(dists)
        ind_tuple = np.unravel_index(ind, chisqs.shape)

        # also find nearest successfully fitted point for later use
        dists = np.asarray((iaia - ind_tuple[0]) ** 2 + (ibib - ind_tuple[1]) ** 2, dtype=float)
        dists[np.logical_not(successful_fits)] = np.nan
        nearest_ind = np.nanargmin(dists)
        nearest_ind_tuple = np.unravel_index(nearest_ind, chisqs.shape)

        # find nearest fitted pairs
        pair_fitted_a = np.logical_and(successful_fits[1:, :], successful_fits[:-1, :])
        # get coordinates for pairs
        bb, aa = np.meshgrid(inds_b, 0.5 * (inds_a[1:] + inds_a[:-1]))
        # find pair with minimum distance
        dists = (aa - ind_tuple[0]) ** 2 + (bb - ind_tuple[1]) ** 2
        dists[pair_fitted_a != 1] = np.nan
        ind_pair_a = np.nanargmin(dists)
        ind_pair_tuple_a = np.unravel_index(ind_pair_a, aa.shape)
        amin = aa[ind_pair_tuple_a]
        a1 = int(np.floor(amin))
        a2 = int(np.ceil(amin))
        bmin = bb[ind_pair_tuple_a]
        # get estimated a vector
        vec_a_guess = fps[a2, bmin, 1:3] - fps[a1, bmin, 1:3]

        # similarly for b
        pair_fitted_b = np.logical_and(successful_fits[:, 1:], successful_fits[:, :-1])
        bb, aa = np.meshgrid(0.5 * (inds_b[1:] + inds_b[:-1]), inds_a)
        # find pair with minimum distance
        dists = (aa - ind_tuple[0]) ** 2 + (bb - ind_tuple[1]) ** 2
        dists[pair_fitted_b != 1] = np.nan
        ind_pair_b = np.nanargmin(dists)
        ind_pair_tuple_b = np.unravel_index(ind_pair_b, aa.shape)
        amin = aa[ind_pair_tuple_b]
        bmin = bb[ind_pair_tuple_b]
        b1 = int(np.floor(bmin))
        b2 = int(np.ceil(bmin))
        # get estimated b vector
        vec_b_guess = fps[amin, b2, 1:3] - fps[amin, b1, 1:3]

        # guess point
        diff_a = ind_tuple[0] - nearest_ind_tuple[0]
        diff_b = ind_tuple[1] - nearest_ind_tuple[1]
        center_guess = fps[nearest_ind_tuple[0], nearest_ind_tuple[1], 1:3] + \
                       vec_a_guess * diff_a + vec_b_guess * diff_b
        xc = int(center_guess[0])
        yc = int(center_guess[1])

        # get roi
        roi = get_centered_rois([yc, xc], [roi_size, roi_size])[0]

        xstart = int(roi[2])
        xend = int(roi[3])
        ystart = int(roi[0])
        yend = int(roi[1])

        # do fitting if end points are reasonable
        end_points_ok = xstart >= 0 and xend < img.shape[1] and ystart >= 0 and yend < img.shape[0]
        if end_points_ok:
            xx, yy = np.meshgrid(range(xstart, xend), range(ystart, yend))
            cell = img[ystart:yend, xstart:xend]
            cell_sd = img_sd[ystart:yend, xstart:xend]

            gauss_model = gauss2d()
            # init_params = gauss_model.estimate_parameters(cell, (yy, xx))
            result = gauss_model.fit(cell, (yy, xx), init_params=None, sd=cell_sd)
            def fit_fn(x, y): return gauss_model.model((y, x), result["fit_params"])
            pfit = result['fit_params']
            chi_sq = result['chi_squared']

            # check if actual center was within max_position_err of the distance we expected
            vec_ref = diff_a * vec_a_ref + diff_b * vec_b_ref
            center_err = np.linalg.norm(np.array([pfit[1], pfit[2]]) - np.array([xc, yc])) / np.linalg.norm(vec_ref)

            if center_err > max_position_err:
                fps[ind_tuple[0], ind_tuple[1], :] = np.nan
                chisqs[ind_tuple] = np.nan

            else:
                fps[ind_tuple[0], ind_tuple[1], :] = pfit
                chisqs[ind_tuple] = chi_sq

                if debug and np.sum(chisqs != 0) < 20:
                    vmin = pfit[5]
                    vmax = pfit[5] + pfit[0] * 1.2

                    fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2)
                    ax1.imshow(cell, vmin=vmin, vmax=vmax, extent=(xx[0, 0], xx[0, -1], yy[-1, 0], yy[0, 0]))
                    ax1.scatter(pfit[1], pfit[2], c='r', marker='x')
                    ax2.imshow(fit_fn(xx, yy), vmin=vmin, vmax=vmax, extent=(xx[0, 0], xx[0, -1], yy[-1, 0], yy[0, 0]))
                    ax3.imshow(cell - fit_fn(xx, yy), extent=(xx[0, 0], xx[0, -1], yy[-1, 0], yy[0, 0]))
                    ax4.imshow(img, vmin=vmin, vmax=vmax)
                    rec = Rectangle((xstart, ystart), xend - xstart, yend - ystart, color='white', fill=0)
                    ax4.add_artist(rec)
                    fig.suptitle('(ia, ib) = (%d, %d), chi sq=%0.2f' % (ind_tuple[0], ind_tuple[1], chi_sq))

        else:
            fps[ind_tuple[0], ind_tuple[1], :] = np.nan
            chisqs[ind_tuple] = np.nan

    # final succesful fits
    with np.errstate(invalid="ignore"):
        successful_fits = np.logical_and(chisqs > 0, chisqs < chi_sq_max)

    return fps, chisqs, successful_fits


def plot_affine_summary(img: np.ndarray,
                        mask: np.ndarray,
                        fps: np.ndarray,
                        chisqs: np.ndarray,
                        succesful_fits: np.ndarray,
                        dmd_centers: np.ndarray,
                        affine_xform: np.ndarray,
                        options: dict,
                        indices_init: Optional[Sequence[Sequence[int]]] = None,
                        vmin_percentile: float = 5,
                        vmax_percentile: float = 99.9,
                        gamma: float = 1.,
                        cmap: str = "inferno",
                        **kwargs) -> Figure:
    """
    Plot results of DMD affine transformation fitting using results of fit_pattern_peaks()

    :param img:
    :param mask:
    :param fps: fit parameters for each peak ny x nx x 7 array
    :param chisqs: chi squared statistic for each peak
    :param bool succesful_fits: boolean giving which fits we consider successful
    :param dmd_centers: ny x nx x 2, center positiosn on DMD
    :param affine_xform: affine transformation
    :param options: {'cam_pix', 'dmd_pix', 'dmd2cam_mag_expected', 'cam_mag'}
    :param indices_init:
    :param vmin_percentile:
    :param vmax_percentile:
    :param gamma: gamma used in displaying image
    :param cmap:
    :param kwargs: passed through to figure
    :return fig: figure handle to summary figure
    """

    if "cam_mag" not in options.keys():
        options.update({"cam_mag": np.nan})

    fps = np.array(fps, copy=True)
    # ensure longer sigma is first
    with np.errstate(invalid="ignore"):
        to_swap = fps[:, :, 4] > fps[:, :, 3]
    sigma_max = np.array(fps[to_swap, 4], copy=True)
    sigma_min = np.array(fps[to_swap, 3], copy=True)

    fps[to_swap, 3] = sigma_max
    fps[to_swap, 4] = sigma_min
    fps[to_swap, 6] += np.pi / 2

    # extract useful info from fits
    xcam = fps[succesful_fits, 1]
    ycam = fps[succesful_fits, 2]
    amps = fps[:, :, 0]
    bgs = fps[:, :, 5]
    with np.errstate(invalid="ignore"):
        sigma_mean_cam = np.sqrt(fps[:, :, 3] * fps[:, :, 4])
        sigma_asymmetry_cam = np.abs(fps[:, :, 3] - fps[:, :, 4]) / (0.5 * (fps[:, :, 3] + fps[:, :, 4]))
        # angles differing by pi represent same ellipse
        angles = np.mod(fps[:, :, 6], np.pi)

    # convert fits to DMD space
    xform_inv = np.linalg.inv(affine_xform)
    centers_fit_xform = xform_points(np.stack((xcam.ravel(), ycam.ravel()), axis=1), xform_inv)
    xcam_xform = centers_fit_xform[:, 0]
    ycam_xform = centers_fit_xform[:, 1]

    # extract parameters from options
    pixel_correction_factor = options['cam_pix'] / options['dmd_pix']
    expected_mag = options['dmd2cam_mag_expected']

    # transform DMD points and mask to image space
    dmd_coords_xform = xform_points(dmd_centers[succesful_fits, :], affine_xform)
    xdmd_xform = dmd_coords_xform[:, 0]
    ydmd_xform = dmd_coords_xform[:, 1]

    # DMD axes image space
    nvec = 100
    ny, nx = mask.shape
    xc_dmd = nx // 2
    yc_dmd = ny // 2

    dmd_axes = np.array([[xc_dmd, yc_dmd],
                         [xc_dmd + nvec, yc_dmd],
                         [xc_dmd, yc_dmd],
                         [xc_dmd, yc_dmd + nvec]])

    dmd_axes_cam = xform_points(dmd_axes, affine_xform)

    # residual position error
    residual_dist_err = np.zeros((fps.shape[0], fps.shape[1])) * np.nan
    residual_dist_err[succesful_fits] = np.sqrt((xdmd_xform - xcam)**2 + (ydmd_xform - ycam)**2)

    # get mask transformed to image space
    swap_xy = np.array([[0, 1, 0],
                        [1, 0, 0],
                        [0, 0, 1]])
    affine_xform_yx = swap_xy.dot(affine_xform.dot(swap_xy))

    xx, yy = np.meshgrid(range(img.shape[0]), range(img.shape[1]))
    mask_xformed = xform_mat(mask,
                             affine_xform_yx,
                             (yy, xx),
                             mode='nearest')

    # #####################################
    # plot results
    # #####################################
    fig = plt.figure(**kwargs)
    grid = fig.add_gridspec(nrows=6,
                            ncols=6, width_ratios=[1.5, 1.5, 1, 1, 1, 1])

    # chi squareds
    ax = fig.add_subplot(grid[:2, 0])
    ax.set_title('$\\chi^2$')
    no_nans = chisqs.ravel()[np.logical_not(np.isnan(chisqs.ravel()))]
    vmin = np.percentile(no_nans, 1)
    vmax = np.percentile(no_nans, 90)
    im = ax.imshow(chisqs, vmin=vmin, vmax=vmax, cmap=cmap)
    plt.colorbar(im)

    # position errors
    ax = fig.add_subplot(grid[:2, 1])
    ax.set_title('position error (pix)')

    no_nans = residual_dist_err.ravel()[np.logical_not(np.isnan(residual_dist_err.ravel()))]
    vmax = np.percentile(no_nans, 90)
    im = ax.imshow(residual_dist_err, vmin=0, vmax=vmax, cmap=cmap)
    plt.colorbar(im)

    # sigmas
    ax = fig.add_subplot(grid[2:4, 0])
    no_nans = sigma_mean_cam.ravel()[np.logical_not(np.isnan(sigma_mean_cam.ravel()))]
    sigma_mean_median = np.median(no_nans)
    vmin = np.percentile(no_nans, 1)
    vmax = np.percentile(no_nans, 90)
    im = ax.imshow(sigma_mean_cam, vmin=vmin, vmax=vmax, cmap=cmap)

    sigma_m = sigma_mean_median * options["cam_pix"] / options["cam_mag"]
    ax.set_title("$\\sqrt{\\sigma_x \\sigma_y}$, "
                 f"median={sigma_mean_median:.2f} pix\n"
                 f"$\\sigma$={sigma_m * 1e9:.1f}nm, "
                 f"FWHM={sigma_m * 2 * np.sqrt(2 * np.log(2)) * 1e9:.1f}nm")
    plt.colorbar(im)

    # amplitudes
    ax = fig.add_subplot(grid[2:4, 1])
    no_nans = amps.ravel()[np.logical_not(np.isnan(amps.ravel()))]
    amp_median = np.median(no_nans)
    vmin = np.percentile(no_nans, 1)
    vmax = np.percentile(no_nans, 90)
    im = ax.imshow(amps, vmin=vmin, vmax=vmax, cmap=cmap)
    ax.set_title(f'amps median={amp_median:.0f}')
    plt.colorbar(im)

    # sigma asymmetry
    ax = fig.add_subplot(grid[4:6, 0])
    no_nans = sigma_asymmetry_cam.ravel()[np.logical_not(np.isnan(sigma_asymmetry_cam.ravel()))]
    sigma_asym_median = np.median(no_nans)
    im = ax.imshow(sigma_asymmetry_cam, vmin=0, vmax=1, cmap=cmap)
    ax.set_title(f'$\\sigma$ asym median={sigma_asym_median:.2f}')
    plt.colorbar(im)

    # angles
    ax = fig.add_subplot(grid[4:6, 1])
    no_nans = angles.ravel()[np.logical_not(np.isnan(angles.ravel()))]
    median_angle = np.median(no_nans * 180 / np.pi)
    im = ax.imshow(angles * 180 / np.pi, vmin=0, vmax=180, cmap=cmap)
    ax.set_title(f'angle, median={median_angle:.1f}$^\\deg$')
    plt.colorbar(im)

    # raw image with fit points overlayed
    ax = fig.add_subplot(grid[:3, 2:4])
    ax.set_title('raw image and fits')

    im = ax.imshow(img,
                   norm=PowerNorm(vmin=np.percentile(img, vmin_percentile),
                                  vmax=np.percentile(img, vmax_percentile),
                                  gamma=gamma),
                   cmap="bone")
    ax.plot(xcam, ycam, 'rx', label="fit points")
    ax.plot(xdmd_xform, ydmd_xform, 'y1', label="affine xform")
    ax.plot(dmd_axes_cam[:2, 0], dmd_axes_cam[:2, 1], 'm', label="dmd axes")
    ax.plot(dmd_axes_cam[2:, 0], dmd_axes_cam[2:, 1], 'm')
    ax.legend()
    ax.set_xticks([])
    ax.set_yticks([])
    plt.colorbar(im)

    # dmd image transformed
    ax = fig.add_subplot(grid[:3, 4:])
    ax.set_title('DMD pattern (camera space)')
    im = ax.imshow(mask_xformed, cmap="bone")
    ax.plot(xcam, ycam, 'rx', label="fit points")
    ax.plot(dmd_axes_cam[:2, 0], dmd_axes_cam[:2, 1], 'm')
    ax.plot(dmd_axes_cam[2:, 0], dmd_axes_cam[2:, 1], 'm')
    ax.set_xticks([])
    ax.set_yticks([])
    plt.colorbar(im)

    # dmd image
    ax = fig.add_subplot(grid[3:, 2:])
    ax.set_title("DMD pattern (DMD space)")
    im = ax.imshow(mask, cmap="bone")
    ax.plot(xcam_xform, ycam_xform, 'rx', label="fit points")
    ax.plot(dmd_axes[:2, 0], dmd_axes[:2, 1], 'm')
    ax.plot(dmd_axes[2:, 0], dmd_axes[2:, 1], 'm')
    if indices_init is not None:
        dmd_centers_init = np.array([dmd_centers[ind[0], ind[1], :] for ind in indices_init])
        ax.plot(dmd_centers_init[:, 0], dmd_centers_init[:, 1], 'gx', label="centers init")
        ax.legend()
    ax.set_xticks([])
    ax.set_yticks([])

    affine_params = xform2params(affine_xform)
    fig.suptitle(f'theta_x={affine_params[1] * 180 / np.pi:.2f}deg, '
                 f'mx={affine_params[0]:.3f}, '
                 f'cx={affine_params[2]:.1f}, '
                 f'pixel corrected mx={affine_params[0] * pixel_correction_factor:.3f}\n'
                 f'theta_y={affine_params[4] * 180 / np.pi:.2f}deg, '
                 f'my={affine_params[3]:.3f}, '
                 f'cy={affine_params[5]:.1f}, '
                 f'pixel corrected my={affine_params[3] * pixel_correction_factor:.3f}\n'
                 f'expected mag={expected_mag:.3f}'
                 )

    return fig


# main function for estimating affine transform
def estimate_xform(img: np.ndarray,
                   mask: np.ndarray,
                   pattern_centers,
                   centers_init: Sequence[Sequence[float]],
                   indices_init: Sequence[Sequence[float]],
                   options: dict,
                   roi_size: int = 25,
                   export_fname: str = "affine_xform",
                   plot: bool = True,
                   export_dir: Optional[Union[str, Path]] = None,
                   debug: bool = False,
                   vmin_percentile: float = 5.,
                   vmax_percentile: float = 99.,
                   gamma: float = 1.,
                   figsize: Sequence[float, float] = (16., 12.),
                   **kwargs) -> tuple[dict, Figure]:
    """
    Estimate affine transformation from DMD space to camera image space from an image.

    :param img: single image to analyze
    :param mask: DMD image that is being imaged
    :param pattern_centers: pattern centers in DMD coordinates
    :param centers_init: [[cy1, cx1], [cy2, cx2], [cy3, cx3]] in image space
    :param indices_init: [[ia1, ib1], [ia2, ib2], [ia3, ib3]]. ia indices are along the y-direction of the DMD pattern,
      starting with 0 which is the topmost when plotted. ib indices are long the x-direction starting with 0 which is the
      leftmost when plotted
    :param options: {'cam_pix', 'dmd_pix', 'dmd2cam_mag_expected'}. Distances are
      in meters.
    :param roi_size: size of ROI to fit each peak
    :param export_fname: file name (not including extension) to use when saving results
    :param plot:
    :param export_dir: directory to save results
    :param debug:
    :param vmin_percentile:
    :param vmax_percentile:
    :param gamma:
    :param figsize:
    :return (data, figure): affine trasnformation data and figure handle
    """

    # fit points
    fps, chisqs, succesful_fits = fit_pattern_peaks(img,
                                                    pattern_centers,
                                                    centers_init,
                                                    indices_init,
                                                    roi_size,
                                                    img_sd=None,
                                                    debug=debug,
                                                    **kwargs)

    # affine_xform = get_affine_xform(fps, succesful_fits, pattern_centers)
    # get affine xform
    # using pixels as coordinates, instead of real distance
    out_pts = fps[succesful_fits, 1:3]
    in_pts = pattern_centers[succesful_fits, 0:2]
    affine_xform, _ = fit_xform_points(in_pts, out_pts)

    # export data
    data = {'affine_xform': affine_xform.tolist(),
            'pattern_centers': pattern_centers.tolist(),
            'fit_params': fps.tolist(),
            'fit_params_description': "[amp, cx, cy, sigma_x, sigma_y, theta, bg]",
            "chi_squareds": chisqs.tolist()}
    data.update(options)

    if export_dir is not None:
        export_dir = Path(export_dir)
        export_dir.mkdir(exist_ok=True)

        fpath = Path(export_dir) / f"{export_fname:s}.json"
        with open(fpath, "w") as f:
            json.dump(data, f, indent="\t")

    # plot data
    if not plot:
        fig = None
    else:
        fig = plot_affine_summary(img,
                                  mask,
                                  fps,
                                  chisqs,
                                  succesful_fits,
                                  pattern_centers,
                                  affine_xform,
                                  options,
                                  indices_init=indices_init,
                                  vmin_percentile=vmin_percentile,
                                  vmax_percentile=vmax_percentile,
                                  gamma=gamma,
                                  figsize=figsize)

        if export_dir is not None:
            fig.savefig(Path(export_dir) / f"{export_fname:s}.png")

        plt.show()

    return data, fig