local_max_peak_finder reports erroneous x/y/z data with cropped datasets

[ttung]

Description

If local_max_peak_finder is called with a cropped dataset, the x/y/z pixel coordinates in the PerImageSliceSpotResults table are the pixel offsets and not the pixel coordinates.

starfish/starfish/core/spots/FindSpots/local_max_peak_finder.py

Lines 252 to 260 in 1e10dae

	res = {Axes.X.value: spot_coords[:, 2],
	Axes.Y.value: spot_coords[:, 1],
	Axes.ZPLANE.value: spot_coords[:, 0],
	Features.SPOT_RADIUS: 1,
	Features.SPOT_ID: np.arange(spot_coords.shape[0]),
	Features.INTENSITY: data_image[spot_coords[:, 0],
	spot_coords[:, 1],
	spot_coords[:, 2]],
	}

[mattcai]

This issue seems to apply to all spot finding algorithms (e.g. FindSpots.BlobDetector and FindSpots.TrackpyLocalMaxPeakFinder). One possible fix is to use transfer_physical_coords_to_intensity_table to convert pixel coordinates from cropped to original. That way the SpotFindingResults from cropped images can still be viewed on the cropped image.

The code below can be used to test the issue:

# Load STARmap
from starfish.image import ApplyTransform, LearnTransform, Filter
from starfish.types import Axes
from starfish import data, FieldOfView
from starfish.spots import FindSpots
from starfish.util.plot import imshow_plane
from starfish.core.spots.DecodeSpots.trace_builders import build_spot_traces_exact_match
from starfish import data, display
from starfish.core.intensity_table.intensity_table_coordinates import transfer_physical_coords_to_intensity_table

experiment = data.STARmap(use_test_data=True)
imgs = experiment['fov_000'].get_image('primary')

# register rounds
projection = imgs.reduce({Axes.CH, Axes.ZPLANE}, func="max")
reference_image = projection.sel({Axes.ROUND: 0})

ltt = LearnTransform.Translation(
    reference_stack=reference_image,
    axes=Axes.ROUND,
    upsampling=1000,
)
transforms = ltt.run(projection)

warp = ApplyTransform.Warp()
imgs = warp.run(
    stack=imgs,
    transforms_list=transforms,
)

# make reference round for decoding
dots = imgs.reduce({Axes.CH, Axes.ROUND}, func="max")

bd = FindSpots.BlobDetector(
    min_sigma=2,
    max_sigma=6,
    num_sigma=20,
    threshold=0.1,
    is_volume=True,
    measurement_type='mean',
)

lmp = FindSpots.LocalMaxPeakFinder(
    min_distance=2,
    stringency=8,
    min_obj_area=6,
    max_obj_area=600,
    is_volume=True
)

tlmpf = FindSpots.TrackpyLocalMaxPeakFinder(
    spot_diameter=11,  # must be odd integer
    min_mass=0.2,
    max_size=8,  # this is max radius
    separation=3,
    preprocess=False,
    percentile=80,  # this has no effect when min_mass, spot_diameter, and max_size are set properly
    verbose=True,
)

# crop imagestacks
crop_selection = {Axes.X: (300, 700), Axes.Y: (300, 700)}
cropped_imgs= imgs.sel(crop_selection)
cropped_dots = dots.sel(crop_selection)

# find spots on cropped
bd_spots = bd.run(image_stack=cropped_imgs, reference_image=cropped_dots)
lmp_spots = lmp.run(image_stack=cropped_imgs, reference_image=cropped_dots)
tlmpf_spots = tlmpf.run(image_stack=cropped_imgs, reference_image=cropped_dots)

# build spot traces
bd_table = build_spot_traces_exact_match(bd_spots)
lmp_table = build_spot_traces_exact_match(lmp_spots)
tlmpf_table = build_spot_traces_exact_match(tlmpf_spots)

# transfer physical coordinates to intensity table
new_bd_table = transfer_physical_coords_to_intensity_table(intensity_table=bd_table, spots=bd_spots)
new_lmp_table = transfer_physical_coords_to_intensity_table(intensity_table=lmp_table, spots=lmp_spots)
new_tlmpf_table = transfer_physical_coords_to_intensity_table(intensity_table=tlmpf_table, spots=tlmpf_spots)

# show x,y pixel coordinates (expected to be between 300 and 700)
print("BlobDetector pixel coordinates")
print(f"Min x: {new_bd_table.to_features_dataframe()['x'].min()}")
print(f"Max x: {new_bd_table.to_features_dataframe()['x'].max()}")
print(f"Min y: {new_bd_table.to_features_dataframe()['y'].min()}")
print(f"Max y: {new_bd_table.to_features_dataframe()['y'].max()}")

print("LocalMaxPeakFinder pixel coordinates")
print(f"Min x: {new_lmp_table.to_features_dataframe()['x'].min()}")
print(f"Max x: {new_lmp_table.to_features_dataframe()['x'].max()}")
print(f"Min y: {new_lmp_table.to_features_dataframe()['y'].min()}")
print(f"Max y: {new_lmp_table.to_features_dataframe()['y'].max()}")

print("TrackpyLocalMaxPeakFinder pixel coordinates")
print(f"Min x: {new_tlmpf_table.to_features_dataframe()['x'].min()}")
print(f"Max x: {new_tlmpf_table.to_features_dataframe()['x'].max()}")
print(f"Min y: {new_tlmpf_table.to_features_dataframe()['y'].min()}")
print(f"Max y: {new_tlmpf_table.to_features_dataframe()['y'].max()}")

# View with napari
%gui qt

# BlobDetector
# display(stack=cropped_dots, spots=new_bd_table)

# LocalMaxPeakFinder
# display(stack=cropped_dots.reduce({Axes.ZPLANE}, func="max"), spots=new_lmp_table)

# TrackpyLocalMaxPeakFinder
# display(stack=cropped_dots, spots=new_tlmpf_table)