Skip to content

Commit

Permalink
allow import of simulation functions and precompute with 1nm resolution
Browse files Browse the repository at this point in the history
  • Loading branch information
@Henley13
Henley13 committed Jun 1, 2021
1 parent e1c104e commit f544f2e
Show file tree
Hide file tree
Showing 3 changed files with 122 additions and 114 deletions.
9 changes: 8 additions & 1 deletion bigfish/detection/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,10 @@
from .spot_modeling import build_reference_spot
from .spot_modeling import modelize_spot
from .spot_modeling import precompute_erf
from .spot_modeling import initialize_grid
from .spot_modeling import gaussian_2d
from .spot_modeling import gaussian_3d


from .cluster_detection import detect_clusters

Expand All @@ -39,7 +43,10 @@
"fit_subpixel",
"build_reference_spot",
"modelize_spot",
"precompute_erf"]
"precompute_erf",
"initialize_grid",
"gaussian_2d",
"gaussian_3d"]

_clusters = [
"detect_clusters"]
Expand Down
44 changes: 22 additions & 22 deletions bigfish/detection/dense_decomposition.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,8 +13,8 @@

import bigfish.stack as stack
from .spot_modeling import build_reference_spot, modelize_spot, precompute_erf
from .spot_modeling import _gaussian_2d, _initialize_grid_2d
from .spot_modeling import _gaussian_3d, _initialize_grid_3d
from .spot_modeling import gaussian_2d, _initialize_grid_2d
from .spot_modeling import gaussian_3d, _initialize_grid_3d

from skimage.measure import regionprops
from skimage.measure import label
Expand Down Expand Up @@ -174,7 +174,7 @@ def decompose_dense(image, spots, voxel_size_z=None, voxel_size_yx=100,
return spots, dense_regions, reference_spot

# precompute gaussian function values
max_grid = max(200, region_size + 1)
max_grid = region_size + 1
precomputed_gaussian = precompute_erf(
voxel_size_z, voxel_size_yx, sigma_z, sigma_yx, max_grid=max_grid)

Expand Down Expand Up @@ -650,17 +650,17 @@ def _gaussian_mixture_3d(image, region, voxel_size_z, voxel_size_yx, sigma_z,
while diff_ssr < 0 or nb_gaussian == limit_gaussian:
position_gaussian = np.argmax(residual)
positions_gaussian.append(list(grid[:, position_gaussian]))
simulation += _gaussian_3d(grid=grid,
mu_z=float(positions_gaussian[-1][0]),
mu_y=float(positions_gaussian[-1][1]),
mu_x=float(positions_gaussian[-1][2]),
sigma_z=sigma_z,
sigma_yx=sigma_yx,
voxel_size_z=voxel_size_z,
voxel_size_yx=voxel_size_yx,
psf_amplitude=amplitude,
psf_background=background,
precomputed=precomputed_gaussian)
simulation += gaussian_3d(grid=grid,
mu_z=float(positions_gaussian[-1][0]),
mu_y=float(positions_gaussian[-1][1]),
mu_x=float(positions_gaussian[-1][2]),
sigma_z=sigma_z,
sigma_yx=sigma_yx,
voxel_size_z=voxel_size_z,
voxel_size_yx=voxel_size_yx,
psf_amplitude=amplitude,
psf_background=background,
precomputed=precomputed_gaussian)
residual = image_region_raw - simulation
new_ssr = np.sum(residual ** 2)
diff_ssr = new_ssr - ssr
Expand Down Expand Up @@ -746,14 +746,14 @@ def _gaussian_mixture_2d(image, region, voxel_size_yx, sigma_yx, amplitude,
while diff_ssr < 0 or nb_gaussian == limit_gaussian:
position_gaussian = np.argmax(residual)
positions_gaussian.append(list(grid[:, position_gaussian]))
simulation += _gaussian_2d(grid=grid,
mu_y=float(positions_gaussian[-1][0]),
mu_x=float(positions_gaussian[-1][1]),
sigma_yx=sigma_yx,
voxel_size_yx=voxel_size_yx,
psf_amplitude=amplitude,
psf_background=background,
precomputed=precomputed_gaussian)
simulation += gaussian_2d(grid=grid,
mu_y=float(positions_gaussian[-1][0]),
mu_x=float(positions_gaussian[-1][1]),
sigma_yx=sigma_yx,
voxel_size_yx=voxel_size_yx,
psf_amplitude=amplitude,
psf_background=background,
precomputed=precomputed_gaussian)
residual = image_region_raw - simulation
new_ssr = np.sum(residual ** 2)
diff_ssr = new_ssr - ssr
Expand Down
183 changes: 92 additions & 91 deletions bigfish/detection/spot_modeling.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -377,7 +377,7 @@ def modelize_spot(reference_spot, voxel_size_z=None, voxel_size_yx=100,
voxel_size_z, psf_z = None, None

# initialize a grid representing the reference spot
grid, centroid_coord = _initialize_grid(
grid, centroid_coord = initialize_grid(
image_spot=reference_spot,
voxel_size_z=voxel_size_z,
voxel_size_yx=voxel_size_yx,
Expand Down Expand Up @@ -445,8 +445,8 @@ def modelize_spot(reference_spot, voxel_size_z=None, voxel_size_yx=100,

# ### Spot modelization: initialization ###

def _initialize_grid(image_spot, voxel_size_z, voxel_size_yx,
return_centroid=False):
def initialize_grid(image_spot, voxel_size_z, voxel_size_yx,
return_centroid=False):
"""Build a grid in nanometer to compute gaussian function values over a
full volume or surface.
Expand Down Expand Up @@ -583,7 +583,6 @@ def _initialize_grid_2d(image_spot, voxel_size_yx, return_centroid=False):

# build meshgrid
yy, xx = np.meshgrid(np.arange(nb_y), np.arange(nb_x), indexing="ij")

yy = yy.astype(np.float32) * float(voxel_size_yx)
xx = xx.astype(np.float32) * float(voxel_size_yx)

Expand Down Expand Up @@ -704,50 +703,50 @@ def _objective_function_3d(voxel_size_z, voxel_size_yx, psf_z, psf_yx,
and psf_yx is not None
and psf_amplitude is None):
def f(grid, mu_z, mu_y, mu_x, psf_amplitude, psf_background):
values = _gaussian_3d(grid=grid,
mu_z=mu_z,
mu_y=mu_y,
mu_x=mu_x,
sigma_z=psf_z,
sigma_yx=psf_yx,
voxel_size_z=voxel_size_z,
voxel_size_yx=voxel_size_yx,
psf_amplitude=psf_amplitude,
psf_background=psf_background)
values = gaussian_3d(grid=grid,
mu_z=mu_z,
mu_y=mu_y,
mu_x=mu_x,
sigma_z=psf_z,
sigma_yx=psf_yx,
voxel_size_z=voxel_size_z,
voxel_size_yx=voxel_size_yx,
psf_amplitude=psf_amplitude,
psf_background=psf_background)
return values

# amplitude is known, we fit sigma, mu and background
elif (psf_amplitude is not None
and psf_z is None
and psf_yx is None):
def f(grid, mu_z, mu_y, mu_x, psf_z, psf_yx, psf_background):
values = _gaussian_3d(grid=grid,
mu_z=mu_z,
mu_y=mu_y,
mu_x=mu_x,
sigma_z=psf_z,
sigma_yx=psf_yx,
voxel_size_z=voxel_size_z,
voxel_size_yx=voxel_size_yx,
psf_amplitude=psf_amplitude,
psf_background=psf_background)
values = gaussian_3d(grid=grid,
mu_z=mu_z,
mu_y=mu_y,
mu_x=mu_x,
sigma_z=psf_z,
sigma_yx=psf_yx,
voxel_size_z=voxel_size_z,
voxel_size_yx=voxel_size_yx,
psf_amplitude=psf_amplitude,
psf_background=psf_background)
return values

# amplitude and sigma are known, we fit mu and background
elif (psf_amplitude is not None
and psf_z is not None
and psf_yx is not None):
def f(grid, mu_z, mu_y, mu_x, psf_background):
values = _gaussian_3d(grid=grid,
mu_z=mu_z,
mu_y=mu_y,
mu_x=mu_x,
sigma_z=psf_z,
sigma_yx=psf_yx,
voxel_size_z=voxel_size_z,
voxel_size_yx=voxel_size_yx,
psf_amplitude=psf_amplitude,
psf_background=psf_background)
values = gaussian_3d(grid=grid,
mu_z=mu_z,
mu_y=mu_y,
mu_x=mu_x,
sigma_z=psf_z,
sigma_yx=psf_yx,
voxel_size_z=voxel_size_z,
voxel_size_yx=voxel_size_yx,
psf_amplitude=psf_amplitude,
psf_background=psf_background)
return values

# we fit mu, sigma, amplitude and background
Expand All @@ -756,16 +755,16 @@ def f(grid, mu_z, mu_y, mu_x, psf_background):
and psf_yx is None):
def f(grid, mu_z, mu_y, mu_x, psf_z, psf_yx, psf_amplitude,
psf_background):
values = _gaussian_3d(grid=grid,
mu_z=mu_z,
mu_y=mu_y,
mu_x=mu_x,
sigma_z=psf_z,
sigma_yx=psf_yx,
voxel_size_z=voxel_size_z,
voxel_size_yx=voxel_size_yx,
psf_amplitude=psf_amplitude,
psf_background=psf_background)
values = gaussian_3d(grid=grid,
mu_z=mu_z,
mu_y=mu_y,
mu_x=mu_x,
sigma_z=psf_z,
sigma_yx=psf_yx,
voxel_size_z=voxel_size_z,
voxel_size_yx=voxel_size_yx,
psf_amplitude=psf_amplitude,
psf_background=psf_background)
return values

else:
Expand All @@ -775,11 +774,11 @@ def f(grid, mu_z, mu_y, mu_x, psf_z, psf_yx, psf_amplitude,
return f


def _gaussian_3d(grid, mu_z, mu_y, mu_x, sigma_z, sigma_yx, voxel_size_z,
voxel_size_yx, psf_amplitude, psf_background,
precomputed=None):
"""Compute the gaussian function over the grid 'xdata' representing a
volume V with shape (V_z, V_y, V_x).
def gaussian_3d(grid, mu_z, mu_y, mu_x, sigma_z, sigma_yx, voxel_size_z,
voxel_size_yx, psf_amplitude, psf_background,
precomputed=None):
"""Compute the gaussian function over the grid representing a volume V
with shape (V_z, V_y, V_x).
# TODO add equations
Expand Down Expand Up @@ -830,9 +829,9 @@ def _gaussian_3d(grid, mu_z, mu_y, mu_x, sigma_z, sigma_yx, voxel_size_z,
table_erf_x = precomputed[2]

# get indices for the tables
i_z = np.around(np.abs(meshgrid_z - mu_z) / 5).astype(np.int64)
i_y = np.around(np.abs(meshgrid_y - mu_y) / 5).astype(np.int64)
i_x = np.around(np.abs(meshgrid_x - mu_x) / 5).astype(np.int64)
i_z = np.abs(meshgrid_z - mu_z).astype(np.int64)
i_y = np.abs(meshgrid_y - mu_y).astype(np.int64)
i_x = np.abs(meshgrid_x - mu_x).astype(np.int64)

# get precomputed values
voxel_integral_z = table_erf_z[i_z, 1]
Expand Down Expand Up @@ -893,58 +892,58 @@ def _objective_function_2d(voxel_size_yx, psf_yx, psf_amplitude=None):
# sigma is known, we fit mu, amplitude and background
if psf_yx is not None and psf_amplitude is None:
def f(grid, mu_y, mu_x, psf_amplitude, psf_background):
values = _gaussian_2d(grid=grid,
mu_y=mu_y,
mu_x=mu_x,
sigma_yx=psf_yx,
voxel_size_yx=voxel_size_yx,
psf_amplitude=psf_amplitude,
psf_background=psf_background)
values = gaussian_2d(grid=grid,
mu_y=mu_y,
mu_x=mu_x,
sigma_yx=psf_yx,
voxel_size_yx=voxel_size_yx,
psf_amplitude=psf_amplitude,
psf_background=psf_background)
return values

# amplitude is known, we fit sigma, mu and background
elif psf_amplitude is not None and psf_yx is None:
def f(grid, mu_y, mu_x, psf_yx, psf_background):
values = _gaussian_2d(grid=grid,
mu_y=mu_y,
mu_x=mu_x,
sigma_yx=psf_yx,
voxel_size_yx=voxel_size_yx,
psf_amplitude=psf_amplitude,
psf_background=psf_background)
values = gaussian_2d(grid=grid,
mu_y=mu_y,
mu_x=mu_x,
sigma_yx=psf_yx,
voxel_size_yx=voxel_size_yx,
psf_amplitude=psf_amplitude,
psf_background=psf_background)
return values

# amplitude and sigma are known, we fit mu and background
elif psf_amplitude is not None and psf_yx is not None:
def f(grid, mu_y, mu_x, psf_background):
values = _gaussian_2d(grid=grid,
mu_y=mu_y,
mu_x=mu_x,
sigma_yx=psf_yx,
voxel_size_yx=voxel_size_yx,
psf_amplitude=psf_amplitude,
psf_background=psf_background)
values = gaussian_2d(grid=grid,
mu_y=mu_y,
mu_x=mu_x,
sigma_yx=psf_yx,
voxel_size_yx=voxel_size_yx,
psf_amplitude=psf_amplitude,
psf_background=psf_background)
return values

# we fit mu, sigma, amplitude and background
else:
def f(grid, mu_y, mu_x, psf_yx, psf_amplitude, psf_background):
values = _gaussian_2d(grid=grid,
mu_y=mu_y,
mu_x=mu_x,
sigma_yx=psf_yx,
voxel_size_yx=voxel_size_yx,
psf_amplitude=psf_amplitude,
psf_background=psf_background)
values = gaussian_2d(grid=grid,
mu_y=mu_y,
mu_x=mu_x,
sigma_yx=psf_yx,
voxel_size_yx=voxel_size_yx,
psf_amplitude=psf_amplitude,
psf_background=psf_background)
return values

return f


def _gaussian_2d(grid, mu_y, mu_x, sigma_yx, voxel_size_yx, psf_amplitude,
psf_background, precomputed=None):
"""Compute the gaussian function over the grid 'xdata' representing a
surface S with shape (S_y, S_x).
def gaussian_2d(grid, mu_y, mu_x, sigma_yx, voxel_size_yx, psf_amplitude,
psf_background, precomputed=None):
"""Compute the gaussian function over the grid representing a surface S
with shape (S_y, S_x).
# TODO add equations
Expand Down Expand Up @@ -987,8 +986,8 @@ def _gaussian_2d(grid, mu_y, mu_x, sigma_yx, voxel_size_yx, psf_amplitude,
table_erf_x = precomputed[1]

# get indices for the tables
i_y = np.around(np.abs(meshgrid_y - mu_y) / 5).astype(np.int64)
i_x = np.around(np.abs(meshgrid_x - mu_x) / 5).astype(np.int64)
i_y = np.abs(meshgrid_y - mu_y).astype(np.int64)
i_x = np.abs(meshgrid_x - mu_x).astype(np.int64)

# get precomputed values
voxel_integral_y = table_erf_y[i_y, 1]
Expand Down Expand Up @@ -1100,7 +1099,7 @@ def _fit_gaussian(f, grid, image_spot, p0, lower_bound=None, upper_bound=None):

def precompute_erf(voxel_size_z=None, voxel_size_yx=100, sigma_z=None,
sigma_yx=200, max_grid=200):
"""Precompute different values for the erf with a resolution of 5 nm.
"""Precompute different values for the erf with a nanometer resolution.
Parameters
----------
Expand Down Expand Up @@ -1131,10 +1130,11 @@ def precompute_erf(voxel_size_z=None, voxel_size_yx=100, sigma_z=None,
sigma_yx=(int, float),
max_grid=int)

# build a grid with a spatial resolution of 5 nm and a size of
# build a grid with a spatial resolution of 1 nm and a size of
# max_grid * resolution nm
yy = np.array([i for i in range(0, int(max_grid * voxel_size_yx), 5)])
xx = np.array([i for i in range(0, int(max_grid * voxel_size_yx), 5)])
max_size_yx = np.ceil(max_grid * voxel_size_yx).astype(np.int64)
yy = np.array([i for i in range(0, max_size_yx)])
xx = np.array([i for i in range(0, max_size_yx)])
mu_y, mu_x = 0, 0

# compute erf values for this grid
Expand All @@ -1155,7 +1155,8 @@ def precompute_erf(voxel_size_z=None, voxel_size_yx=100, sigma_z=None,
return table_erf_y, table_erf_x

else:
zz = np.array([i for i in range(0, int(max_grid * voxel_size_z), 5)])
max_size_z = np.ceil(max_grid * voxel_size_z).astype(np.int64)
zz = np.array([i for i in range(0, max_size_z)])
mu_z = 0
erf_z = _rescaled_erf(low=zz - voxel_size_z / 2,
high=zz + voxel_size_z / 2,
Expand Down

0 comments on commit f544f2e

Please sign in to comment.