import numpy as np
from skimage.transform import warp
from ..stitch.common import translation
from ..stitch.info import average_overlap
[docs]
def equalize_intensity(images, params, *,
roi_fraction_up_down=(0.3, 0.3),
roi_fraction_left_right=(0.3, 0.3)):
"""Often, images might have slight variations from camera to camera.
This attempts to estimate these intensity variations by looking at the
region around the overlap. Given an image with an overlap of
$(N_y, N_x)$ pixels, this algorithm will analyze a number of pixels equal to
$(N_y * overlapFraction_y, N_x * overlapFraction_x)$ and ensure that the
average intensity in that region is the same.
Parameters
----------
images: mcam_data
must be at least 4D (grayscale), with ``(image_y, image_x)`` as the
leading dimensions.
params: dictionary
Stitching parameters from the fourier stitcher.
roi_fraction_up_down: (float, float)
Tuple describing the fraction of the overlap in each dimension to use
for intensity calibration when comparing the images in the Y dimension.
roi_fraction_left_right: (float, float)
Tuple describing the fraction of the overlap in each dimension to use
for intensity calibration when comparing the images in the X dimension.
Returns
-------
images_out:
mcam_data of the same shape, but float-like, with the intensity
normalized.
intensity_variations:
The relative intensity to which each image was multiplied.
Example
-------
This example will load data in, stitch the data, then use the stitching
parameters to normalize the intensity of the image so as to normalize the
imaging lighting intensity between acquisitions.
This normalization will work best with images of specular targets
illuminated with incoherent light.
>>> from owl import mcam_data
>>> from owl.color import rgb2gray
>>> from owl.analysis import equalize_intensity
>>> from owl.visualize import view_large_image
>>> data = mcam_data.load('your_data_set', output_mode='rgb')
>>> data_gray = rgb2gray(data)
>>> data_float = data.astype('float32') / 255
>>> params = find_stitching_parameters(data_gray)
>>> data_equalized, variations = equalize_intensity(data_float, params)
>>> image_corrected = generate_composite(data_equalized, params)
>>> view_large_image(image_corrected)
"""
left_over_right = pairwise_intensity_ratio(
images, params, axis=1, roi_fraction=roi_fraction_left_right)
# Select only the first column, but make sure we keep the axis
top_over_bottom = pairwise_intensity_ratio(
images[:, :1], params, axis=0, roi_fraction=roi_fraction_up_down)
top_over_bottom = np.cumprod(top_over_bottom, axis=0)
left_over_right = np.cumprod(left_over_right, axis=1)
intensity_variations = left_over_right * top_over_bottom[:, :1]
# Make the largest intensity 1, otherwise we will saturate the image
intensity_variations /= np.max(intensity_variations)
# Broadcast correctly
images_out = images * (
intensity_variations[(Ellipsis,) + (np.newaxis,) * (images.ndim - 2)])
return (images_out, intensity_variations)
def pairwise_intensity_ratio(data_in, params, *,
axis=0, roi_fraction=(0.3, 0.3)):
# Cast ints
dtype = data_in.dtype if np.issubdtype(data_in.dtype, np.floating) else np.float32
transforms = params['transforms']
tributary_transforms = params['tributary_transforms']
tributary_shapes = params['tributary_shapes']
overlap = average_overlap(data_in.shape[2:4], global_transforms=transforms)
# Copy in in the first row of images
intensity_ratio = np.ones(shape=data_in.shape[:2], dtype=dtype)
start_j = 1 if axis == 0 else 0
start_i = 1 if axis == 1 else 0
for j in range(start_j, data_in.shape[0]):
for i in range(start_i, data_in.shape[1]):
# We define the output shape as something that looks
if axis == 1:
output_shape = (int(tributary_shapes[j][i][0] * roi_fraction[0]),
int(overlap[1] * roi_fraction[1] * data_in.shape[2 + 1]))
else: # axis == 0:
output_shape = (int(overlap[0] * roi_fraction[0] * data_in.shape[2 + 0]),
int(tributary_shapes[j][i][1] * roi_fraction[1]))
image_translation = np.asarray(output_shape, dtype=dtype) // 2
the_transform = np.linalg.inv(tributary_transforms[j][i]) @ transforms[j][i]
image_translation[axis] = -output_shape[axis] // 2
warp_translation = the_transform @ translation(image_translation)
this_img_warped = warp(data_in[j, i],
warp_translation,
output_shape=output_shape)
# Other image
j_other = j - 1 if axis == 0 else j
i_other = i - 1 if axis == 1 else i
the_transform = np.linalg.inv(
tributary_transforms[j_other][i_other]) @ transforms[j_other][i_other]
image_translation[axis] += tributary_shapes[j_other][i_other][axis]
warp_translation = the_transform @ translation(image_translation)
# This one should be data out, because it should be the data
# that has already been rescaled
other_img_warped = warp(data_in[j_other, i_other],
warp_translation,
output_shape=output_shape)
this_average = np.mean(this_img_warped)
other_average = np.mean(other_img_warped)
intensity_ratio[j, i] = other_average / this_average
return intensity_ratio