import os
import re
import shutil
import tempfile
from fractions import Fraction
from pathlib import Path
from warnings import warn
import numpy as np
import tifffile
from skimage import img_as_ubyte
from skimage.transform import warp
from owl import mcam_data
from owl.mcam_data import get_bayer_pattern, to_rgb
from owl.stitch.commandline_utils import SilentOutput, directed_print, run_command
from owl.stitch.stitch_utils import (
calculate_output_shape,
lines_intersect,
locate_composite_edges,
locate_transformed_corners_edges,
)
from ..mcam_data._export import _load_metadata
from ..stitch.common import exif_orientation_to_rotation_and_mirror, translation
from ..util import make_timestamp
from .overlap import estimate_overlap
def _tqdm(iterable, *args, **kwargs):
return iterable
def simple_merge(mcam_dataset, global_transforms, downsample=1, blend=False, tqdm=_tqdm):
"""Creates a composite image by of a given mcam_dataset.
Parameters
----------
mcam_dataset : xarray Dataset
Dataset containing the color or grayscale mcam_data.
global_transforms : array of floats
An array of floats of shape (y_cameras, x_cameras, 3, 3) made up of the
3x3 global homography matrices used to stitch merge image.
downsample : int, default 1
Value to downsample the image by. Image is downsampled using Bi-linear
interpolation.
Returns
-------
composite : array of float
Final stitched image as unsigned integers. The shape of this array will
be the minimum size to hold all pixels without clipping as determined
by an individual image's shape and the given global transforms. If three
or four channel color images are given a three or four channel colored
image will be returned, respectively. If grayscale images are given then
a grayscale image will be returned
"""
image_array_tmp = mcam_dataset.images
global_transforms_tmp = global_transforms.copy()
if downsample != 1:
scale_down_matrix = np.array([[1 / downsample, 0, 0],
[0, 1 / downsample, 0],
[0, 0, 1]])
for i in np.ndindex(global_transforms.shape[:2]):
global_transforms_tmp[i] = scale_down_matrix @ global_transforms[i]
output_shape = calculate_output_shape(image_array_tmp.shape[2:4], global_transforms_tmp)
output_shape = tuple(int(np.ceil(o)) for o in output_shape)
if len(image_array_tmp.shape) == 5:
output_shape += (image_array_tmp.shape[4],)
if blend:
h, w = image_array_tmp.shape[2:4]
x_grid, y_grid = np.meshgrid(np.arange(-(w / 2), (w / 2), 1),
np.arange(-(h / 2), (h / 2), 1))
single_mask = np.sqrt(x_grid ** 2 + y_grid ** 2)
single_mask = 1 / single_mask
single_mask -= single_mask.min()
else:
single_mask = np.ones_like(image_array_tmp[0, 0], dtype='float64')
composite = np.zeros(output_shape, dtype='float64')
full_mask = np.zeros_like(composite)
for c, i in tqdm(enumerate(np.ndindex(global_transforms_tmp.shape[:2]))):
warped_image = warp(image_array_tmp[i],
np.linalg.inv(global_transforms_tmp[i]),
output_shape=output_shape
)
single_mask_warped = warp(single_mask,
np.linalg.inv(global_transforms_tmp[i]),
output_shape=output_shape
)
composite += warped_image * single_mask_warped
full_mask += single_mask_warped
full_mask[full_mask == 0] = 1
composite /= full_mask
return img_as_ubyte(composite)
def _get_pto_data(pto_filename, imagename_format):
"""
Pulls all data from the hugin template file and organizes in key value pairs
using the same labels as used in hugin template "i" lines.
http://hugin.sourceforge.net/docs/manual/PTOptimizer.html
Parameters
----------
pto_filename : PathLike
Path to the desired hugin template (.pto) in which to take the position data from.
Returns
-------
images_params : List
List of dicts containing geometric and photometric transformation data.
Length equal to number of images in the template file
"""
with open(pto_filename, encoding='utf-8') as f:
lines = f.readlines()
images_params = []
for line in lines:
if line.startswith('i'):
image_params = {}
image_path = Path(line.split('"')[1])
image_filename = image_path.stem
cam_index = _camera_index_parser(image_filename, imagename_format)
image_params['cam_index'] = cam_index
components = line.split(' ')
for cmp in components:
if not (cmp.startswith(' ') or cmp.startswith('n')):
# need to account for scientific notation...
key = cmp.rstrip('-.1234567890e')
value = cmp[len(key):]
# ...but some keys are supposed to end with e
if value.startswith('e'):
value = value[1:]
key += 'e'
# if the key ends with a `=` then this value is linked to
# the value of a similar key of the image whose index is
# given after the `=`
if key.endswith('='):
key = key[:-1]
value = images_params[int(value)][key]
image_params[key] = value
images_params.append(image_params)
return images_params
def _camera_index_parser(image_filename, imagename_format):
delimiters = re.split('{image_[yx].*?}', Path(imagename_format).stem)
image_filename_cropped = image_filename[len(delimiters[0]):]
if len(delimiters[-1]) > 0:
image_filename_cropped = image_filename_cropped[:-len(delimiters[-1])]
y_index = int(image_filename_cropped.split(delimiters[1])[0])
x_index = int(image_filename_cropped.split(delimiters[1])[1])
return y_index, x_index
def _get_image_shape_from_ptofile(pto_filename, imagename_format):
image_params = _get_pto_data(pto_filename, imagename_format)
width = int(image_params[0]['w'])
height = int(image_params[0]['h'])
return height, width
def hugin_homography(image_shape, rot, scale, shift):
# hugin uses the Mosaic lens model (https://wiki.panotools.org/Stitching_a_photo-mosaic)
# Scale is represented by changing the Z position.
# rotation and scaling are performed about the center of the image
# http://hugin.sourceforge.net/docs/manual/Image_positioning_model.html
image_shape = np.asarray(image_shape)
translate_to_center = translation(-image_shape / 2)
rot_matrix = np.array([[np.cos(rot), -np.sin(rot), 0],
[np.sin(rot), np.cos(rot), 0], # noqa
[0, 0, 1]]) # noqa
scale_matrix = np.array([[scale, 0, 0], # noqa
[0, scale, 0], # noqa
[0, 0, 1]]) # noqa
shift_matrix = translation(shift)
# we do not apply the inverse of `translate_to_center` here so that we can
# accurately reflect the current position of images in hugin's coordinate space.
# Hugin places images without shift into it's coordinate space with the center
# of the image at the origin, while our in out software images without shift
# are place with their corner on the origin.
return shift_matrix @ scale_matrix @ rot_matrix @ translate_to_center
def hugin_homography_array(image_shape, rot_array, scale_array, shift_array):
global_transforms = []
N_cameras = rot_array.shape
for i in np.ndindex(N_cameras):
rot = rot_array[i]
scale = scale_array[i]
shift = shift_array[i]
t = hugin_homography(image_shape, rot, scale, shift)
global_transforms.append(t)
global_transforms = np.reshape(global_transforms, N_cameras + (3, 3))
return global_transforms
def zero_transforms_to_origin(global_transforms, image_shape):
# shift transforms so that no image extends to negative pixels
(top_point,
_,
left_point,
_) = locate_composite_edges(global_transforms, image_shape)
correction_transform = np.linalg.inv(translation((top_point, left_point)))
for i in np.ndindex(global_transforms.shape[:2]):
global_transforms[i] = correction_transform @ global_transforms[i]
return global_transforms
def _check_pto_path(temp_pto_filename, cam_coords, blender, mcam_data_dir):
temp_dir = Path(temp_pto_filename).parent
corrected_temp_pto_filename = temp_dir / 'corrected_template.pto'
number_cameras = len(cam_coords)
if os.path.exists(temp_pto_filename):
i_count = 0
number_cameras_flag = 0
def _create_i_line(line):
components = line.split(' ')
for c in range(len(components)):
if components[c].startswith('n"'):
image_file_path = Path(components[c][2:-2])
image_file_name = image_file_path.stem + image_file_path.suffix
components[c] = (
components[c][:2] +
str(mcam_data_dir / image_file_name) + components[c][-2:]
)
line = (' ').join(components)
return line
def _create_p_line(line, blender):
if blender == 'nona':
tiff_format = 'TIFF'
elif blender == 'enblend':
tiff_format = 'TIFF_m'
elif blender == 'none':
tiff_format = ''
components = line.split(' ')
for c in range(len(components)):
if components[c].startswith('n"'):
components[c] = components[c][:2] + tiff_format
if components[c].startswith('w'):
width = int(components[c][1:])
if components[c].startswith('h'):
height = int(components[c][1:])
line = (' ').join(components)
output_shape = (height, width)
return line, output_shape
# Mark: 2020/12/02
# Putting with statements on multiple lines requires Python 3.9
# https://bugs.python.org/issue12782
with open(temp_pto_filename, 'r', encoding='utf-8') as fin, open(corrected_temp_pto_filename, 'w', encoding='utf-8') as fout: # noqa
# check that contents of template are as expected
for line in fin:
if line.startswith('i'):
i_count += 1
line = _create_i_line(line)
if line.startswith('p'):
line, output_shape = _create_p_line(line, blender)
fout.write(line)
if i_count != number_cameras:
number_cameras_flag = number_cameras
# Code that uses `corrected_temp_pto_path` will want to format the pto path into
# a shell command, so we want to make sure it is formatted as a str for that.
return str(corrected_temp_pto_filename), number_cameras_flag, output_shape
def _hugin_create_pto(pto_filename, mcam_data_dir, cam_coords,
image_shape, imagename_format,
blender, estimated_overlap):
numRows = 1 + cam_coords[-1][0] - cam_coords[0][0]
numColumns = 1 + cam_coords[-1][1] - cam_coords[0][1]
TrX_shift = (1 - estimated_overlap[1]) * np.pi / 180
TrY_shift = (1 - estimated_overlap[0]) * np.pi / 180 * (image_shape[0] / image_shape[1])
with open(pto_filename, 'w+', encoding='utf-8') as f:
if blender == 'nona':
line = 'p f0 n"TIFF c:NONE r:CROP"'
elif blender == 'enblend':
line = 'p f0 n"TIFF_m c:NONE r:CROP"'
elif blender == 'none':
line = ''
f.write(line + '\n')
f.write('# output image lines\n')
# v1 is also used in `pto_file_from_calibrated_images`
s = 'i w{w} h{h} f0 y0 p0 r0 TrX{TrX} TrY{TrY} v1 n"{img_path}"'
for (j, i) in cam_coords:
img_path = str(
mcam_data_dir / imagename_format.format(image_y=j, image_x=i)
)
line = s.format(w=image_shape[1],
h=image_shape[0],
TrX=str(TrX_shift * i - TrX_shift * (numColumns - 1) / 2),
TrY=str(TrY_shift * j - TrY_shift * (numRows - 1) / 2),
img_path=img_path)
f.write(line + '\n')
def _add_optimization_parameters_to_pto_file(pto_filename, cam_coords,
enable_rotation=True, enable_scaling=True):
# optimization parameters
with open(pto_filename, 'a', encoding='utf-8') as f:
components = ['v ']
opt_params = ['TrX', 'TrY', 'Eev']
if enable_rotation:
opt_params.append('r')
if enable_scaling:
opt_params.append('TrZ')
print('opt_params:', opt_params)
[components.append(str(i) + ' ') for i in range(1, len(cam_coords))]
for op in opt_params:
line = op.join(components)
f.write(line + '\n')
def _trim_control_points(pto_filename, cp_threshhold_dist, cam_coords, text_output):
global_transforms, _ = hugin_global_transforms(pto_filename)
editted_cam_coords = []
for cc in cam_coords:
editted_cam_coords.append((cc[0] - cam_coords[0][0], cc[1] - cam_coords[0][1]))
with open(pto_filename, encoding='utf-8') as f:
lines = f.readlines()
trimmed_points = 0
total_points = 0
with open(pto_filename, 'w', encoding='utf-8') as f:
for line in lines:
# control point trimming
if line.startswith('c '):
total_points += 1
comp = line.split(' ')
img1 = int(comp[1][1:])
img2 = int(comp[2][1:])
index1 = editted_cam_coords[img1]
index2 = editted_cam_coords[img2]
gt1 = global_transforms[index1]
gt2 = global_transforms[index2]
# [x_pos, y_pos, 1] to match with global_transforms format
cp1 = [float(comp[3][1:]), float(comp[4][1:]), 1.0]
cp2 = [float(comp[5][1:]), float(comp[6][1:]), 1.0]
transformed_cp1 = gt1 @ cp1
transformed_cp2 = gt2 @ cp2
dist = np.sqrt((transformed_cp1[0] - transformed_cp2[0]) ** 2
+ (transformed_cp1[1] - transformed_cp2[1]) ** 2)
if dist > cp_threshhold_dist:
trimmed_points += 1
continue
f.write(line)
directed_print(f'removed {trimmed_points} of {total_points} control points.\n',
stdout=text_output)
def _my_geocpset(img1, img2, cam_coords, pto_path, global_transforms, image_shape):
editted_cam_coords = []
for cc in cam_coords:
editted_cam_coords.append((cc[0] - cam_coords[0][0], cc[1] - cam_coords[0][1]))
index1 = editted_cam_coords[img1]
index2 = editted_cam_coords[img2]
gt1 = global_transforms[index1]
gt2 = global_transforms[index2]
height, width = image_shape
top_left = np.array([0, 0, 1])[..., None]
top_right = np.array([width, 0, 1])[..., None]
bot_left = np.array([0, height, 1])[..., None]
bot_right = np.array([width, height, 1])[..., None]
corners = [top_left, top_right, bot_left, bot_right]
(img1_top_edge,
img1_bottom_edge,
img1_left_edge,
img1_right_edge) = locate_transformed_corners_edges(gt1, corners)
(img2_top_edge,
img2_bottom_edge,
img2_left_edge,
img2_right_edge) = locate_transformed_corners_edges(gt2, corners)
box1_corners = ((img1_top_edge, img1_left_edge),
(img1_bottom_edge, img1_left_edge),
(img1_top_edge, img1_right_edge),
(img1_bottom_edge, img1_right_edge))
box2_corners = ((img2_top_edge, img2_left_edge),
(img2_bottom_edge, img2_left_edge),
(img2_top_edge, img2_right_edge),
(img2_bottom_edge, img2_right_edge))
if not lines_intersect(box1_corners, box2_corners):
raise RuntimeError(f'image {index1} does not overlap with image {index2}. \
Original alignment is bad. \
Check estimated_overlap or calibrated alignment')
else:
overlap_top_edge = max(img1_top_edge, img2_top_edge)
overlap_bottom_edge = min(img1_bottom_edge, img2_bottom_edge)
overlap_left_edge = max(img1_left_edge, img2_left_edge)
overlap_right_edge = min(img1_right_edge, img2_right_edge)
top_left_control_point = np.array((overlap_left_edge,
overlap_top_edge, 1))[..., None]
bottom_left_control_point = np.array((overlap_left_edge,
overlap_bottom_edge, 1))[..., None]
top_right_control_point = np.array((overlap_right_edge,
overlap_top_edge, 1))[..., None]
bottom_right_control_point = np.array((overlap_right_edge,
overlap_bottom_edge, 1))[..., None]
center_control_point = (top_left_control_point + bottom_left_control_point
+ top_right_control_point + bottom_right_control_point) / 4
control_points = [top_left_control_point,
bottom_left_control_point,
top_right_control_point,
bottom_right_control_point,
center_control_point]
with open(pto_path, 'a', encoding='utf-8') as f:
for control_point in control_points:
img1_control_point = np.linalg.inv(gt1) @ control_point
img2_control_point = np.linalg.inv(gt2) @ control_point
# some out of range points will be generated, but hugin will ignore these points
line = f'c n{img1} N{img2} x{img1_control_point[0][0]} y{img1_control_point[1][0]}'\
+ f' X{img2_control_point[0][0]} Y{img2_control_point[1][0]} t0\n'
f.write(line)
def _resolve_unconnected_images(pto_path, cam_coords, image_shape):
global_transforms, _ = hugin_global_transforms(pto_path)
with open(pto_path, encoding='utf-8') as f:
lines = f.readlines()
images_connectivity = {}
for i in range(len(cam_coords)):
images_connectivity[i] = []
for line in lines:
# record connectivity between images
if line.startswith('c '):
comp = line.split(' ')
img1 = int(comp[1][1:])
img2 = int(comp[2][1:])
if img2 not in images_connectivity[img1]:
images_connectivity[img1].append(img2)
anchor = 0
visited = np.zeros(shape=len(cam_coords), dtype=bool)
visited[anchor] = True
columns = cam_coords[-1][1] - cam_coords[0][1] + 1
# ensure all images are connected to the image above and to the left
for img2 in images_connectivity:
if img2 >= columns:
img1 = img2 - columns
if img2 not in images_connectivity[img1]:
_my_geocpset(img1,
img2,
cam_coords,
pto_path,
global_transforms,
image_shape)
if img2 % columns != 0:
img1 = img2 - 1
if img2 not in images_connectivity[img1]:
_my_geocpset(img1,
img2,
cam_coords,
pto_path,
global_transforms,
image_shape)
def _run_nona(pto_filename, mcam_data_dir, stitch_filename,
cam_coords, imagename_format, use_gpu, text_output):
image_filename = []
for j, i in cam_coords:
image_filename.append('"' + str(
mcam_data_dir / imagename_format.format(image_y=j, image_x=i)
) + '"')
sep = ' '
images_str = sep.join(image_filename)
cmd = 'nona -v --compression=lzw '
if use_gpu:
cmd += '--gpu '
cmd += f'-o "{stitch_filename}" "{pto_filename}" {images_str}'
run_command(cmd, stdout=text_output)
def _run_enblend(*, stitch_filename, load_masks_filename, save_masks_filename,
output_shape, primary_seam_generator, text_output):
image_layers = stitch_filename + '0*'
cmd = f'enblend --compression=lzw --verbose --primary-seam-generator={primary_seam_generator}'
if save_masks_filename is not None:
cmd += f' --save-masks="{save_masks_filename}"'
if load_masks_filename is not None:
cmd += f' --load-masks="{load_masks_filename}"'
# This flag is needed if the output shape is a different aspect ratio from the composite.
# Without it the composite create by enblend will be fitted to the image and ignore the
# desired output shape.
if output_shape is not None:
cmd += f' -f {output_shape[1]}x{output_shape[0]}'
# adding double quotes on image_layers makes enblend
# unable to find images aligned by nona
cmd += f' --output="{stitch_filename}" {image_layers}'
run_command(cmd, stdout=text_output)
def _pixel_shift_to_TrY_TrX(width, shift):
hfov = 1
radius = 180 * width / (hfov * np.pi)
TrX = float(shift[1]) / radius
TrY = float(shift[0]) / radius
return TrY, TrX
def _calculate_current_center(pto_filename, imagename_format):
global_transforms, _ = hugin_global_transforms(pto_filename,
imagename_format=imagename_format,
correct_transforms=False)
image_shape = _get_image_shape_from_ptofile(pto_filename, imagename_format)
(top_point,
bot_point,
left_point,
right_point) = locate_composite_edges(global_transforms, image_shape)
return (top_point + bot_point) / 2, (left_point + right_point) / 2
def _center_pano(pto_filename, imagename_format, desired_center, text_output):
# place the center of the composite at hugin's origin
composite_center = _calculate_current_center(pto_filename,
imagename_format)
_, width = _get_image_shape_from_ptofile(pto_filename,
imagename_format)
shift = (desired_center[0] + composite_center[0],
desired_center[1] + composite_center[1])
# TrY and TrX are the y and x translation variables in the pto
# file's units (radii of the pano sphere)
(TrY_shift,
TrX_shift) = _pixel_shift_to_TrY_TrX(width, shift)
flags = f'--translate={-TrX_shift},{-TrY_shift},0 '
cmd = f'pano_modify -o "{pto_filename}" {flags} "{pto_filename}"'
run_command(cmd, stdout=text_output)
def _resize_and_crop_pano(pto_filename, output_shape, hfov, text_output):
flags = f'--canvas={output_shape[1]}x{output_shape[0]} --fov={hfov} '
cmd = f'pano_modify -o "{pto_filename}" {flags} "{pto_filename}"'
run_command(cmd, stdout=text_output)
def _move_files_from_temp(pto_filename, temp_pto_filename):
out_pto_filename = Path(pto_filename)
if out_pto_filename.exists():
new_out_pto_stem = (out_pto_filename.stem + make_timestamp() + out_pto_filename.suffix)
new_out_pto_filename = out_pto_filename.parent / new_out_pto_stem
msg = f'template file "{out_pto_filename}" already exits. \
Saving template at "{new_out_pto_filename}"'
warn(msg, stacklevel=2)
out_pto_filename = new_out_pto_filename
out_pto_filename = shutil.copyfile(temp_pto_filename, out_pto_filename)
return out_pto_filename
def _copy_stitch(temp_stitch_filename, stitch_filename, metadata):
# make sure the desired path for stitched image doesn't exit,
# if it does add a timestamp
stitch_filename = Path(stitch_filename)
stitch_filename = str(stitch_filename)
shutil.move(temp_stitch_filename, stitch_filename)
_add_metadata_to_image(stitch_filename, metadata)
return stitch_filename
def _add_metadata_to_image(stitch_filename, metadata):
if 'pixel_width' in metadata:
meter_per_pixel = np.asarray(metadata['pixel_width'])
else:
return # don't change metadata if we don't have good info to write
_write_resolution_to_tiff(stitch_filename, meter_per_pixel)
def _write_resolution_to_tiff(stitch_filename, meter_per_pixel):
cm_per_pixel = Fraction(meter_per_pixel * 100)
# values in resolution must not exceed 4294967295 (set by tifffile).
# This correction assumes that `cm_per_pixel` will always be less than 1.
cm_per_pixel = cm_per_pixel.limit_denominator(4294967295)
pixels_per_cm = (1 / cm_per_pixel).as_integer_ratio()
resolution_unit = 3 # centimeter
# r+ means read and write. it will fail with FileNotFoundError if the file does
# not exist.
# b means that the file should be interpreted as a binary file
# as opposed to a text file. This avoids strange encoding errors
# That might appear from treating datafiles as "text"
with tifffile.TiffFile(stitch_filename, mode='r+b') as tiff:
tags = tiff.pages[0].tags
tags['XResolution'].overwrite(pixels_per_cm)
tags['YResolution'].overwrite(pixels_per_cm)
tags['ResolutionUnit'].overwrite(resolution_unit)
def _overlap_to_global_transforms(N_cameras, image_shape, estimated_overlap):
if None in estimated_overlap:
estimated_overlap = estimate_overlap(image_shape, overlap=estimated_overlap)
y_shift = (1 - estimated_overlap[0]) * image_shape[0]
x_shift = (1 - estimated_overlap[1]) * image_shape[1]
global_transforms = np.zeros(shape=N_cameras + (3, 3), dtype=np.float64)
global_transforms[:] = np.eye(3)
for j, i in np.ndindex(N_cameras):
global_transforms[j, i, 0, 2] = i * x_shift
global_transforms[j, i, 1, 2] = j * y_shift
return global_transforms
def _correct_pixel_width(metadata, output_shape_scaling):
# correct pixel width for downsampled/binned images
if 'pixel_width' not in metadata:
return
binning = metadata.x.data[1] - metadata.x.data[0]
metadata.pixel_width.data *= binning
# correct pixel width for changes in composite shape
metadata.pixel_width.data /= output_shape_scaling
return
def _calculate_hfov(full_res_output_shape, output_shape, image_shape):
# when determining fov we need to consider aspect ration changes
if full_res_output_shape[1] / full_res_output_shape[0] < output_shape[1] / output_shape[0]:
vfov = np.arctan(np.tan(np.pi / 360) * full_res_output_shape[0] / image_shape[1])
hfov = np.arctan(np.tan(vfov) * output_shape[1] / output_shape[0]) * 360 / np.pi
else:
hfov = np.arctan(np.tan(np.pi / 360) *
full_res_output_shape[1] / image_shape[1]) * 360 / np.pi
return hfov
def _calc_output_shape_scaling(output_shape, full_res_output_shape):
if full_res_output_shape[1] / full_res_output_shape[0] < output_shape[1] / output_shape[0]:
output_shape_scaling = output_shape[0] / full_res_output_shape[0]
else:
output_shape_scaling = output_shape[1] / full_res_output_shape[1]
return output_shape_scaling
def _bbox_from_indices(bbox_indices, global_transforms):
index_0_vector = np.array([bbox_indices[0][3], bbox_indices[0][2], 1])
index_0_composite_space = global_transforms[bbox_indices[0][:2]] @ index_0_vector
index_1_vector = np.array([bbox_indices[1][3], bbox_indices[1][2], 1])
index_1_composite_space = global_transforms[bbox_indices[1][:2]] @ index_1_vector
y_pos = min(index_0_composite_space[1], index_1_composite_space[1])
x_pos = min(index_0_composite_space[0], index_1_composite_space[0])
height = abs(index_0_composite_space[1] - index_1_composite_space[1])
width = abs(index_0_composite_space[0] - index_1_composite_space[0])
return y_pos, x_pos, height, width
def _apply_exif_orientation(pto_filename, exif_orientation, text_output):
angle_radian, _ = exif_orientation_to_rotation_and_mirror[exif_orientation]
roll = -angle_radian * 180 / np.pi
# Currently this method does not handle the possible mirroring from exif orientation flags
cmd = f'pano_modify --rotate=0,0,{roll} -o "{pto_filename}" "{pto_filename}"'
run_command(cmd, stdout=text_output)
[docs]
def hugin_stitching(mcam_data_path,
*,
save_directory=None,
stitch_filename=None,
pto_filename=None,
template_pto_filename=None,
load_masks_filename=None,
save_masks_filename=None,
blender='enblend',
estimated_overlap=(None, 0.35),
selection_slice=None,
cp_threshhold_dist=100,
ignore_calibration=False,
attempt_custom_alignment=False,
use_gpu=False,
output_shape=None,
bbox_indices=None,
primary_seam_generator='nft',
text_output=None,
**kwargs
):
"""Stitch mcam images using hugin.
stitches mcam images utilizing the hugin toolbox and exports the final
stitched image as a tiff file.
Data can be stitched by giving only the mcam_data directory path as shown:
>>> from pathlib import Path
>>> from owl.stitch.hugin import hugin_stitching
>>> mcam_data_filename = Path('path_to_walkthrough_images')
>>> stitch_filename, pto_filename = hugin_stitching(mcam_data_filename)
Set pass a pto file path to template_pto_path to used premade template instead
of creating a new one
>>> stitch_filename, pto_filename = hugin_stitching(mcam_data_path,
template_pto_filename=pto_filename)
Give a location to save the masks and the blending masks will be saved
during blending
>>> save_masks_filename = mcam_data_filename / 'mask%n.tif'
>>> stitch_filename, pto_filename = hugin_stitching(mcam_data_path,
save_masks_filename=save_masks_filename)
Give a location to save the masks and the blending masks will be saved
during blending
>>> load_masks_filename = save_masks_filename
>>> stitch_filename, pto_filename = hugin_stitching(mcam_data_path,
template_pto_filename=pto_filename,
load_masks_filename=load_masks_filename)
Designate the nona blender (default blender is enblend) - this is quicker,
but seams are more obvious
>>> stitch_filename, pto_filename = hugin_stitching(mcam_data_path,
blender='nona')
Parameters
----------
mcam_data_path: PathLike
The path to either the directory containing all of the exported images to
stitch, or the single file path to the mcam_data.
save_directory: PathLike
If provided, the stitched image and the pto file will be created in the
provided path. By default, if an exported dataset is provided, the
save_directory will be that of the dataset. If a single ``.nc`` file is
provided, a new directory will be created with the same name as the
stem of the ``.nc`` file.
stitch_filename: PathLike, optional
Path to file that will contain the stitched tiff (extension will be appended).
Default is stitched.tif in the ``save_directory``.
pto_filename: PathLike, optional
Path to the pto file containing all the stitching information
(must include pto extension).
Default is template.pto in the ``save_directory``.
template_pto_filename: PathLike
When paseed a pto path the given path will be used to find a prewriten pto file that
will be used to align the images before blending. If it is not found or the
file does not match the N_cameras the code will raise an error. If a template
pto file is used, then a copy of the pto file modified to reflect the used
blender is saved as `template.pto` in the same directory as the stitched
image, or at the file path given to pto_path.
load_masks_filename: String or None
If given a string the masks generated during the blending process are saved
at the location given with the given file extension (tif is suggested).
Include "%n" o add a counter to the name so that all masks have unique names.
Only applies if blender='enblend'. The default is None.
save_masks_filename: String or None
If given a string the blender will attempt to load masks from the filename given.
Use "%n" to reference indexed masks (similar to load_masks_path). Only applies
if blender='enblend'. Defaults to False.
blender: 'enblend', 'nona', OR 'none'
This is the blending engine use for the stitched image. Enblend will use
Dijkstra's shortest path algorithm to minimize the difference between
overlapping pixels along the seam, and nona will selects seams based on the
watershed algorithm. Enblend provides less visible seams, but nona is much
faster. If the blender is set to 'none' only the template will be generated,
without creating a stitched image. The default is 'enblend'.
estimated_overlap: (float, float), optional
Percentage of image size is expected to overlap with neighboring images.
(vertical, horizontal). The default is (0.13, 0.35).
selection_slice: slice
Section of mcam_data array to stitch. If section is used with template
the template must have been made using the same section. If section is
None, full array will be used. When selecting a section of a saved dataset
without camera index (0, 0), index as if the lowest indexed camera is (0, 0).
cp_threshhold_dist: Num
The distance in pixels from the initial positioning that a control point pair
can be. All points pairs a greater distance will be removed.
ignore_calibration: Boolean
If True, ignores the calibrationg stitching transforms when generating
initial alingments and instead using the given estimated overlap.
attempt_custom_alignment: Boolean
If True will attempt to align images based on found control points. If False
will stitch images based on initial position.
use_gpu: bool
If True use the gpu to speed up the performance of `nona`. setting this to True
with `blender='enblend'` decreases the time very slightly, but with `blender='nona'`
the time to stitch decreases about 30%.
output_shape: (int, int)
A tuple holding the desired dimensions of the composite image in (HEIGHT, WIDTH).
bbox_indices: ((image_y0, image_x0, pix_y0, pix_x0), (image_y1, image_x1, pix_y1, pix_x1))
Two diagnolly opposite points to define the bounding box. Points are in data
array space meaning a point is defined by their camera index and pixel index.
primary_seam_generator: str
The primary seam generator used by enblend. Set to either 'graph-cut' ('gc') or
'nearest-feature-transform' ('nft').
Returns
-------
stitch_filename: Path or None
The path to the stitched image. If the blender is set to 'none' this value will be None.
pto_filename: Path
The path to the template file containing the data necessary to create the global transforms.
"""
mcam_data_path = Path(mcam_data_path).resolve().absolute()
if mcam_data_path.suffix == '.nc':
return _hugin_stitching_single_file(mcam_data_filename=mcam_data_path,
save_directory=save_directory,
stitch_filename=stitch_filename,
pto_filename=pto_filename,
template_pto_filename=template_pto_filename,
load_masks_filename=load_masks_filename,
save_masks_filename=save_masks_filename,
blender=blender,
estimated_overlap=estimated_overlap,
selection_slice=selection_slice,
cp_threshhold_dist=cp_threshhold_dist,
ignore_calibration=ignore_calibration,
attempt_custom_alignment=attempt_custom_alignment,
use_gpu=use_gpu,
output_shape=output_shape,
bbox_indices=bbox_indices,
primary_seam_generator=primary_seam_generator,
text_output=text_output,
**kwargs,
)
if mcam_data_path.is_dir():
return _hugin_stitching_exported_data(mcam_data_dir=mcam_data_path,
save_directory=save_directory,
stitch_filename=stitch_filename,
pto_filename=pto_filename,
template_pto_filename=template_pto_filename,
load_masks_filename=load_masks_filename,
save_masks_filename=save_masks_filename,
blender=blender,
estimated_overlap=estimated_overlap,
selection_slice=selection_slice,
cp_threshhold_dist=cp_threshhold_dist,
ignore_calibration=ignore_calibration,
attempt_custom_alignment=attempt_custom_alignment,
use_gpu=use_gpu,
output_shape=output_shape,
bbox_indices=bbox_indices,
primary_seam_generator=primary_seam_generator,
text_output=text_output,
**kwargs,
)
def _hugin_stitching_exported_data(mcam_data_dir,
*,
save_directory=None,
stitch_filename=None,
pto_filename=None,
template_pto_filename=None,
load_masks_filename=None,
save_masks_filename=None,
blender='enblend',
estimated_overlap=(None, 0.35),
selection_slice=None,
cp_threshhold_dist=100,
ignore_calibration=False,
attempt_custom_alignment=False,
use_gpu=False,
output_shape=None,
primary_seam_generator='nft',
bbox_indices=None,
text_output=None,
include_timestamp=True,
enable_rotation=True,
enable_scaling=True,
):
mcam_data_dir = Path(mcam_data_dir)
if save_directory is None:
save_directory = mcam_data_dir
else:
save_directory = Path(save_directory)
filename_stem = mcam_data_dir.stem
if stitch_filename is None:
stitch_filename = save_directory / f'{filename_stem}_stitched.tif'
if pto_filename is None:
pto_filename = save_directory / f'{filename_stem}_template.pto'
Path(stitch_filename).parent.mkdir(parents=True, exist_ok=True)
Path(pto_filename).parent.mkdir(parents=True, exist_ok=True)
if include_timestamp:
stitch_filename = Path(stitch_filename)
pto_filename = Path(pto_filename)
timestamp = make_timestamp()
new_stitch_stem = (stitch_filename.stem + '_' + timestamp + stitch_filename.suffix)
new_stitch_filename = stitch_filename.parent / new_stitch_stem
stitch_filename = new_stitch_filename
new_pto_stem = (pto_filename.stem + '_' + timestamp + pto_filename.suffix)
new_pto_filename = pto_filename.parent / new_pto_stem
pto_filename = new_pto_filename
stitch_filename = str(stitch_filename)
pto_filename = str(pto_filename)
metadata = _load_metadata(mcam_data_dir, metadata_filename=None)
imagename_format = metadata.attrs['imagename_format']
N_cameras = (metadata.sizes['image_y'], metadata.sizes['image_x'])
cam_coords_array = np.empty(shape=N_cameras, dtype=object)
for j, i in np.ndindex(cam_coords_array.shape):
cam_coords_array[j, i] = (metadata['image_y'].data[j], metadata['image_x'].data[i])
cam_coords_array = cam_coords_array[selection_slice]
cam_coords = list(cam_coords_array.flatten())
image_shape = (metadata.sizes['y'], metadata.sizes['x'])
if 'stitch_global_homographies' not in metadata or ignore_calibration:
# translations will be reduced by the binning factor when creating the pto
# file so rather than keeping track of whether these homographies were created
# from estimated overlap or not just scale them up by the binning factor so
# then after binning correction they are correct.
binning = float(metadata.y[1] - metadata.y[0])
unbinned_image_shape = (image_shape[0] * binning,
image_shape[1] * binning)
global_transforms = _overlap_to_global_transforms(N_cameras,
unbinned_image_shape,
estimated_overlap)
dims = ('image_y', 'image_x', 'stitch_j', 'stitch_i')
metadata = metadata.assign({'stitch_global_homographies': (dims, global_transforms)})
metadata = metadata.assign_coords({'stitch_j': np.arange(global_transforms.shape[-2]),
'stitch_i': np.arange(global_transforms.shape[-1])})
metadata['calibration_image_height'] = image_shape[0]
metadata['calibration_image_width'] = image_shape[1]
with tempfile.TemporaryDirectory() as temp_dir:
temp_dir = Path(temp_dir)
temp_pto_filename = str(temp_dir / 'template.pto')
temp_stitch_filename = str(temp_dir / 'stitched.tif')
if template_pto_filename is not None:
if os.path.exists(template_pto_filename):
temp_pto_filename = shutil.copyfile(template_pto_filename, temp_pto_filename)
else:
raise FileNotFoundError(f'No file at "{template_pto_filename}"')
(temp_pto_filename,
number_cameras_flag,
output_shape) = _check_pto_path(temp_pto_filename,
cam_coords,
blender,
mcam_data_dir)
_, full_res_output_shape = hugin_global_transforms(temp_pto_filename,
imagename_format=imagename_format)
output_shape_scaling = output_shape[1] / full_res_output_shape[1]
if number_cameras_flag > 0:
raise OSError(f'given template file "{template_pto_filename}" does not contain '
f'transforms for expected number of images {number_cameras_flag}')
else:
temp_pto_filename = pto_file_from_calibrated_images(
metadata,
pto_filename=temp_pto_filename,
blender=blender,
cam_coords=cam_coords,
imagename_format=imagename_format,
mcam_data_dir=mcam_data_dir,
text_output=text_output)
if attempt_custom_alignment:
# add parameters to optimize
_add_optimization_parameters_to_pto_file(pto_filename=temp_pto_filename,
cam_coords=cam_coords,
enable_rotation=enable_rotation,
enable_scaling=enable_scaling)
# initial control point finder
cmd = f'cpfind --prealigned -o "{temp_pto_filename}" "{temp_pto_filename}"'
run_command(cmd, stdout=text_output)
_trim_control_points(temp_pto_filename, cp_threshhold_dist, cam_coords, text_output)
# add control points to unconnected image sets based on estimated overlap positions
_resolve_unconnected_images(temp_pto_filename, cam_coords, image_shape)
# optimize geometery based on optimization parameters
cmd = f'autooptimiser -n -o "{temp_pto_filename}" "{temp_pto_filename}"'
run_command(cmd, stdout=text_output)
if 'exif_orientation' in metadata:
exif_orientation = int(metadata['exif_orientation'])
_apply_exif_orientation(temp_pto_filename,
exif_orientation,
text_output=text_output)
(global_transforms,
full_res_output_shape) = hugin_global_transforms(temp_pto_filename,
imagename_format=imagename_format)
if bbox_indices is None:
bbox = (0, 0) + full_res_output_shape
else:
bbox = _bbox_from_indices(bbox_indices, global_transforms)
desired_center = ((bbox[0] + bbox[2] / 2) - full_res_output_shape[0] / 2,
(bbox[1] + bbox[3] / 2) - full_res_output_shape[1] / 2)
# center pano
_center_pano(temp_pto_filename,
imagename_format=imagename_format,
desired_center=desired_center,
text_output=text_output)
if output_shape is None:
output_shape = (round(bbox[2]), round(bbox[3]))
output_shape_scaling = _calc_output_shape_scaling(output_shape,
bbox[2:])
hfov = _calculate_hfov(bbox[2:], output_shape, image_shape)
# crop and resize pano
_resize_and_crop_pano(temp_pto_filename,
output_shape=output_shape,
hfov=hfov,
text_output=text_output)
# signal pto file to compress output tiff
cmd = f'pano_modify --ldr-compression=LZW -o "{temp_pto_filename}" ' \
f'"{temp_pto_filename}"'
run_command(cmd, stdout=text_output)
if blender != 'none':
_correct_pixel_width(metadata, output_shape_scaling)
# transform and blend images
_run_nona(temp_pto_filename, mcam_data_dir,
temp_stitch_filename, cam_coords, imagename_format, use_gpu,
text_output=text_output)
if blender == 'enblend':
_run_enblend(stitch_filename=temp_stitch_filename,
load_masks_filename=load_masks_filename,
save_masks_filename=save_masks_filename,
output_shape=output_shape,
primary_seam_generator=primary_seam_generator,
text_output=text_output)
# crop stitched image and move image to desired file
stitch_filename = _copy_stitch(temp_stitch_filename, stitch_filename, metadata)
stitch_filename = Path(stitch_filename)
else:
stitch_filename = None
# move pto file to desired location
out_pto_filename = _move_files_from_temp(pto_filename, temp_pto_filename)
directed_print('Stitching Complete', stdout=text_output)
return stitch_filename, out_pto_filename
def _hugin_stitching_single_file(mcam_data_filename,
*,
save_directory=None,
**kwargs
):
if save_directory is None:
save_directory = mcam_data_filename.parent / mcam_data_filename.stem
else:
save_directory = Path(save_directory)
# Mark: OK I really want to remove chunking in the mcam_data.load
# but this is the one place in production code where we are relying
# on the ability to convert between raw datasets and color ones.
import warnings
with warnings.catch_warnings():
warnings.simplefilter("ignore")
dataset = mcam_data.load(mcam_data_filename)
images = dataset.images
if 'rgb' not in images.coords:
bayer_pattern = get_bayer_pattern(dataset)
images_rgb = to_rgb(images, bayer_pattern=bayer_pattern)
dataset['images'] = images_rgb
with tempfile.TemporaryDirectory() as temp_dir:
temp_images_filename = mcam_data.export(dataset,
temp_dir + f'/{mcam_data_filename.stem}',
include_timestamp=False)
return _hugin_stitching_exported_data(temp_images_filename,
save_directory=save_directory,
**kwargs
)
def pto_file_from_calibrated_images(mcam_dataset, pto_filename,
blender='enblend', cam_coords=None,
imagename_format=None, mcam_data_dir='',
text_output=None):
"""Generate a pto file from the mcam_data's calibrated stitching parameters
Parameters
----------
mcam_dataset : xarrray Dataset
mcam_data containing calibratation data.
pto_path : PathLike
Desired location for the generated pto file. If None the location will be
in mcam_data_dir with the file name `template.pto`.
blender : 'enblend' or 'nona'
The desired blender to be designated by the pto file. 'Nona' will be much
faster, but 'enblend' will produce a more seamles image.
cam_coords : list of tuple of int
The list of all sensors to be used for generating the pto file in the format
[(sensor0_y, sensor0_x), (sensor1_y, sensor1_x), ...] with the top left
sensor listed first and the bottom right sensor listed last. The listed
sensors must come from a connected block. If None is given `cam_coords`
will be generated by the full mcam_data array.
imagename_format : string
The format of the saved images filenames. If none is given will use the
default filename format of `cam{image_y}_{image_x}.bmp`.
mcam_data_dir : PathLike
The directory that holds the images to be stitched. By default, the images are
assumed to be in the same directory as the pto file.
Returns
-------
pto_filename : PathLike
Path to the generated pto file.
"""
if imagename_format is None:
imagename_format = 'cam{image_y}_{image_x}.bmp'
pto_filename = Path(pto_filename)
mcam_data_dir = Path(mcam_data_dir)
pto_filename.parent.mkdir(parents=True, exist_ok=True)
global_transforms = mcam_dataset['stitch_global_homographies'].data
image_shape = (mcam_dataset['calibration_image_height'].data,
mcam_dataset['calibration_image_width'].data)
output_shape = calculate_output_shape(image_shape, global_transforms)
output_shape = (round(output_shape[0]), round(output_shape[1]))
if cam_coords is None:
cam_coords = [i for i in np.ndindex(mcam_dataset.images.shape[:2])]
with open(pto_filename, 'w+', encoding='utf-8') as f:
line = 'p f0 '
if blender == 'nona':
line += 'n"TIFF c:NONE r:CROP"'
elif blender == 'enblend':
line += 'n"TIFF_m c:NONE r:CROP"'
elif blender == 'none':
line = ''
f.write(line + '\n')
f.write('# output image lines\n')
# explanation of these parameters can be found at
# http://hugin.sourceforge.net/docs/nona/nona.txt
s = 'i w{w} h{h} f0 v{hfov} Ra0 Rb0 Rc0 Rd0 Re0 Eev0 Er1 Eb1 r{rot} \
p0 y0 TrX{TrX} TrY{TrY} TrZ{TrZ} Tpy0 Tpp0 j0 a0 b0 c0 d0 e0 g0 \
t0 Va1 Vb0 Vc0 Vd0 Vx0 Vy0 Vm5 n"{img_filename}"'
imagename_format = imagename_format
hfov = 1
for index in cam_coords:
image = mcam_dataset.sel({'image_y': index[0],
'image_x': index[1]})
j = int(image.image_y.data)
i = int(image.image_x.data)
width = len(image.x)
height = len(image.y)
radius = 180 * width / (hfov * np.pi)
(scale,
rotation,
y_translation,
x_translation) = _get_homography_components(image)
rot = np.degrees(rotation)
TrX = float(x_translation) / radius
TrY = float(y_translation) / radius
TrZ = float(scale) - 1
img_filename = mcam_data_dir / imagename_format.format(image_y=j, image_x=i)
line = s.format(w=width,
h=height,
rot=rot,
TrX=TrX,
TrY=TrY,
TrZ=TrZ,
hfov=hfov,
img_filename=img_filename)
f.write(line + '\n')
# reselect pano to properly size image
_center_pano(pto_filename,
imagename_format=imagename_format,
desired_center=(0, 0),
text_output=text_output)
hfov = _calculate_hfov(output_shape, output_shape, image_shape)
_resize_and_crop_pano(pto_filename,
output_shape=output_shape,
hfov=hfov,
text_output=text_output)
return pto_filename
def _hugin_components(global_transform, height, width):
"""
global_transform =
[[scale * cos(rot), -scale * sin(rot), scale * -(Cx * cos(rot) - Cy * sin(rot)) + Cx + Tx],
[scale * sin(rot), scale * cos(rot), scale * -(Cx * sin(rot) + Cy * cos(rot)) + Cy + Ty],
[ 0, 0, 1]]
"""
scale = np.sqrt(global_transform[0, 0] ** 2 + global_transform[1, 0] ** 2)
rotation = np.arcsin(global_transform[1, 0] / scale)
Cx = width / 2
Cy = height / 2
x_translation = global_transform[0, 2] - \
(scale * -(Cx * np.cos(rotation) - Cy * np.sin(rotation)) + Cx)
y_translation = global_transform[1, 2] - \
(scale * -(Cx * np.sin(rotation) + Cy * np.cos(rotation)) + Cy)
return scale, rotation, y_translation, x_translation
def _get_homography_components(image):
global_transform = crop_bin_global_transforms(image)
height = int(image['calibration_image_height'])
width = int(image['calibration_image_width'])
(scale,
rotation,
y_translation,
x_translation) = _hugin_components(global_transform, height, width)
return scale, rotation, y_translation, x_translation
def crop_bin_global_transforms(dataset, resize_binned_data=False):
"""Adjust global transforms in a dataset to account for image cropping and binning.
Computes the shift and scale based on the spatial coordinate differences in the
dataset, and applies them to the global homography transforms to align them with
the full-resolution coordinate system. Optionally preserves the binned coordinate
space.
Parameters
----------
dataset : xarray.Dataset
Dataset containing the `stitch_global_homographies` variable and spatial
coordinates `y` and `x`.
resize_binned_data : bool, optional
If True, keeps the transforms aligned to the binned (downsampled) data.
If False (default), adjusts transforms back to full-resolution space.
Returns
-------
out_global_transforms : numpy.ndarray
Array of shape (..., 3, 3) containing the transformed global homography matrices.
"""
y0 = int(dataset.y[0])
x0 = int(dataset.x[0])
# scaling to account for binning - assume uniform for all images
y1 = int(dataset.y[1])
x1 = int(dataset.x[1])
# While a binning parameter might be helpful to support different scaling in
# x and y, something that the current binning parameter does not support
# as of version 0.14.0
scale_y = y1 - y0
scale_x = x1 - x0
global_transforms = dataset['stitch_global_homographies'].data
return shift_and_scale_global_transforms(global_transforms,
(y0, x0),
(scale_y, scale_x),
resize_binned_data)
def shift_and_scale_global_transforms(global_transforms, shift, scale, resize_binned_data=False):
"""Shift and scale global transform matrices to account for cropping and binning.
Applies a uniform translation shift and scaling operation to a set of global
transformation matrices. Optionally reverses the scaling transformation
to keep the output in the original image resolution space.
Parameters
----------
global_transforms : numpy.ndarray
Array of shape (..., 3, 3) containing transform matrices for each image
in the dataset.
shift : tuple of float
(y, x) tuple specifying the translation shift to apply.
scale : tuple of float
(y_scale, x_scale) tuple specifying the scaling factors due to binning.
resize_binned_data : bool, optional
If True, keeps the output in the binned coordinate space.
If False (default), undoes the binning so output is in original resolution.
Returns
-------
out_global_transforms : numpy.ndarray
Array of shape (..., 3, 3) with the adjusted global transform matrices.
"""
# shifting to account for cropped images - assume uniform for all images
crop_shift = translation(shift)
binning_scale = np.array([[scale[1], 0, 0], # noqa
[0, scale[0], 0], # noqa
[0, 0, 1]]) # noqa
out_global_transforms = (
global_transforms @
crop_shift @
binning_scale
)
if not resize_binned_data:
out_global_transforms = np.linalg.inv(binning_scale) @ out_global_transforms
return out_global_transforms
def hugin_stitching_stack(mcam_dataset_stack_filename,
*,
stack_indices=None,
stitch_filenames=None,
save_directory=None,
include_timestamp=True,
tqdm=None,
reuse_masks=True,
**kwargs):
"""Create multiple composite images from a dataset stack.
Calls hugin_stitching to create a composite at each index of dataset stack.
See Also
--------
hugin_stitching
Parameters
----------
mcam_dataset_stack_filename: PathLike
The path to MCAM dataset stack.
stitch_filenames: list(PathLike), optional
Path to files that will contain the stitched tiffs (extension will be appended).
Default to stitched[index].tif within a directory named after the path holding the
MCAM dataset stack.
stitch_indices: list(int)
The indices of the subsection of layers to be stitched.
save_directory: Path Like
The directory in which to save the stitched images and the template file. Files
will be saved based on the default outlined above.
Returns
-------
stitch_filenames: list(Path)
A list of paths to the stitched images.
pto_filename: Path
The path to the template file containing the stitching data for all composites.
"""
mcam_dataset_stack_filename = Path(mcam_dataset_stack_filename)
dataset_stack = mcam_data.load(mcam_dataset_stack_filename)
stack_size = dataset_stack.images.shape[0]
stack_dim = dataset_stack.images.dims[0]
if stack_indices is None:
stack_indices = np.arange(stack_size)
zero_pad = len(str(max(stack_indices)))
stem = mcam_dataset_stack_filename.stem
if save_directory is None:
save_directory = mcam_dataset_stack_filename.parent / stem
save_directory = Path(save_directory)
if stitch_filenames is None:
stitch_filenames = [save_directory / f'{stem}_stitched_{i:0{zero_pad}d}.tif'
for i in stack_indices]
if 'pto_filename' not in kwargs:
kwargs['pto_filename'] = save_directory / 'template.pto'
if include_timestamp:
timestamp = make_timestamp()
stitch_filenames = [Path(s).parent / f'{Path(s).stem}_{timestamp}{Path(s).suffix}'
for s in stitch_filenames]
pto_filename = Path(kwargs['pto_filename'])
kwargs['pto_filename'] = \
pto_filename.parent / f'{pto_filename.stem}_{timestamp}{pto_filename.suffix}'
kwargs['include_timestamp'] = False
if len(stitch_filenames) != len(stack_indices):
raise ValueError('Given `stitch_filenames` and `stack_indices` must have the same length.')
with tempfile.TemporaryDirectory() as temp_dir:
temp_dir = Path(temp_dir)
dataset = dataset_stack.isel({stack_dim: stack_indices[0]})
if mcam_data.get_chroma(dataset) == 'bayer':
dataset = mcam_data.bayer_dataset_to_rgb(dataset)
temp_images_directory = mcam_data.export(dataset, temp_dir / 'tmp')
if reuse_masks:
save_masks_filename = kwargs.get('save_masks_filename', None)
if save_masks_filename is None:
save_masks_filename = temp_dir / 'mask%n.bmp'
kwargs.update({'save_masks_filename': save_masks_filename})
if tqdm is None:
tqdm = _tqdm
elif 'text_output' not in kwargs:
kwargs['text_output'] = SilentOutput()
stitch_filename = stitch_filenames[0]
stitch_filename, template_pto_filename = hugin_stitching(
temp_images_directory,
stitch_filename=stitch_filename,
**kwargs
)
# in case stitch_filename need to be renamed in `hugin_stitching`
stitch_filenames[0] = stitch_filename
shutil.rmtree(temp_images_directory)
pto_filename = temp_dir / 'tmp1.pto'
kwargs.update({'template_pto_filename': template_pto_filename,
'attempt_custom_alignment': False,
'pto_filename': pto_filename})
if reuse_masks:
kwargs['load_masks_filename'] = save_masks_filename
kwargs['save_masks_filename'] = None
for i, stack_index in tqdm(enumerate(stack_indices[1:]), total=len(stack_indices) - 1):
dataset = dataset_stack.isel({stack_dim: stack_index})
if mcam_data.get_chroma(dataset) == 'bayer':
dataset = mcam_data.bayer_dataset_to_rgb(dataset)
temp_images_directory = mcam_data.export(dataset, temp_dir / 'tmp')
stitch_filename = stitch_filenames[i + 1]
kwargs.update({'stitch_filename': stitch_filename})
(stitch_filename,
pto_filename) = hugin_stitching(temp_images_directory,
**kwargs)
# in case stitch_filename need to be renamed in `hugin_stitching`
stitch_filenames[i + 1] = stitch_filename
pto_filename.unlink()
shutil.rmtree(temp_images_directory)
return stitch_filenames, template_pto_filename