import os
import cv2
import h5py
import imagesize
import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import shapely
from scipy.optimize import least_squares
from scipy.sparse import lil_matrix
from tqdm import tqdm
from .data_formats_io import get_dm4_with_metadata, get_emd_with_metadata
from .image_aligning import close_translation_by_phase_correlation
from .image_processing import img_padding_attenuation, img_to_uint8
[docs]
class image_stack:
"""class to handle stacks of images"""
def __init__(self, mode="storage", no_print=False):
"""
Object to handle large stacks (potentially larger than size of RAM)
for background subtraction, shift correction or stitching
Args:
mode (string):
'storage' only single images are loaded into RAM
when needed
('storage' supports images as .tif, .png, dm4, emd, ...
or as multiple files in .h5 or as datacube in .h5)
The mode 'memory' can be used, for small image series,
that are already being loaded fully into RAM
Defaults to 'storage'
no_print (bool):
set to 'True' for silent execution
Defaults to 'False'
"""
self.mode = mode
self.no_print = no_print
self.directory = None
self.modifiable_file = None
self.overwrite = True
self.units_per_pixel = 1
self.warning_mod_is_none = "please first provide a modifiable dataset" \
" by using 'create_h5cube_duplicate_for_modifying'"
if not no_print:
if mode == "storage":
print("For .h5 provide the path to the h5-file "
"(that serves as h5-directory) via the 'set_directory_path' method")
print("For .tif, .png, ... image files either provide a list of " \
"relative or absolute filepaths via 'set_img_list'")
print("or provide the folder containing only the contributing images " \
"via 'set_directory_path'")
if mode == "memory":
print("please provide the images, as a list of images " \
"via 'set_img_list'")
[docs]
def info(self):
"""print info"""
print("General conditions: ")
print("All images of the series are assumed, to have no rotation and " \
"same pixel size and same dimensions.")
# print("Some functions additionally require same dimensions for all images
# (for example, all image shapes: 512x256).")
[docs]
def set_directory_path(self, path_string):
"""
Set the path to either an h5-file or a folder containing the images
Args:
path_string (string):
looking like, for example: 'folderA/data.h5' or 'folderB'
"""
self.directory = path_string
if path_string[-3:] == ".h5":
self.h5_mode = True
print("please provide the h5 image file paths as a list via 'set_img_list' "
"or use 'imgs_from_datacube'")
else:
rel_img_list = sorted(os.listdir(path_string))
self.img_list = []
for i in range(len(rel_img_list)):
self.img_list.append(path_string + "/" + rel_img_list[i])
self.h5_mode = False
[docs]
def imgs_from_datacube(self, dataset_name):
"""
Images are assumed to be represented by the first index
(for example: 16 x 1024 x 1024)
Args:
dataset_name (string):
name (h5 internal path) of the dataset
could look like for example: 'sampleA/data' or just 'images'
"""
self.mode = "h5_datacube"
self.dataset_name = dataset_name
with h5py.File(self.directory, 'r') as h5:
datashape = h5[dataset_name].shape
self.img_list = np.arange(datashape[0])
[docs]
def change_source_imgs(self, mode="memory", no_print=False):
"""
Change source images of the stack
Args:
mode (string):
'storage' only single images are loaded into RAM
when needed
('storage' supports images as .tif, .png, dm4, emd, ...
or as multiple files in .h5 or as datacube in .h5)
The mode 'memory' can be used, for small image series,
that are already being loaded fully into RAM
Defaults to 'storage'
no_print (bool):
set to 'True' for silent execution
Defaults to 'False'
"""
self.mode = mode
self.directory = None
if not no_print:
if mode == "storage":
print("For .h5 provide the path to the h5-file "
"(that serves as h5-directory) via the 'set_directory_path' method")
print("For .tif, .png, ... image files either provide a list of " \
"relative or absolute filepaths via 'set_img_list'")
print("or provide the folder containing only the contributing images " \
"via 'set_directory_path'")
if mode == "memory":
print("please provide the images, as a list of images " \
"via 'set_img_list'")
[docs]
def set_img_list(self, img_list):
"""
Set the image list
Args:
img_list (list of strings or list of arrays):
corresponding to either paths within the file-system or the h5-hierachy
or in case of mode='memory' a list containing the actual images
"""
if self.directory is None:
self.h5_mode = False
self.img_list = img_list
else:
self.img_list = img_list
[docs]
def get_img(self, index,dset_name=None):
"""
Load the requested image into RAM and return it
Args:
index (int):
index of the image in 'img_list'
Returns
-------
img (array):
requested image
"""
if dset_name:
with h5py.File(self.modifiable_file, 'r') as h5:
img = h5[dset_name][index, :, :]
return img
if self.mode == "memory":
img = self.img_list[index]
elif self.mode == "storage":
if self.h5_mode:
with h5py.File(self.directory, 'r') as h5:
img = h5[self.img_list[index]][()]
else:
if self.img_list[index][-3:]=="dm4":
data,metadata=get_dm4_with_metadata(self.img_list[index])
img=data["data"]
elif self.img_list[index][-3:]=="emd":
img,metadata=get_emd_with_metadata(self.img_list[index])
else:
img = cv2.imread(self.img_list[index], 0)
elif self.mode == "h5_datacube":
with h5py.File(self.directory, 'r') as h5:
img = h5[self.dataset_name][index, :, :]
return img
[docs]
def set_positions(self, positions):
"""
Set the estimated x and y positions of the images for stitching
Args:
postions (array_like or list of tuples):
postions containing x,y
either as array with shape (N,2) or list of tuples
"""
self.positions = positions
[docs]
def set_units_per_pixel(self, units_per_pixel):
"""
Set the units per pixel for example (0.3 micrometer per pixel),
corresponding to positions ('set_positions')
Args:
units_per_pixel (float):
scalar value
"""
self.units_per_pixel = units_per_pixel
[docs]
def sniff_dimensions(self): # ,dimensions=None,no_print=False):
"""
Determine image dimensions from the first image
"""
# if dimensions is not None:
# self.dimensions=dimensions
# else:
self.dimensions = np.zeros([len(self.img_list), 2], dtype=int)
if self.h5_mode:
img = self.get_img(0)
imgshape = img.shape
else:
if self.mode == "storage":
if self.img_list[0][-3:]=="dm4":
data,metadata=get_dm4_with_metadata(self.img_list[0])
imgshape=data["data"].shape
elif self.img_list[0][-3:]=="emd":
img,metadata=get_emd_with_metadata(self.img_list[0])
imgshape= img.shape
else:
imgshape = imagesize.get(self.img_list[0])
else:
imgshape = self.img_list[0].shape
for i in range(len(self.img_list)):
self.dimensions[i] = imgshape
# maybe still usefull, when different dimensions need to be considered
# for i in tqdm(range(len(self.img_list)),disable=no_print):
# if self.h5_mode:
# img=self.get_img(i)
# self.dimensions[i]=img.shape
# else:
# self.dimensions[i]=imagesize.get(self.img_list[i])
[docs]
def overwrite_modifiable_datasets(self,overwrite=True):
"""
set the mode for overwriting (just inside the modifiable dataset)
Args:
overwrite (bool):
set to 'False' to make sure that datasets in the
modifiabel h5-file are conserved
Defaults to 'True'
"""
self.overwrite=overwrite
[docs]
def set_print_mode(self, no_print=False):
"""
Prevent or enable print output
Args:
no_print (bool):
set to 'True' for silent execution
Defaults to 'False'
"""
self.no_print = no_print
[docs]
def set_img(self, index, img, dset_name):
"""
Exchange the image with given index in the modifiable dataset
Args:
index (int):
index corresponding to the image be changed
img (array):
new image that exchanges the old one
dsetname (string):
name of/path to the dataset inside the h5-file
"""
if self.modifiable_file is None:
print(self.warning_mod_is_none)
return 0
with h5py.File(self.modifiable_file, "r+") as f:
f[dset_name][index, :, :] = img
[docs]
def create_h5cube_duplicate_for_modifying(self, filename="temp.h5",
dset_name="data", new_dimensions=None,
force_dtype=False):
"""
Create an h5 file containing the data as 3D-dataset
Args:
filename (string):
name of the h5-file
dset_name (string):
name or internal path to the dataset
new_dimensions (tuple):
if the original dimensions are a multiple of the new_dimensions
cv2.INTER_AREA corresponds to an effective mean value
if None, the original dimensions are used
Defaults to None
force_dtype (bool):
set to 'True' to not change the datatype,
otherwise Float32 is chosen for best compatibility with functions
acting on the data
(datatype of a created dataset can not be changed -> another duplicate
dataset necessary for change)
Defaults to 'False'
"""
self.modifiable_file = filename
if force_dtype:
img = self.get_img(0)
self.dtype = img.dtype
else:
self.dtype = "f4"
if len(np.unique(self.dimensions)) > 2:
# with h5py.File(filename, "w") as f:
# for i in tqdm(range(len(self.img_list)),disable=no_print):
# img=self.get_img(i)
# f["data"][str(i)]=img
print("all images must have equal dimensions")
return 0
else:
if new_dimensions:
factor=self.dimensions[0][0]/new_dimensions[0]
for i in range(len(self.img_list)):
self.dimensions[i] = new_dimensions
self.units_per_pixel *=factor
newshape = [len(self.img_list), self.dimensions[0][0],
self.dimensions[0][1]]
with h5py.File(filename, "w") as f:
f.create_dataset(dset_name, newshape, self.dtype, # dtype="f4")
chunks=(1, self.dimensions[0][0], self.dimensions[0][1]))
for i in tqdm(range(len(self.img_list)), disable=self.no_print):
img = self.get_img(i)
if new_dimensions:
f[dset_name][i, :, :] = cv2.resize(img, new_dimensions,
interpolation=cv2.INTER_AREA)
#img_rebin_by_mean(img,self.dimensions[i])
else:
f[dset_name][i, :, :] = img
[docs]
def subtract_dark_field(self, dark_field, dset_name="data", new_dset_name=None):
"""
substract value or array from whole stack
Args:
dark_field (array_like or float):
scalar value or image for subtraction
dset_name (string):
source of data
new_dset_name (string):
target of data
"""
if self.modifiable_file is None:
print(self.warning_mod_is_none)
return 0
with h5py.File(self.modifiable_file, "r+") as f:
if new_dset_name is None:
new_dset_name = dset_name
else:
newshape = f[dset_name].shape
if self.overwrite:
if new_dset_name in f:
del f[new_dset_name]
f.create_dataset(new_dset_name, newshape, self.dtype,
chunks=(1, self.dimensions[0][0], self.dimensions[0][1]))
for i in range(len(self.img_list)):
f[new_dset_name][i, :, :] = np.maximum(
f[dset_name][i, :, :] - dark_field, 0)
[docs]
def stats(self, histogram=True, histo_levels=100, dset_name=None, mask=None,
original=False):
"""
get statistics from a stack or a tiled image,
obtain minimum, maximum and histogram
Args:
histogram (bool):
set to 'False' for faster calculation
when only minimum and maximum are needed
Defaults to 'True'
histo_levels (int>1):
number of levels of the histogram
(leading to histo_levels + 1 bin borders)
Defaults to 100
Returns:
minimum (float):
global minimum of all images of stack or full tiled image
maximum (float):
global maximum of all images of stack or full tiled image
bins (array_like):
borders of the bins --> histo_levels +1 values
histogram (array_like):
count of each bin
"""
stack_max = -np.inf
stack_min = np.inf
if not original and self.modifiable_file and not dset_name:
dset_name="data"
if self.modifiable_file is None:
original=True
if self.mode=="memory" or original:
if dset_name:
print(dset_name+ " as dset_name not supported in mode 'memory', "
"wihtout modifiable file or setting 'original' True")
return 0
for i in tqdm(range(len(self.img_list)),
disable=self.no_print, desc="min/max"):
img = self.get_img(i)
if histogram:
bin_edges = np.linspace(stack_min, stack_max, histo_levels + 1)
cumulative = np.zeros(histo_levels)
for i in tqdm(range(len(self.img_list)),
disable=self.no_print, desc="histo"):
hist, new_bins = np.histogram(np.ravel(self.get_img(i)),
bins=bin_edges)
cumulative += hist
return stack_min, stack_max, bin_edges, cumulative
else:
return stack_min, stack_max
with h5py.File(self.modifiable_file, "r") as f:
if len(f[dset_name].shape) == 3:
for i in tqdm(range(len(self.img_list)),
disable=self.no_print, desc="min/max"):
img = f[dset_name][i,:,:]
stack_max = max(np.max(img), stack_max)
stack_min = min(np.min(img), stack_min)
if histogram:
bin_edges = np.linspace(stack_min, stack_max, histo_levels + 1)
cumulative = np.zeros(histo_levels)
for i in tqdm(range(len(self.img_list)),
disable=self.no_print, desc="histo"):
hist, new_bins = np.histogram(np.ravel(f[dset_name][i,:,:]),
bins=bin_edges)
cumulative += hist
return stack_min, stack_max, bin_edges, cumulative
else:
return stack_min, stack_max
else:
chunk_rows, chunk_cols = f[dset_name].chunks
n_rows, n_cols = f[dset_name].shape
for row_start in tqdm(range(0, n_rows, chunk_rows),
desc="Rows", disable=self.no_print):
row_end = min(row_start + chunk_rows, n_rows)
# loop over column chunks
for col_start in range(0, n_cols, chunk_cols):
col_end = min(col_start + chunk_cols, n_cols)
# slice the 2D chunk
tile = f[dset_name][row_start:row_end, col_start:col_end]
if mask:
tile = tile[f[mask][row_start:row_end, col_start:col_end]]
if len(tile)!=0:
stack_max = max(np.max(tile), stack_max)
stack_min = min(np.min(tile), stack_min)
if histogram:
bin_edges = np.linspace(stack_min, stack_max, histo_levels + 1)
cumulative = np.zeros(histo_levels)
for row_start in tqdm(range(0, n_rows, chunk_rows),
desc="Rows", disable=self.no_print):
row_end = min(row_start + chunk_rows, n_rows)
# loop over column chunks
for col_start in range(0, n_cols, chunk_cols):
col_end = min(col_start + chunk_cols, n_cols)
# slice the 2D chunk
tile = f[dset_name][row_start:row_end, col_start:col_end]
if mask:
tile = tile[f[mask][row_start:row_end,
col_start:col_end]]
hist, new_bins = np.histogram(tile, bins=bin_edges)
else:
hist, new_bins = np.histogram(np.ravel(tile),
bins=bin_edges)
cumulative += hist
return stack_min, stack_max, bin_edges, cumulative
else:
return stack_min, stack_max
[docs]
def clip(self, vmin=-np.inf, vmax=np.inf, dset_name="data", new_dset_name=None):
"""
clip values of a stack or a tiled image,
Args:
vmin (float):
minimum clip value
Deault vmin negative infinity (no clippling)
vmax (float):
maximum clip value
Deault vmax positive infinity (no clippling)
dset_name (string):
dataset name for source of data
new_dset_name (string):
dataset name for target of data
if new_dset_name is None the data in dset_name will be overwritten
Defaults to None
"""
if self.modifiable_file is None:
print(self.warning_mod_is_none)
return 0
newshape = [len(self.img_list), self.dimensions[0][0], self.dimensions[0][1]]
with h5py.File(self.modifiable_file, "r+") as f:
if len(f[dset_name].shape) == 3:
if new_dset_name is None:
new_dset_name = dset_name
else:
if self.overwrite:
if new_dset_name in f:
del f[new_dset_name]
f.create_dataset(new_dset_name, newshape, dtype="f4",
chunks=(1, self.dimensions[0][0], self.dimensions[0][1]))
# f.create_dataset("data_corrected", newshape,dtype="f4",
# chunks=(1,self.dimensions[0][0],self.dimensions[0][1]))
for i in tqdm(range(len(self.img_list)), disable=self.no_print):
f[new_dset_name][i, :, :] = np.clip(f[dset_name][i, :, :],
vmin, vmax) # clipped
else:
# image=f[dset_name]
chunk_rows, chunk_cols = f[dset_name].chunks
n_rows, n_cols = f[dset_name].shape
if new_dset_name is None:
new_dset_name = dset_name
else:
if self.overwrite:
if new_dset_name in f:
del f[new_dset_name]
f.create_dataset(
new_dset_name,
shape=(n_rows, n_cols),
dtype="float32",
chunks=(chunk_rows, chunk_cols),
# fillvalue=0
)
for row_start in tqdm(range(0, n_rows, chunk_rows),
desc="Rows", disable=self.no_print):
row_end = min(row_start + chunk_rows, n_rows)
# loop over column chunks
for col_start in range(0, n_cols, chunk_cols):
col_end = min(col_start + chunk_cols, n_cols)
# slice the 2D chunk
f[new_dset_name][row_start:row_end, col_start:col_end] = \
np.clip(f[dset_name][row_start:row_end, col_start:col_end],
vmin, vmax)
[docs]
def normalize(self, old_min=None, old_max=None, new_min=0, new_max=1,
dset_name="data", new_dset_name=None):
"""
clip values of a stack or a tiled image,
Args:
old_min (float):
offset subtracted from data for normalization
if None -> old_min is set to global minimum
Deaults to None
old_max (float):
maximum value from data, that will correspond to new_max
after normalization
if None -> old_max is set to global maximum
Deaults to None
new_min (float):
Deaults to 0
new_max (float):
Deaults to 1
dset_name (string):
dataset name for source of data
new_dset_name (string):
dataset name for target of data
if new_dset_name is None the data in dset_name will be overwritten
Defaults to None
"""
if old_min is None or old_max is None:
old_min_auto,old_max_auto=self.stats(histogram=False,dset_name=dset_name)
if not old_min:
old_min = old_min_auto
if not old_max:
old_max = old_max_auto
if self.modifiable_file is None:
print(self.warning_mod_is_none)
return 0
old_delta = old_max - old_min
new_delta = new_max - new_min
factor = new_delta / old_delta
newshape = [len(self.img_list), self.dimensions[0][0], self.dimensions[0][1]]
with h5py.File(self.modifiable_file, "r+") as f:
if len(f[dset_name].shape) == 3:
if new_dset_name is None:
new_dset_name = dset_name
else:
if self.overwrite:
if new_dset_name in f:
del f[new_dset_name]
f.create_dataset(new_dset_name, newshape, dtype="f4",
chunks=(1, self.dimensions[0][0], self.dimensions[0][1]))
# f.create_dataset("data_corrected", newshape,dtype="f4"
# chunks=(1,self.dimensions[0][0],self.dimensions[0][1]))
for i in tqdm(range(len(self.img_list)), disable=self.no_print):
f[new_dset_name][i, :, :] = (f[dset_name][i, :, :] - old_min) * \
factor + new_min
else:
# image=f[dset_name]
chunk_rows, chunk_cols = f[dset_name].chunks
n_rows, n_cols = f[dset_name].shape
if new_dset_name is None:
new_dset_name = dset_name
else:
if self.overwrite:
if new_dset_name in f:
del f[new_dset_name]
f.create_dataset(
new_dset_name,
shape=(n_rows, n_cols),
dtype="float32",
chunks=(chunk_rows, chunk_cols),
# fillvalue=0
)
for row_start in tqdm(range(0, n_rows, chunk_rows),
desc="Rows", disable=self.no_print):
row_end = min(row_start + chunk_rows, n_rows)
# loop over column chunks
for col_start in range(0, n_cols, chunk_cols):
col_end = min(col_start + chunk_cols, n_cols)
# slice the 2D chunk
f[new_dset_name][row_start:row_end, col_start:col_end] = \
(f[dset_name][row_start:row_end, col_start:col_end] \
- old_min) * factor + new_min
[docs]
def flat_field_generation(self, percentile_steps=19, subdiv=4, dset_name="data"):
"""
Calculate a local flat field background of a stack by averaging percentiles
Args:
percentile_steps (int > 1):
number of percentiles used
subdiv (int > 0):
subdiv=1 corresponds to no subdivsion, all data is loaded
subdiv=2 loads only half of the data into RAM
subdiv=3 corresponds to a third ...
adjust according to the size of data and RAM
Defaults to 4
dset_name (string):
dataset name for source of data
Returns:
flat_field (array_like):
resulting flat field image
"""
if self.modifiable_file is None:
print(self.warning_mod_is_none)
return 0
step = 100 / (percentile_steps + 1)
steps = np.linspace(step, 100 - step, percentile_steps)
dims = self.dimensions[0]
flats = np.zeros([percentile_steps, dims[0], dims[1]])
subdiv_steps = int(dims[0] // 4)
if self.mode == "memory":
flats = np.percentile(np.array(self.img_list), steps, axis=0)
else:
with h5py.File(self.modifiable_file, "r") as f:
start_index = 0
for i in range(1, subdiv):
if not self.no_print:
print(str(i) + " out of " + str(subdiv) + " iterations")
flats[:, start_index:i * subdiv_steps, :] = np.percentile(
f[dset_name][:, start_index:i * subdiv_steps, :], steps, axis=0)
if not self.no_print:
print(str(subdiv) + " out of " + str(subdiv) + " iterations")
flats[:, (subdiv - 1) * subdiv_steps:, :] = np.percentile(
f[dset_name][:, (subdiv - 1) * subdiv_steps:, :], steps, axis=0)
# flats=np.percentile(img_series,steps,axis=0)
self.flat_field = np.mean(flats, axis=0)
self.flats = flats
return self.flat_field
[docs]
def dust_from_flat_field(self, flat_field, threshold_block_size=17,
morph_closing_size=5, dark_background=True, dilate=0):
"""
extract positions of dust or dead pixels from flat field image
Args:
flat_field (array_like):
flat field image
threshold_block_size (int, uneven):
block size for binarization using adaptive gaussian thresholding (cv2)
morph_closing_size (int, uneven):
dark_background (bool):
dilate (int >= 0):
Returns:
image_check:
check
"""
th = cv2.adaptiveThreshold(img_to_uint8(flat_field), 255,
cv2.ADAPTIVE_THRESH_GAUSSIAN_C, cv2.THRESH_BINARY,
threshold_block_size, 2)
kernel = np.ones((morph_closing_size, morph_closing_size), np.uint8)
closed = cv2.morphologyEx(th, cv2.MORPH_CLOSE, kernel)
if dark_background:
num, img = cv2.connectedComponents(255 - closed)
else:
num, img = cv2.connectedComponents(closed)
if not self.no_print:
print(str(num) + " dust particles found")
img_d = cv2.dilate(img.astype(np.uint8), np.ones([3, 3], dtype=np.uint8))
last = np.copy(img)
for _ in range(dilate): #_ =i
last = np.copy(img_d)
img_d = cv2.dilate(img_d, np.ones([3, 3], dtype=np.uint8))
rings = img_d - last
self.dust_dict = dict()
for i in range(num - 1):
self.dust_dict[i] = list()
# dust indices
self.dust_dict[i].append(np.where(last == i + 1))
# ring indices
self.dust_dict[i].append(np.where(rings == i + 1))
#self.dust_dict
return last
[docs]
def dust_removal(self, img, mode="median"):
# needs temp data
result = np.copy(img)
for i in self.dust_dict:
if len(self.dust_dict[i][1][0]) < 1:
print("Warning: empty ring")
if mode == "median":
ring_values = np.median(img[self.dust_dict[i][1]])
elif mode == "mean":
ring_values = np.mean(img[self.dust_dict[i][1]])
elif mode == "normal":
ring_values = np.mean(img[self.dust_dict[i][1]])
ring_std = np.std(img[self.dust_dict[i][1]], mean=ring_values)
ring_values = np.random.normal(ring_values, ring_std,
len(img[self.dust_dict[i][1]]))
result[self.dust_dict[i][0]] = ring_values
return result
[docs]
def dust_removal_all(self,dset_name="data", new_dset_name=None):
"""
Apply dust removal on stack
"""
if self.modifiable_file is None:
print(self.warning_mod_is_none)
return 0
if new_dset_name is None:
new_dset_name=dset_name
self.flat_field = self.dust_removal(self.flat_field)
with h5py.File(self.modifiable_file, "r+") as f:
for i in range(len(self.img_list)):
f[new_dset_name][i, :, :] = self.dust_removal(f[dset_name][i, :, :])
[docs]
def flat_field_correction(self, flat_field=None,
dset_name="data", new_dset_name=None):
"""
Apply flat field correction by division
"""
if self.modifiable_file is None:
print(self.warning_mod_is_none)
return 0
if flat_field is None:
flat_field = self.flat_field
if np.min(flat_field) <= 0:
flat_field = flat_field.astype(np.double)
flat_field[flat_field <= 0] = np.inf
flat_field[flat_field == np.inf] = np.min(flat_field)
newshape = [len(self.img_list), self.dimensions[0][0], self.dimensions[0][1]]
with h5py.File(self.modifiable_file, "r+") as f:
if new_dset_name is None:
new_dset_name = dset_name
else:
if self.overwrite:
if new_dset_name in f:
del f[new_dset_name]
f.create_dataset(new_dset_name, newshape, dtype="f4",
chunks=(1, self.dimensions[0][0], self.dimensions[0][1]))
# f.create_dataset("data_corrected", newshape,dtype="f4",
# chunks=(1,self.dimensions[0][0],self.dimensions[0][1]))
for i in tqdm(range(len(self.img_list)), disable=self.no_print):
f[new_dset_name][i, :, :] = f[dset_name][i, :, :] / flat_field
[docs]
def convert_to_uint16(self, real_zero=False,dset_name="data"):
"""
save uint16 dataset in h5
"""
img_maxs = np.zeros(len(self.img_list))
img_mins = np.zeros(len(self.img_list))
with h5py.File(self.modifiable_file, "r") as f:
for i in range(len(self.img_list)):
corr = f[dset_name][i, :, :]
img_maxs[i] = np.max(corr)
img_mins[i] = np.min(corr)
img_min = np.min(img_mins)
img_max = np.max(img_maxs)
div = img_max - img_min
newname = self.modifiable_file[:-3] + "_very_temporary.h5"
newshape = [len(self.img_list),
self.dimensions[0][0],
self.dimensions[0][1]]
with h5py.File(newname, "w") as tf:
tf.create_dataset("data", newshape, np.uint16,
chunks=(1, self.dimensions[0][0], self.dimensions[0][1]))
for i in range(len(self.img_list)):
corr = ((f[dset_name][i, :, :] - img_min) / div) * 65535
tf["data"][i, :, :] = corr.astype(np.uint16)
if real_zero:
newname = self.modifiable_file[:-3] + "_real_zero.h5"
with h5py.File(newname, "w") as tf:
tf.create_dataset("data", newshape, np.uint16,
chunks=(1, self.dimensions[0][0], self.dimensions[0][1]))
for i in range(len(self.img_list)):
corr = (f[dset_name][i, :, :] / img_max) * 65535
tf["data"][i, :, :] = corr.astype(np.uint16)
os.remove(self.modifiable_file)
os.rename(self.modifiable_file[:-3] + "_very_temporary.h5",
self.modifiable_file)
[docs]
def make_polygons(self, units_per_pixel=None, positions=None,
dimensions=None, orientation=2):
"""
create polygon representation
"""
if units_per_pixel is None:
units_per_pixel = self.units_per_pixel
if positions is None:
positions = self.positions
if dimensions is None:
dimensions = self.dimensions
N = len(positions)
anchor_points = np.zeros([N, 2])
polygons = []
im_size = units_per_pixel * dimensions
side0comp = np.ones(N)
side1comp = np.zeros(N)
for i in range(N):
deltax = im_size[i, 1] * np.array([side0comp[i], side1comp[i]])
deltay = im_size[i, 0] * np.array([-side1comp[i], side0comp[i]])
p0 = positions[i] + deltax / 2 - deltay / 2
p1 = p0 - deltax
p2 = p1 + deltay
p3 = p2 + deltax
polygons.append(shapely.geometry.Polygon([p0, p1, p2, p3]))
if orientation == 0:
anchor_points[i] = p0
elif orientation == 1:
anchor_points[i] = p1
elif orientation == 2:
anchor_points[i] = p2
elif orientation == 3:
anchor_points[i] = p3
self.polygons = polygons
self.anchor_points = anchor_points
return polygons, anchor_points
[docs]
def connection_groups(self, polygons=None, units_per_pixel=None,
minimal_number_of_pixels=64, inverse=False):
"""
create connection graph of overlapping polygons and return connected groups
"""
if polygons is None:
polygons = self.polygons
if units_per_pixel is None:
units_per_pixel = self.units_per_pixel
N = len(polygons)
square_units = units_per_pixel * units_per_pixel
G = nx.Graph()
G.add_nodes_from(np.arange(N))
overlaps=[]
for i in range(N - 1):
for j in range(i + 1, N):
inter = shapely.intersection(polygons[i], polygons[j])
area = inter.area / square_units
cond1 = area > minimal_number_of_pixels
if cond1:
overlaps.append(int(np.round(area)))
if inverse:
G.add_edge(i, j, weight=1 / area)
else:
G.add_edge(i, j, weight=overlaps[-1])
con_groups = []
for i in nx.connected_components(G):
con_groups.append(np.array(list(G.subgraph(i).edges)))
self.G = G
if len(con_groups) > 1:
print("not all images can be connected via overlap")
print("sorting out of non-connected images needed")
self.con_group = con_groups[0]
self.max_overlap=np.max(overlaps)
return con_groups
[docs]
def plot_connection_network(self, figsize=None, relative=True):
"""
plot connection graph
"""
if figsize is None:
figsize=[15,15]
plt.figure(figsize=figsize)
if relative:
plt.title("percentage of maximum neighbor-overlap")
else:
plt.title("neighbor-overlap in pixels")
nx.draw(self.G, pos=self.positions, with_labels=True,)
edge_labels = nx.get_edge_attributes(self.G, "weight")
keylist = list(edge_labels.keys())
maxval = 0
for i in keylist:
area = edge_labels[i]
maxval = np.maximum(area, maxval)
for i in keylist:
if relative:
edge_labels[i] = np.round(edge_labels[i] / maxval * 100, 1)
else:
edge_labels[i] = int(np.round(edge_labels[i]))
nx.draw_networkx_edge_labels(self.G, self.positions, edge_labels, label_pos=0.3)
plt.show()
[docs]
def get_outer_polygon_limits(self, polygons):
"""
Obtain bounding box containing all polygons
Args:
polygons (list):
List of polygons
Returns:
points (array_like):
corner points of bounding box as x and y coordinates with shape (4,2)
"""
minx = np.inf
miny = np.inf
maxx = -np.inf
maxy = -np.inf
for polygon in polygons:
# if not isinstance(i,float):
x, y = polygon.exterior.xy
pminx, pmaxx = np.min(x), np.max(x)
pminy, pmaxy = np.min(y), np.max(y)
minx = min(pminx, minx)
miny = min(pminy, miny)
maxx = max(pmaxx, maxx)
maxy = max(pmaxy, maxy)
points = np.zeros([4, 2])
points[0] = minx, miny
points[1] = maxx, miny
points[2] = maxx, maxy
points[3] = minx, maxy
return points
[docs]
def real_to_pixel(self, index, points, orientation=2):
"""
coordinate transformation from real space to pixel space
"""
unit_per_pixel = self.units_per_pixel
anchor_point = self.anchor_points[index]
points = np.vstack((points, anchor_point))
# print("sdfsdf")
# print(points)
points = points[:-1] - points[-1]
points /= unit_per_pixel
points = points[:, ::-1] # *-1
# if orientation==0:
# points[:,1] *= -1
if orientation == 2:
points[:, 0] *= -1
return points # np.round(points).astype(int)
[docs]
def crop_from_points(self, img, points, shape=None):
"""
crop a section of an image conating all points (bounding box around points)
Args:
img (array_like):
points (array_like,list):
shape (tuple,optional):
Returns:
cropped_img (array_like):
"""
vmin, hmin = np.min(points, axis=0)
vmax, hmax = np.max(points, axis=0)
vmin_int = int(np.round(vmin))
vmax_int = int(np.round(vmax))
hmin_int = int(np.round(hmin))
hmax_int = int(np.round(hmax))
# print(shape)
# print("----------")
# print(vmax-vmin)
if shape: # is not None:
if vmax_int - vmin_int > shape[0]: # bigger than intended
lower_residual = vmin_int - vmin # negative makes smaller
upper_residual = vmax - vmax_int # positive makes bigger
if lower_residual < upper_residual:
vmin_int += 1
else:
vmax_int -= 1
elif vmax_int - vmin_int < shape[0]: # smaller than intended
lower_residual = vmin_int - vmin # negative makes smaller
upper_residual = vmax - vmax_int # positive makes bigger
if lower_residual > upper_residual:
vmin_int -= 1
else:
vmax_int += 1
if hmax_int - hmin_int > shape[1]: # bigger than intended
lower_residual = hmin_int - hmin # negative makes smaller
upper_residual = hmax - hmax_int # positive makes bigger
if lower_residual < upper_residual:
hmin_int += 1
else:
hmax_int -= 1
elif hmax_int - hmin_int < shape[1]: # smaller than intended
lower_residual = hmin_int - hmin # negative makes smaller
upper_residual = hmax - hmax_int # positive makes bigger
if lower_residual > upper_residual:
hmin_int -= 1
else:
hmax_int += 1
return img[vmin_int:vmax_int, hmin_int:hmax_int]
[docs]
def pixel_to_real(self, index, points):
"""
coordinate transformation from pixel space to real space
"""
anchor_point = self.anchor_points[index]
unit_per_pixel = self.units_per_pixel
vec1 = np.array([unit_per_pixel, 0])
vec0 = np.array([0, unit_per_pixel])
real_points = np.empty(points.shape)
for i in range(len(points)):
real_points[i] = points[i, 1] * vec1 - points[i, 0] * vec0
real_points[i] += anchor_point
return real_points
[docs]
@staticmethod
def sort_tuples_close_repetition(list_of_tuples):
"""
sort a list of 2-tuples containing integers from 0 to N, such that tuples
with one identical element are in neighbouring positions
Args:
list_of_tuples (list or array_like):
shape: Lx2
Returns:
sorted_tuples (list):
shape: Lx2
"""
N=np.max(list_of_tuples)
tuple_sequence = []
tuple_array = list_of_tuples#np.array(self.G.edges())
booleans = np.ones(len(tuple_array), dtype=bool)
for i in range(N+1):#len(self.img_list)):
for j in range(len(tuple_array)):
if booleans[j]:
if tuple_array[j][0] == i:
tuple_sequence.append(tuple_array[j])
booleans[j] = False
elif tuple_array[j][1] == i:
tuple_sequence.append(tuple_array[j][::-1])
booleans[j] = False
return tuple_sequence
[docs]
def check_pairs(self, max_transl_pix=None, sigma=1, check_data=False,
dset_name=None):
"""
obtain translation vectors from pairwise phase correlation
"""
shifts = []
if check_data:
croplist1 = []
croplist2 = []
self.edge_sequence=self.sort_tuples_close_repetition(np.array(self.G.edges()))
weights=[]
for edge in self.edge_sequence:
weights.append(self.G.get_edge_data(*edge)["weight"]/self.max_overlap)
self.edge_weights=weights
pairs = self.edge_sequence
old_index1 = -1
if dset_name:
with h5py.File(self.modifiable_file, "r") as f:
for pair in tqdm(pairs, disable=self.no_print):
index1 = pair[0]
index2 = pair[1]
poly1 = self.polygons[index1]
poly2 = self.polygons[index2]
overlap = shapely.intersection(poly1, poly2)
overlap_coord = np.array(
overlap.oriented_envelope.exterior.xy).T[:-1]
im1_overlap_coord = self.real_to_pixel(index1,
overlap_coord)
im2_overlap_coord = self.real_to_pixel(index2,
overlap_coord)
im2 = f[dset_name][index2,:,:]#self.get_img(index2)
im2_crop = self.crop_from_points(im2, im2_overlap_coord)
if old_index1 != index1:
im1 = f[dset_name][index1,:,:]#self.get_img(index1)
old_index1 = index1
im1_crop = self.crop_from_points(im1, im1_overlap_coord)
if check_data:
croplist1.append(im1_crop)
croplist2.append(im2_crop)
pix_transl, certainty = close_translation_by_phase_correlation(
im1_crop[:, :], im2_crop[:, :],
sigma=sigma,
max_transl=max_transl_pix)
real_transl = self.units_per_pixel * pix_transl[::-1]
real_transl[1] *= -1
shifts.append(self.positions[index2] - self.positions[index1] +
real_transl)
pcm_distances = dict()
pcm_weights = dict()
edge_tuples = list(map(tuple, pairs))
for i, edge in enumerate(edge_tuples):
pcm_distances[edge] = shifts[i]
pcm_weights[edge]=weights[i]
self.pcm_weights=pcm_weights
self.pcm_distances = pcm_distances
self.edge_tuples = edge_tuples
if check_data:
return pcm_distances, croplist1, croplist2
else:
return pcm_distances
else:
for pair in tqdm(pairs, disable=self.no_print):
index1 = pair[0]
index2 = pair[1]
poly1 = self.polygons[index1]
poly2 = self.polygons[index2]
overlap = shapely.intersection(poly1, poly2)
overlap_coord = np.array(overlap.oriented_envelope.exterior.xy).T[:-1]
im1_overlap_coord = self.real_to_pixel(index1,
overlap_coord)
im2_overlap_coord = self.real_to_pixel(index2,
overlap_coord)
im2 = self.get_img(index2)
im2_crop = self.crop_from_points(im2, im2_overlap_coord)
if old_index1 != index1:
im1 = self.get_img(index1)
old_index1 = index1
im1_crop = self.crop_from_points(im1, im1_overlap_coord)
if check_data:
croplist1.append(im1_crop)
croplist2.append(im2_crop)
pix_transl, certainty = close_translation_by_phase_correlation(
im1_crop[:, :], im2_crop[:, :],
sigma=sigma,
max_transl=max_transl_pix)
real_transl = self.units_per_pixel * pix_transl[::-1]
real_transl[1] *= -1
shifts.append(self.positions[index2] - self.positions[index1] +
real_transl)
pcm_distances = dict()
pcm_weights = dict()
edge_tuples = list(map(tuple, pairs))
for i, edge in enumerate(edge_tuples):
pcm_distances[edge] = shifts[i]
pcm_weights[edge]=weights[i]
self.pcm_weights=pcm_weights
self.pcm_distances = pcm_distances
self.edge_tuples = edge_tuples
if check_data:
return pcm_distances, croplist1, croplist2
else:
return pcm_distances
@staticmethod
def _residuals(posflat, edge_tuples, pcm_distances,pcm_weights):
n = int(len(posflat) // 2)
pts = posflat.reshape(n, 2)
r = []
for i, j in edge_tuples:
r.append(np.linalg.norm(pts[j] - pts[i] - pcm_distances[i, j]) *
pcm_weights[i,j])
return np.array(r)
[docs]
def optimize_positions(self, verbose=True):
"""
global optimization of image positions in real space to minimize deviations
from ideal image translations given by phase correlation
"""
n_res = len(self.edge_tuples)
n_var = len(self.positions) * 2
S = lil_matrix((n_res, n_var))
counter = 0
for i, j in self.edge_tuples:
S[counter, 2 * i] = 1
S[counter, 2 * i + 1] = 1
S[counter, 2 * j] = 1
S[counter, 2 * j + 1] = 1
counter += 1
x0 = np.ravel(self.positions)
if verbose:
verbose_level = 2
else:
verbose_level = 0
result = least_squares(
self._residuals,
x0,
jac_sparsity=S,
args=(self.edge_tuples, self.pcm_distances,self.pcm_weights),
method='trf',
verbose=verbose_level
)
final = result.x.reshape(len(self.positions), 2)
moved_polygons = [
shapely.affinity.translate(poly, xoff=px - poly.centroid.x,
yoff=py - poly.centroid.y)
for poly, (px, py) in zip(self.polygons, final,strict=True) ]
self.moved_polygons=moved_polygons
return moved_polygons
[docs]
def refine_positions(self,moved_polygons,plot_connection=True):
positions=np.array([(poly.centroid.x,poly.centroid.y)
for poly in moved_polygons])
self.positions=positions
self.make_polygons()
self.connection_groups()
if plot_connection:
self.plot_connection_network(relative=False)
[docs]
def map_from_polygons_h5(self, polygons=None, h5file=None, blending="average",
custom_mask=None, boolean_mask=False, border_value=0,
dset_name=None):
"""
create image map by tiling all single images and fusing them together
the data is written into dataset "map" and "division_mask", "boolean_mask"
Args:
polygons (list of poly):
List of polygons, used for fusing the image map
if None, "moved_polygons" is used
h5file (string):
flie name / file path to write the new dataset to
if None, "modifiable_file' is used
Defaults to None
blending (string):
blending between image tiles in overlap area, choose between:
'average' : for using the average value at each pixel of different tiles
'linear' : for a weighted average, with weights linearly decreasing from
1 to 0 going from the center towards the edge of each image
'quadratic' : for a weighted average, with weights decreasing from
1 to 0 like 1-x^2, going from the center towards the edge of each image
'maximum' : use maximum value for pixel in overlapping region
'minimum' : use minimum value for pixel in overlapping region
'hard_cut' : overwrite pixel values in overlapping regions
'custom_single' : give a set of custom weights to weighted average
Defaults to 'average'
custom_mask (array_like):
set of weights in the shape of a single image, only used when
'blending' is set to 'custom_single'
Defaults to None
boolean_mask (bool):
set to 'True' to return a mask with values being 'True' for areas
covered by the stitched image and 'False' for areas around it
Defaults to 'False'
border_value (float):
value used for map borders, when 'blending' is 'minimum' or 'maximum'
Defaults to 0
"""
if h5file is None and self.modifiable_file is None:
print(self.warning_mod_is_none)
return 0
if h5file is None:
h5file = self.modifiable_file
start_index = 0
outer_realspace = self.get_outer_polygon_limits(polygons)
pixelouter = self.real_to_pixel(start_index, outer_realspace)
pixelouter = np.round(pixelouter).astype(int)
offset_x_y = -np.min(pixelouter, axis=0)
image_dims = np.max(pixelouter, axis=0) + offset_x_y
division_needed = blending not in ("hard_cut", "minimum", "maximum")
if blending=="minimum":
fillvalue=np.inf
elif blending=="maximum":
fillvalue=-np.inf
else:
fillvalue=0
with h5py.File(h5file, "a") as f:
if division_needed:
if self.overwrite:
if "division_mask" in f:
del f["division_mask"]
division_mask = f.create_dataset(
"division_mask",
shape=image_dims,
dtype="float32",
chunks=(512, 512),
fillvalue=fillvalue
)
if boolean_mask:
if self.overwrite:
if "boolean_mask" in f:
del f["boolean_mask"]
bmask = f.create_dataset(
"boolean_mask",
shape=image_dims,
dtype="bool",
chunks=(512, 512),
fillvalue=False
)
if self.overwrite:
if "map" in f:
del f["map"]
image = f.create_dataset(
"map",
shape=image_dims,
dtype="float32",
chunks=(512, 512),
fillvalue=0
)
chunk_rows, chunk_cols = image.chunks
imshape = self.dimensions[0]
if blending in ("minimum", "maximum"):
stack = np.empty([imshape[0], imshape[1], 2])
elif blending == "linear":
blending_helper = np.ones(imshape)
blending_helper = img_padding_attenuation(blending_helper,
imshape // 2, mode="linear")
elif blending == "quadratic":
blending_helper = np.ones(imshape)
blending_helper = img_padding_attenuation(blending_helper,
imshape, mode="linear")
for index, poly in enumerate(tqdm(polygons, disable=self.no_print)):
np_realspace = np.array(poly.oriented_envelope.exterior.xy).T[:-1]
np_pixelspace = np.round(self.real_to_pixel(start_index,
np_realspace)).astype(int)
image_pixelspace = np_pixelspace + offset_x_y
if dset_name:
img = f[dset_name][index,:,:]
else:
img = self.get_img(index)
img_start = np.min(image_pixelspace, axis=0)
# img_end=np.max(image_pixelspace,axis=0)
img_end = img_start + img.shape
if blending == "hard_cut":
image[img_start[0]:img_end[0], img_start[1]:img_end[1]] = img
if boolean_mask:
bmask[img_start[0]:img_end[0], img_start[1]:img_end[1]] = True
elif blending == "maximum":
stack[:, :, 0] = \
image[img_start[0]:img_end[0], img_start[1]:img_end[1]]
stack[:, :, 1] = img
image[img_start[0]:img_end[0], img_start[1]:img_end[1]] = np.max(
stack, axis=-1)
if boolean_mask:
bmask[img_start[0]:img_end[0], img_start[1]:img_end[1]] = True
elif blending=="minimum":
stack[:,:,0]=image[img_start[0]:img_end[0],img_start[1]:img_end[1]]
stack[:,:,1]=img
image[img_start[0]:img_end[0],img_start[1]:img_end[1]]=np.min(
stack,axis=-1)
if boolean_mask:
bmask[img_start[0]:img_end[0],img_start[1]:img_end[1]]=True
elif blending == "average":
image[img_start[0]:img_end[0], img_start[1]:img_end[1]] += img
division_mask[img_start[0]:img_end[0], img_start[1]:img_end[1]] += 1
if boolean_mask:
bmask[img_start[0]:img_end[0], img_start[1]:img_end[1]] = True
elif blending == "linear":
image[img_start[0]:img_end[0], img_start[1]:img_end[1]] += \
img * blending_helper
division_mask[img_start[0]:img_end[0], img_start[1]:img_end[1]] += \
blending_helper
if boolean_mask:
bmask[img_start[0]:img_end[0], img_start[1]:img_end[1]] = True
elif blending == "quadratic":
blending_helper = np.ones(img.shape)
imshape = np.array(img.shape)
blending_helper = img_padding_attenuation(blending_helper,
imshape, mode="linear")
image[img_start[0]:img_end[0], img_start[1]:img_end[1]] += \
img * blending_helper
division_mask[img_start[0]:img_end[0], img_start[1]:img_end[1]] += \
blending_helper
if boolean_mask:
bmask[img_start[0]:img_end[0], img_start[1]:img_end[1]] = True
elif blending == "custom_single":
image[img_start[0]:img_end[0], img_start[1]:img_end[1]] += \
img * custom_mask
division_mask[img_start[0]:img_end[0], img_start[1]:img_end[1]] += \
custom_mask
if boolean_mask:
bmask[img_start[0]:img_end[0], img_start[1]:img_end[1]] = True
if blending!="hard_cut":
n_rows, n_cols = image.shape
for row_start in tqdm(range(0, n_rows, chunk_rows),
desc="Rows", disable=self.no_print):
row_end = min(row_start + chunk_rows, n_rows)
# loop over column chunks
for col_start in range(0, n_cols, chunk_cols):
col_end = min(col_start + chunk_cols, n_cols)
if division_needed:
divider = division_mask[row_start:row_end,
col_start:col_end]
divider[divider == 0] = 1
# slice the 2D chunk
image[row_start:row_end, col_start:col_end] = \
image[row_start:row_end,
col_start:col_end] / divider
elif blending == "minimum":
tile=image[row_start:row_end, col_start:col_end]
tile[tile==np.inf]=border_value
image[row_start:row_end, col_start:col_end]=tile
elif blending == "maximum":
tile=image[row_start:row_end, col_start:col_end]
tile[tile==-np.inf]=border_value
image[row_start:row_end, col_start:col_end]=tile
[docs]
def map_from_polygons(self, polygons, blending="average", custom_mask=None,
boolean_mask=False, border_value=0,dset_name=None):
"""
create image map by tiling all single images and fusing them together
Args:
polygons (list of poly):
List of polygons, used for fusing the image map
if None, "moved_polygons" is used
blending (string):
blending between image tiles in overlap area, choose between:
'average' : for using the average value at each pixel of different tiles
'linear' : for a weighted average, with weights linearly decreasing from
1 to 0 going from the center towards the edge of each image
'quadratic' : for a weighted average, with weights decreasing from
1 to 0 like 1-x^2, going from the center towards the edge of each image
'maximum' : use maximum value for pixel in overlapping region
'minimum' : use minimum value for pixel in overlapping region
'hard_cut' : overwrite pixel values in overlapping regions
'custom_single' : give a set of custom weights to weighted average
Defaults to 'average'
custom_mask (array_like):
set of weights in the shape of a single image, only used when
'blending' is set to 'custom_single'
Defaults to None
boolean_mask (bool):
set to 'True' to return a mask with values being 'True' for areas
covered by the stitched image and 'False' for areas around it
Defaults to 'False'
border_value (float):
value used for map borders, when 'blending' is 'minimum' or 'maximum'
Defaults to 0
"""
start_index = 0
outer_realspace = self.get_outer_polygon_limits(polygons)
pixelouter = self.real_to_pixel(start_index, outer_realspace)
pixelouter = np.round(pixelouter).astype(int)
offset_x_y = -np.min(pixelouter, axis=0)
image_dims = np.max(pixelouter, axis=0) + offset_x_y
division_needed = blending not in ("hard_cut", "minimum", "maximum")
if division_needed:
division_mask = np.zeros(image_dims) # ,dtype=np.uint8)
image = np.zeros(image_dims)
if boolean_mask:
bmask = np.zeros(image_dims, dtype=bool)
if blending == "minimum":
image += np.inf
elif blending == "maximum":
image += -np.inf
for index, poly in enumerate(tqdm(polygons, disable=self.no_print)):
np_realspace = np.array(poly.oriented_envelope.exterior.xy).T[:-1]
np_pixelspace = np.round(self.real_to_pixel(start_index,
np_realspace)).astype(int)
image_pixelspace = np_pixelspace + offset_x_y
if dset_name:
with h5py.File(self.modifiable_file, "r") as f:
img=f[dset_name][index,:,:]
else:
img = self.get_img(index)
img_start = np.min(image_pixelspace, axis=0)
# img_end=np.max(image_pixelspace,axis=0)
img_end = img_start + img.shape
if blending == "hard_cut":
image[img_start[0]:img_end[0], img_start[1]:img_end[1]] = img
if boolean_mask:
bmask[img_start[0]:img_end[0], img_start[1]:img_end[1]] = True
elif blending == "maximum":
stack = np.empty([img.shape[0], img.shape[1], 2])
stack[:, :, 0] = image[img_start[0]:img_end[0], img_start[1]:img_end[1]]
stack[:, :, 1] = img
image[img_start[0]:img_end[0], img_start[1]:img_end[1]] = \
np.max(stack, axis=-1)
if boolean_mask:
bmask[img_start[0]:img_end[0], img_start[1]:img_end[1]] = True
elif blending == "minimum":
stack = np.empty([img.shape[0], img.shape[1], 2])
stack[:, :, 0] = image[img_start[0]:img_end[0], img_start[1]:img_end[1]]
stack[:, :, 1] = img
image[img_start[0]:img_end[0], img_start[1]:img_end[1]] = \
np.min(stack, axis=-1)
if boolean_mask:
bmask[img_start[0]:img_end[0], img_start[1]:img_end[1]] = True
elif blending == "linear":
blending_helper = np.ones(img.shape)
imshape = np.array(img.shape)
blending_helper = img_padding_attenuation(blending_helper,
imshape // 2, mode="linear")
image[img_start[0]:img_end[0], img_start[1]:img_end[1]] += \
img * blending_helper
division_mask[img_start[0]:img_end[0], img_start[1]:img_end[1]] += \
blending_helper
if boolean_mask:
bmask[img_start[0]:img_end[0], img_start[1]:img_end[1]] = True
elif blending == "quadratic":
blending_helper = np.ones(img.shape)
imshape = np.array(img.shape)
blending_helper = img_padding_attenuation(blending_helper,
imshape, mode="linear")
image[img_start[0]:img_end[0], img_start[1]:img_end[1]] += \
img * blending_helper
division_mask[img_start[0]:img_end[0], img_start[1]:img_end[1]] += \
blending_helper
if boolean_mask:
bmask[img_start[0]:img_end[0], img_start[1]:img_end[1]] = True
elif blending == "average":
image[img_start[0]:img_end[0], img_start[1]:img_end[1]] += img
division_mask[img_start[0]:img_end[0], img_start[1]:img_end[1]] += 1
if boolean_mask:
bmask[img_start[0]:img_end[0], img_start[1]:img_end[1]] = True
elif blending == "custom_single":
image[img_start[0]:img_end[0], img_start[1]:img_end[1]] += \
img * custom_mask
division_mask[img_start[0]:img_end[0], img_start[1]:img_end[1]] += \
custom_mask
if boolean_mask:
bmask[img_start[0]:img_end[0], img_start[1]:img_end[1]] = True
elif blending == "custom_multi":
image[img_start[0]:img_end[0], img_start[1]:img_end[1]] += \
img * custom_mask[index]
division_mask[img_start[0]:img_end[0], img_start[1]:img_end[1]] += \
custom_mask[index]
if boolean_mask:
bmask[img_start[0]:img_end[0], img_start[1]:img_end[1]] = True
if blending == "minimum":
image[image == np.inf] = border_value
elif blending == "maximum":
image[image == -np.inf] = border_value
if division_needed:
division_mask[division_mask == 0] = 1
image[:, :] /= division_mask
if boolean_mask:
return image, bmask
else:
return image