photogrammetry

Code for handling and processing photogrammetry data

Dataset dataclass

Container class for photogrammetry dataset. Provides a dict like interface for accessing points by their labels.

Attributes:

Name	Type	Description
data_dict	dict[str, NDArray[float32]]	Dict of photogrammetry points. You should genrally not touch this directly.

Source code in lat_alignment/photogrammetry.py

@dataclass
class Dataset:
    """
    Container class for photogrammetry dataset.
    Provides a dict like interface for accessing points by their labels.

    Attributes
    ----------
    data_dict : dict[str, NDArray[np.float32]]
        Dict of photogrammetry points.
        You should genrally not touch this directly.
    """

    data_dict: dict[str, NDArray[np.float32]]

    def _clear_cache(self):
        self.__dict__.pop("points", None)
        self.__dict__.pop("labels", None)
        self.__dict__.pop("codes", None)
        self.__dict__.pop("code_labels", None)
        self.__dict__.pop("target", None)
        self.__dict__.pop("target_labels", None)

    def __setattr__(self, name, value):
        if name == "data_dict":
            self._clear_cache()
        return super().__setattr__(name, value)

    def __setitem__(self, key, item):
        self._clear_cache()
        self.data_dict[key] = item

    def __getitem__(self, key):
        return self.data_dict[key]

    def __repr__(self):
        return repr(self.data_dict)

    def __len__(self):
        return len(self.data_dict)

    def __delitem__(self, key):
        self._clear_cache()
        del self.data_dict[key]

    def __contains__(self, item):
        return item in self.data_dict

    def __iter__(self):
        return iter(self.data_dict)

    @cached_property
    def points(self) -> NDArray[np.float32]:
        """
        Get all points in the dataset as an array.
        This is cached.
        """
        return np.array(list(self.data_dict.values()))

    @cached_property
    def labels(self) -> NDArray[np.str_]:
        """
        Get all labels in the dataset as an array.
        This is cached.
        """
        return np.array(list(self.data_dict.keys()))

    @cached_property
    def codes(self) -> NDArray[np.float32]:
        """
        Get all coded points in the dataset as an array.
        This is cached.
        """
        msk = np.char.find(self.labels, "CODE") >= 0
        return self.points[msk]

    @cached_property
    def code_labels(self) -> NDArray[np.str_]:
        """
        Get all coded labels in the dataset as an array.
        This is cached.
        """
        msk = np.char.find(self.labels, "CODE") >= 0
        return self.labels[msk]

    @cached_property
    def targets(self) -> NDArray[np.float32]:
        """
        Get all target points in the dataset as an array.
        This is cached.
        """
        msk = np.char.find(self.labels, "TARGET") >= 0
        return self.points[msk]

    @cached_property
    def target_labels(self) -> NDArray[np.str_]:
        """
        Get all target labels in the dataset as an array.
        This is cached.
        """
        msk = np.char.find(self.labels, "TARGET") >= 0
        return self.labels[msk]

    def copy(self) -> Self:
        """
        Make a deep copy of the dataset.

        Returns
        -------
        copy : Dataset
            A deep copy of this dataset.
        """
        return deepcopy(self)

code_labels cached property

Get all coded labels in the dataset as an array. This is cached.

codes cached property

Get all coded points in the dataset as an array. This is cached.

labels cached property

Get all labels in the dataset as an array. This is cached.

points cached property

Get all points in the dataset as an array. This is cached.

target_labels cached property

Get all target labels in the dataset as an array. This is cached.

targets cached property

Get all target points in the dataset as an array. This is cached.

copy()

Make a deep copy of the dataset.

Returns:

Name	Type	Description
copy	Dataset	A deep copy of this dataset.

Source code in lat_alignment/photogrammetry.py

def copy(self) -> Self:
    """
    Make a deep copy of the dataset.

    Returns
    -------
    copy : Dataset
        A deep copy of this dataset.
    """
    return deepcopy(self)

align_photo(dataset, reference, kill_refs, element='primary', scale=True, *, plot=True, max_dist=100.0)

Align photogrammetry data and then put it into mirror coordinates.

Parameters:

Name	Type	Description	Default
dataset	Dataset	The photogrammetry data to align.	required
reference	dict	Reference dictionary. Should contain a key called coords that specifies the coordinate system that the reference points are in. The rest of the keys should be optical elements (ie: "primary") pointing to a list of reference points to use. Each point given should be a tuple with two elements. The first element is a tuple with the (x, y, z) coordinates of the point in the global coordinate system. The second is a list of nearby coded targets that can be used to identify the point.	required
kill_refs	bool	If True remove reference points from the dataset.	required
element	str	The element that these points belong to. Should be either: 'primary', 'secondary', 'bearing', or 'receiver'.	'primary'
plot	bool	If True show a diagnostic plot of how well the reference points are aligned.	True
max_dist	float	Max distance in mm that the reference poing can be from the target point used to locate it.	100

Returns:

Name	Type	Description
aligned	Dataset	The photogrammetry data aligned to the reference points.
alignment	tuple[NDArray[float32], NDArray[float32]]	The transformation that aligned the points. The first element is a rotation matrix and the second is the shift.

Source code in lat_alignment/photogrammetry.py

def align_photo(
    dataset: Dataset,
    reference: dict,
    kill_refs: bool,
    element: str = "primary",
    scale: bool = True,
    *,
    plot: bool = True,
    max_dist: float = 100.0,
) -> tuple[
    Dataset,
    tuple[NDArray[np.float32], NDArray[np.float32]],
]:
    """
    Align photogrammetry data and then put it into mirror coordinates.

    Parameters
    ----------
    dataset : Dataset
        The photogrammetry data to align.
    reference : dict
        Reference dictionary.
        Should contain a key called `coords` that specifies the
        coordinate system that the reference points are in.
        The rest of the keys should be optical elements (ie: "primary")
        pointing to a list of reference points to use.
        Each point given should be a tuple with two elements.
        The first element is a tuple with the (x, y, z) coordinates
        of the point in the global coordinate system.
        The second is a list of nearby coded targets that can be used
        to identify the point.
    kill_refs : bool
        If True remove reference points from the dataset.
    element: str, default: 'primary'
        The element that these points belong to.
        Should be either: 'primary', 'secondary', 'bearing', or 'receiver'.
    plot : bool, default: True
        If True show a diagnostic plot of how well the reference points
        are aligned.
    max_dist : float, default: 100
        Max distance in mm that the reference poing can be from the target
        point used to locate it.

    Returns
    -------
    aligned : Dataset
        The photogrammetry data aligned to the reference points.
    alignment : tuple[NDArray[np.float32], NDArray[np.float32]]
        The transformation that aligned the points.
        The first element is a rotation matrix and
        the second is the shift.
    """
    logger.info("\tAligning with reference points for %s", element)
    elements = ["primary", "secondary", "bearing", "receiver"]
    if element not in elements and element != "all":
        raise ValueError(f"Invalid element: {element}")
    if len(reference) == 0:
        raise ValueError("Invalid or empty reference")
    if element not in reference and element != "all":
        raise ValueError("Element not found in reference dict")
    if "coords" not in reference:
        raise ValueError("Reference coordinate system not specified")
    if element == "primary":
        transform = partial(
            coord_transform, cfrom=reference["coords"], cto="opt_primary"
        )
    elif element == "secondary":
        transform = partial(
            coord_transform, cfrom=reference["coords"], cto="opt_secondary"
        )
    else:
        transform = partial(
            coord_transform, cfrom=reference["coords"], cto="opt_global"
        )
    if element == "all":
        all_refs = []
        for el in elements:
            if el not in reference:
                continue
            all_refs += reference[el]
        reference["all"] = all_refs

    # Lets find the points we can use
    ref = []
    pts = []
    invars = []
    for rpoint, codes in reference[element]:
        codes = np.array(codes)
        have = np.isin(codes, dataset.code_labels)
        if np.sum(have) == 0:
            continue
        coded = dataset[codes[have][0]]
        # Find the closest point
        dist = np.linalg.norm(dataset.targets - coded, axis=-1)
        if np.min(dist) > max_dist:
            continue
        label = dataset.target_labels[np.argmin(dist)]
        ref += [rpoint]
        pts += [dataset[label]]
        invars += [label]
    if len(ref) < 4:
        raise ValueError(f"Only {len(ref)} reference points found! Can't align!")
    logger.debug(
        "\t\tFound %d reference points in measurements with labels:\n\t\t\t%s",
        len(pts),
        str(invars),
    )
    pts = np.vstack(pts)
    ref = np.vstack(ref)
    pts = np.vstack((pts, np.mean(pts, 0)))
    ref = np.vstack((ref, np.mean(ref, 0)))
    ref = transform(ref)
    logger.debug("\t\tReference points in element coords:\n%s", str(ref))
    scale_fac = 1
    if scale:
        triu_idx = np.triu_indices(len(pts), 1)
        scale_fac = np.nanmedian(make_edm(ref)[triu_idx] / make_edm(pts)[triu_idx])
        pts *= scale_fac
    logger.debug("\t\tScale factor of %f applied", scale_fac)

    rot, sft = get_rigid(pts, ref, method="mean")
    pts_t = apply_transform(pts, rot, sft)

    if plot:
        fig = plt.figure()
        ax = fig.add_subplot(projection="3d")
        ax.scatter(pts[:, 0], pts[:, 1], pts[:, 2], color="g")
        ax.scatter(pts_t[:, 0], pts_t[:, 1], pts_t[:, 2], color="b")
        ax.scatter(ref[:, 0], ref[:, 1], ref[:, 2], color="r", marker="X")
        plt.show()
    logger.info(
        "\t\tRMS of reference points after alignment: %f",
        np.sqrt(np.mean((pts_t - ref) ** 2)),
    )
    coords_transformed = apply_transform(dataset.points * scale_fac, rot, sft)
    labels = dataset.labels

    if kill_refs:
        msk = ~np.isin(dataset.labels, invars)
        labels = labels[msk]
        coords_transformed = coords_transformed[msk]

    data = {label: coord for label, coord in zip(labels, coords_transformed)}
    transformed = Dataset(data)

    logger.debug("\t\tShift is %s mm", str(sft))
    logger.debug("\t\tRotation is %s deg", str(np.rad2deg(decompose_rotation(rot))))
    scale_fac = np.eye(3) * scale_fac
    rot @= scale_fac

    if plot:
        fig = plt.figure()
        ax = fig.add_subplot(projection="3d")
        ax.scatter(
            transformed.targets[:, 0],
            transformed.targets[:, 1],
            transformed.targets[:, 2],
            marker="x",
        )
        plt.show()

    return transformed, (rot, sft)