steinbock.measurement

cellprofiler special

cellprofiler

create_and_save_measurement_pipeline(measurement_pipeline_file, num_channels)

Source code in steinbock/measurement/cellprofiler/cellprofiler.py

def create_and_save_measurement_pipeline(
    measurement_pipeline_file: Union[str, PathLike], num_channels: int
):
    with _measurement_pipeline_file_template.open(mode="r") as f:
        s = f.read()
    s = s.replace("{{NUM_CHANNELS}}", str(num_channels))
    with Path(measurement_pipeline_file).open(mode="w") as f:
        f.write(s)

run_object_measurement(cellprofiler_binary, measurement_pipeline_file, input_dir, output_dir, cellprofiler_plugin_dir=None)

Source code in steinbock/measurement/cellprofiler/cellprofiler.py

def run_object_measurement(
    cellprofiler_binary: str,
    measurement_pipeline_file: Union[str, PathLike],
    input_dir: Union[str, PathLike],
    output_dir: Union[str, PathLike],
    cellprofiler_plugin_dir: Union[str, PathLike, None] = None,
):
    args = [
        cellprofiler_binary,
        "-c",
        "-r",
        "-p",
        str(measurement_pipeline_file),
        "-i",
        str(input_dir),
        "-o",
        str(output_dir),
    ]
    if cellprofiler_plugin_dir is not None:
        args.append("--plugins-directory")
        args.append(str(cellprofiler_plugin_dir))
    return run_captured(args)

data

Measurement

MAX

MEAN

MEDIAN

MIN

STD

SUM

VAR

measure_intensites(img, mask, channel_names, measurement)

Source code in steinbock/measurement/data.py

def measure_intensites(
    img: np.ndarray,
    mask: np.ndarray,
    channel_names: Sequence[str],
    measurement: Measurement,
) -> pd.DataFrame:
    object_ids = np.unique(mask[mask != 0])
    data = {
        channel_name: measurement.value(img[i], labels=mask, index=object_ids)
        for i, channel_name in enumerate(channel_names)
    }
    return pd.DataFrame(
        data=data,
        index=pd.Index(object_ids, dtype=io.mask_dtype, name="Object"),
    )

measure_intensities_from_disk(img_files, mask_files, channel_names, measurement)

Source code in steinbock/measurement/data.py

def measure_intensities_from_disk(
    img_files: Sequence[Union[str, PathLike]],
    mask_files: Sequence[Union[str, PathLike]],
    channel_names: Sequence[str],
    measurement: Measurement,
) -> Generator[Tuple[Path, Path, pd.DataFrame], None, None]:
    for img_file, mask_file in zip(img_files, mask_files):
        intensities = measure_intensites(
            io.read_image(img_file),
            io.read_mask(mask_file),
            channel_names,
            measurement,
        )
        yield Path(img_file), Path(mask_file), intensities
        del intensities

measure_regionprops(img, mask, skimage_regionprops)

Source code in steinbock/measurement/data.py

def measure_regionprops(
    img: np.ndarray, mask: np.ndarray, skimage_regionprops: Sequence[str]
) -> pd.DataFrame:
    data = regionprops_table(
        mask,
        intensity_image=np.moveaxis(img, 0, -1),
        properties=skimage_regionprops,
    )
    object_ids = data.pop("label")
    return pd.DataFrame(
        data=data,
        index=pd.Index(object_ids, dtype=io.mask_dtype, name="Object"),
    )

measure_regionprops_from_disk(img_files, mask_files, skimage_regionprops)

Source code in steinbock/measurement/data.py

def measure_regionprops_from_disk(
    img_files: Sequence[Union[str, PathLike]],
    mask_files: Sequence[Union[str, PathLike]],
    skimage_regionprops: Sequence[str],
) -> Generator[Tuple[Path, Path, pd.DataFrame], None, None]:
    skimage_regionprops = list(skimage_regionprops)
    if "label" not in skimage_regionprops:
        skimage_regionprops.insert(0, "label")
    for img_file, mask_file in zip(img_files, mask_files):
        regionprops = measure_regionprops(
            io.read_image(img_file),
            io.read_mask(mask_file),
            skimage_regionprops,
        )
        yield Path(img_file), Path(mask_file), regionprops
        del regionprops

graphs

construct_centroid_dist_graph(mask, metric, dmax=None, kmax=None)

Source code in steinbock/measurement/graphs.py

def construct_centroid_dist_graph(
    mask: np.ndarray,
    metric: str,
    dmax: Optional[float] = None,
    kmax: Optional[int] = None,
) -> pd.DataFrame:
    props = regionprops(mask)
    labels = np.array([p.label for p in props])
    centroids = np.array([p.centroid for p in props])
    condensed_dists = pdist(centroids, metric=metric)
    if kmax is not None:
        k = min(kmax, len(props) - 1)
        dist_mat = squareform(condensed_dists, checks=False).astype(float)
        np.fill_diagonal(dist_mat, np.inf)
        knn_indices = np.argpartition(dist_mat, k - 1)[:, :k]
        if dmax is not None:
            knn_dists = np.take_along_axis(dist_mat, knn_indices, -1)
            indices1, indices2 = np.nonzero(knn_dists <= dmax)
            indices2 = knn_indices[(indices1, indices2)]
        else:
            indices1 = np.repeat(np.arange(len(props)), k)
            indices2 = np.ravel(knn_indices)
        distances = dist_mat[(indices1, indices2)]
    elif dmax is not None:
        (condensed_indices,) = np.nonzero(condensed_dists <= dmax)
        indices1, indices2 = _to_triu_indices(condensed_indices, len(props))
        distances = condensed_dists[condensed_indices]
        indices1, indices2, distances = (
            np.concatenate((indices1, indices2)),
            np.concatenate((indices2, indices1)),
            np.concatenate((distances, distances)),
        )
    else:
        raise ValueError("Specify either dmax or kmax (or both)")
    return pd.DataFrame(
        data={
            "Object": np.asarray(labels[indices1], dtype=io.mask_dtype),
            "Neighbor": np.asarray(labels[indices2], dtype=io.mask_dtype),
            "Distance": np.asarray(distances, dtype=float),
        }
    )

construct_euclidean_border_dist_graph(mask, dmax=None, kmax=None)

Source code in steinbock/measurement/graphs.py

def construct_euclidean_border_dist_graph(
    mask: np.ndarray,
    dmax: Optional[float] = None,
    kmax: Optional[int] = None,
) -> pd.DataFrame:
    labels1 = []
    labels2 = []
    distances = []
    unique_labels = np.unique(mask)
    unique_labels = unique_labels[unique_labels != 0]
    for label in unique_labels:
        if dmax is not None:
            ys, xs = np.nonzero(mask == label)
            dmax_int = int(np.ceil(dmax))
            ymin, ymax = np.amin(ys) - dmax_int, np.amax(ys) + dmax_int
            xmin, xmax = np.amin(xs) - dmax_int, np.amax(xs) + dmax_int
            mask_or_patch = mask[
                max(0, ymin) : min(mask.shape[0], ymax),
                max(0, xmin) : min(mask.shape[1], xmax),
            ]
            neighbor_labels = np.unique(mask_or_patch)
            neighbor_labels = neighbor_labels[neighbor_labels != 0]
        else:
            mask_or_patch = mask
            neighbor_labels = unique_labels
        dists = distance_transform_edt(mask_or_patch != label)
        neighbor_labels = neighbor_labels[neighbor_labels != label]
        neighbor_dists = np.array(
            [
                np.amin(dists[mask_or_patch == neighbor_label])
                for neighbor_label in neighbor_labels
            ]
        )
        if dmax is not None:
            neighbor_labels = neighbor_labels[neighbor_dists <= dmax]
            neighbor_dists = neighbor_dists[neighbor_dists <= dmax]
        if kmax is not None and len(neighbor_labels) > kmax:
            knn_indices = np.argpartition(neighbor_dists, kmax - 1)[:kmax]
            neighbor_labels = neighbor_labels[knn_indices]
            neighbor_dists = neighbor_dists[knn_indices]
        labels1 += [label] * len(neighbor_labels)
        labels2 += neighbor_labels.tolist()
        distances += neighbor_dists.tolist()
    return pd.DataFrame(
        data={
            "Object": np.asarray(labels1, dtype=io.mask_dtype),
            "Neighbor": np.asarray(labels2, dtype=io.mask_dtype),
            "Distance": np.asarray(distances, dtype=float),
        }
    )

construct_graph(mask, graph_type, metric=None, dmax=None, kmax=None)

Source code in steinbock/measurement/graphs.py

def construct_graph(
    mask: np.ndarray,
    graph_type: str,
    metric: Optional[str] = None,
    dmax: Optional[float] = None,
    kmax: Optional[int] = None,
) -> pd.DataFrame:
    if graph_type == "centroid":
        if metric is None:
            raise ValueError("Metric is required")
        return construct_centroid_dist_graph(
            mask, metric, dmax=dmax, kmax=kmax
        )
    if graph_type == "border":
        if metric not in (None, "euclidean"):
            raise ValueError("Metric has to be Euclidean for border distances")
        return construct_euclidean_border_dist_graph(
            mask, dmax=dmax, kmax=kmax
        )
    if graph_type == "expand":
        if metric not in (None, "euclidean"):
            raise ValueError("Metric has to be Euclidean for pixel expansion")
        if dmax is None:
            raise ValueError("Maximum distance required for pixel expansion")
        if kmax is not None:
            raise ValueError("Pixel expansion does not support kNN graphs")
        mask = expand_mask_euclidean(mask, dmax)
        return construct_euclidean_border_dist_graph(mask, dmax=1.0)
    raise ValueError(f"Unknown graph type: {graph_type}")

construct_graphs_from_disk(mask_files, graph_type, metric=None, dmax=None, kmax=None)

Source code in steinbock/measurement/graphs.py

def construct_graphs_from_disk(
    mask_files: Sequence[Union[str, PathLike]],
    graph_type: str,
    metric: Optional[str] = None,
    dmax: Optional[float] = None,
    kmax: Optional[int] = None,
) -> Generator[Tuple[Path, pd.DataFrame], None, None]:
    for mask_file in mask_files:
        mask = io.read_mask(mask_file)
        graph = construct_graph(
            mask,
            graph_type,
            metric=metric,
            dmax=dmax,
            kmax=kmax,
        )
        yield Path(mask_file), graph
        del graph

expand_mask_euclidean(mask, dmax)

Source code in steinbock/measurement/graphs.py

def expand_mask_euclidean(mask: np.ndarray, dmax: float) -> np.ndarray:
    dists, (i, j) = distance_transform_edt(mask == 0, return_indices=True)
    return np.where(dists <= dmax, mask[i, j], mask)