similarity

Module for computing phase similarity between complex interferogram pixels.

Uses metric from [Wang and Chen, 2022] for similarity.

Wang and Chen, 2022

Wang K. and Chen J., 2022. Accurate Persistent Scatterer Identification Based on Phase Similarity of Radar Pixels. IEEE Transactions on Geoscience and Remote Sensing. 60, pp.1--13. 10.1109/TGRS.2022.3210868

create_similarities(ifg_file_list, output_file, search_radius=7, sim_type='median', block_shape=(512, 512), num_threads=5, add_overviews=True, nearest_n=None)

Create a similarity raster from as stack of ifg files.

Parameters:

Name	Type	Description	Default
ifg_file_list	Sequence[PathOrStr]	Paths to input interferograms	required
output_file	PathOrStr	Output raster path	required
search_radius	int	Maximum radius to search for pixels, by default 7	7
sim_type	str	Type of similarity function to run, by default "median" Choices: "median", "max"	'median'
block_shape	tuple[int, int]	Size of blocks to process at one time from ifg_file_list by default (512, 512)	(512, 512)
num_threads	int	Number of parallel blocks to process, by default 5	5
add_overviews	bool	Whether to create overviews in output_file by default True	True
nearest_n	int	If provided, reform the nearest N interferograms before computing similarity.	None

Source code in src/dolphin/similarity.py

def create_similarities(
    ifg_file_list: Sequence[PathOrStr],
    output_file: PathOrStr,
    search_radius: int = 7,
    sim_type: Literal["median", "max"] = "median",
    block_shape: tuple[int, int] = (512, 512),
    num_threads: int = 5,
    add_overviews: bool = True,
    nearest_n: int | None = None,
):
    """Create a similarity raster from as stack of ifg files.

    Parameters
    ----------
    ifg_file_list : Sequence[PathOrStr]
        Paths to input interferograms
    output_file : PathOrStr
        Output raster path
    search_radius : int, optional
        Maximum radius to search for pixels, by default 7
    sim_type : str, optional
        Type of similarity function to run, by default "median"
        Choices: "median", "max"
    block_shape : tuple[int, int], optional
        Size of blocks to process at one time from `ifg_file_list`
        by default (512, 512)
    num_threads : int, optional
        Number of parallel blocks to process, by default 5
    add_overviews : bool, optional
        Whether to create overviews in `output_file` by default True
    nearest_n : int, optional
        If provided, reform the nearest N interferograms before computing similarity.

    """
    from dolphin._overviews import Resampling, create_image_overviews
    from dolphin.io import BackgroundRasterWriter, VRTStack, process_blocks
    from dolphin.timeseries import get_incidence_matrix

    if Path(output_file).exists():
        logger.info(f"{output_file} exists, skipping")
        return

    if sim_type == "median":
        sim_function = median_similarity
    elif sim_type == "max":
        sim_function = max_similarity
    else:
        raise ValueError(f"Unrecognized {sim_type = }")

    nodata_block = np.full(block_shape, fill_value=np.nan, dtype="float32")

    if nearest_n is not None:
        incidence_matrix = get_incidence_matrix(
            _create_nearest_n_pairs(len(ifg_file_list) + 1, n=nearest_n)
        )
        assert incidence_matrix.shape[1] == len(ifg_file_list)
    else:
        incidence_matrix = None

    def calc_sim(readers, rows, cols):
        block = readers[0][:, rows, cols]
        if np.sum(block) == 0 or np.isnan(block).all():
            return nodata_block[rows, cols], rows, cols

        if incidence_matrix is not None:
            block = _calc_nearest_diffs(block, incidence_matrix)

        out_avg = sim_function(ifg_stack=block, search_radius=search_radius)
        logger.debug(f"{rows = }, {cols = }, {block.shape = }, {out_avg.shape = }")
        return out_avg, rows, cols

    out_dir = Path(output_file).parent
    reader = VRTStack(ifg_file_list, outfile=out_dir / "sim_inputs.vrt")

    writer = BackgroundRasterWriter(
        output_file,
        like_filename=ifg_file_list[0],
        dtype="float32",
        driver="GTiff",
        nodata=np.nan,
    )
    process_blocks(
        [reader],
        writer,
        func=calc_sim,
        block_shape=block_shape,
        overlaps=(search_radius, search_radius),
        num_threads=num_threads,
    )
    writer.notify_finished()

    if add_overviews:
        logger.info("Creating overviews for unwrapped images")
        create_image_overviews(Path(output_file), resampling=Resampling.AVERAGE)

get_circle_idxs(max_radius, min_radius=0, sort_output=True)

Get the relative indices of neighboring pixels in a circle.

Adapted from c++ version of psps package: https://github.com/UT-Radar-Interferometry-Group/psps/blob/a15d458817fe7d06a6edaa0b3208ea78bc4782e7/src/cpp/similarity.cpp#L16

Source code in src/dolphin/similarity.py

def get_circle_idxs(
    max_radius: int, min_radius: int = 0, sort_output: bool = True
) -> np.ndarray:
    """Get the relative indices of neighboring pixels in a circle.

    Adapted from c++ version of `psps` package:
    https://github.com/UT-Radar-Interferometry-Group/psps/blob/a15d458817fe7d06a6edaa0b3208ea78bc4782e7/src/cpp/similarity.cpp#L16
    """
    # using the mid-point circle drawing algorithm to search for neighboring PS pixels
    # # code adapted from "https://www.geeksforgeeks.org/mid-point-circle-drawing-algorithm/"
    visited = np.zeros((max_radius, max_radius), dtype=bool)
    visited[0][0] = True

    indices = []
    for r in range(1, max_radius):
        x = r
        y = 0
        p = 1 - r
        if r > min_radius:
            indices.append([r, 0])
            indices.append([-r, 0])
            indices.append([0, r])
            indices.append([0, -r])

        visited[r][0] = True
        visited[0][r] = True
        # flag > 0 means there are holes between concentric circles
        flag = 0
        while x > y:
            # do not need to fill holes
            if flag == 0:
                y += 1
                if p <= 0:
                    # Mid-point is inside or on the perimeter
                    p += 2 * y + 1
                else:
                    # Mid-point is outside the perimeter
                    x -= 1
                    p += 2 * y - 2 * x + 1

            else:
                flag -= 1

            # All the perimeter points have already been visited
            if x < y:
                break

            while not visited[x - 1][y]:
                x -= 1
                flag += 1

            visited[x][y] = True
            visited[y][x] = True
            if r > min_radius:
                indices.append([x, y])
                indices.append([-x, -y])
                indices.append([x, -y])
                indices.append([-x, y])

                if x != y:
                    indices.append([y, x])
                    indices.append([-y, -x])
                    indices.append([y, -x])
                    indices.append([-y, x])

            if flag > 0:
                x += 1

    if sort_output:
        # Sorting makes it run faster, better data access patterns
        return np.array(sorted(indices))
    else:
        # Indices run from middle outward
        return np.array(indices)

max_similarity(ifg_stack, search_radius, mask=None)

Compute the maximum similarity of each pixel and its neighbors.

Resulting similarity matches Equation (6) of [Wang and Chen, 2022]

Wang and Chen, 2022

Wang K. and Chen J., 2022. Accurate Persistent Scatterer Identification Based on Phase Similarity of Radar Pixels. IEEE Transactions on Geoscience and Remote Sensing. 60, pp.1--13. 10.1109/TGRS.2022.3210868

Parameters:

Name	Type	Description	Default
ifg_stack	ArrayLike	3D stack of complex interferograms, or floating point phase. Shape is (n_ifg, rows, cols)	required
search_radius	int	maximum radius (in pixels) to search for neighbors when comparing each pixel.	required
mask	ArrayLike | None	Array of mask from True/False indicating whether to include the pixel (True) or ignore it (False).	None

Returns:

Type	Description
ndarray	2D array (shape (rows, cols)) of the maximum similarity for any neighbor at a pixel.

Source code in src/dolphin/similarity.py

def max_similarity(
    ifg_stack: ArrayLike, search_radius: int, mask: ArrayLike | None = None
):
    """Compute the maximum similarity of each pixel and its neighbors.

    Resulting similarity matches Equation (6) of [@Wang2022AccuratePersistentScatterer]

    Parameters
    ----------
    ifg_stack : ArrayLike
        3D stack of complex interferograms, or floating point phase.
        Shape is (n_ifg, rows, cols)
    search_radius: int
        maximum radius (in pixels) to search for neighbors when comparing each pixel.
    mask: ArrayLike (optional)
        Array of mask from True/False indicating whether to include the pixel (True)
        or ignore it (False).

    Returns
    -------
    np.ndarray
        2D array (shape (rows, cols)) of the maximum similarity for any neighbor
        at a pixel.

    """
    return _create_loop_and_run(
        ifg_stack=ifg_stack,
        search_radius=search_radius,
        mask=mask,
        func=np.nanmax,
    )

median_similarity(ifg_stack, search_radius, mask=None)

Compute the median similarity of each pixel and its neighbors.

Resulting similarity matches Equation (5) of [Wang and Chen, 2022]

Wang and Chen, 2022

Wang K. and Chen J., 2022. Accurate Persistent Scatterer Identification Based on Phase Similarity of Radar Pixels. IEEE Transactions on Geoscience and Remote Sensing. 60, pp.1--13. 10.1109/TGRS.2022.3210868

Parameters:

Name	Type	Description	Default
ifg_stack	ArrayLike	3D stack of complex interferograms, or floating point phase. Shape is (n_ifg, rows, cols)	required
search_radius	int	maximum radius (in pixels) to search for neighbors when comparing each pixel.	required
mask	ArrayLike | None	Array of mask from True/False indicating whether to include the pixel (True) or ignore it (False).	None

Returns:

Type	Description
ndarray	2D array (shape (rows, cols)) of the median similarity at each pixel.

Source code in src/dolphin/similarity.py

def median_similarity(
    ifg_stack: ArrayLike, search_radius: int, mask: ArrayLike | None = None
):
    """Compute the median similarity of each pixel and its neighbors.

    Resulting similarity matches Equation (5) of [@Wang2022AccuratePersistentScatterer]

    Parameters
    ----------
    ifg_stack : ArrayLike
        3D stack of complex interferograms, or floating point phase.
        Shape is (n_ifg, rows, cols)
    search_radius: int
        maximum radius (in pixels) to search for neighbors when comparing each pixel.
    mask: ArrayLike (optional)
        Array of mask from True/False indicating whether to include the pixel (True)
        or ignore it (False).

    Returns
    -------
    np.ndarray
        2D array (shape (rows, cols)) of the median similarity at each pixel.

    """
    return _create_loop_and_run(
        ifg_stack=ifg_stack,
        search_radius=search_radius,
        mask=mask,
        func=np.nanmedian,
    )

phase_similarity(x1, x2)

Compute the similarity between two complex 1D vectors.

Source code in src/dolphin/similarity.py

@numba.njit(nogil=True)
def phase_similarity(x1: ArrayLike, x2: ArrayLike):
    """Compute the similarity between two complex 1D vectors."""
    n = len(x1)
    out = 0.0
    for i in range(n):
        out += np.real(x1[i] * np.conj(x2[i]))
    return out / n