diff --git a/extract_from_platybrowser.py b/extract_from_platybrowser.py new file mode 100644 index 0000000000000000000000000000000000000000..0e3e7686a02802b3d9b42342b18ccbaa1c4a62d1 --- /dev/null +++ b/extract_from_platybrowser.py @@ -0,0 +1,28 @@ +from scripts.export.extract_subvolume import make_cutout + + +if __name__ == '__main__': + levels = get_res_level() + res_string = " ".join("%i: %s" % (lev, res) for lev, res in enumerate(levels)) + scale_help = ("The resolution level for the cutout." + "The levels range from 0 to 6 with resolutions in micrometer: " + "%s" % res_string) + + parser = argparse.ArgumentParser("Make a cutout from the platynereis EM data-set and save it as tif stack.") + parser.add_argument("scale_level", type=int, help=scale_help) + parser.add_argument("lower_corner", type=str, help="Lower corner of the bounding box to cut out") + parser.add_argument("upper_corner", type=str, help="Upper corner of the bounding box to cut out") + parser.add_argument("save_file", type=str, help="Where to save the cutout.") + args = parser.parse_args() + + print("Converting coordinates from physical coordinates %s:%s" % (args.lower_corner, + args.upper_corner)) + bb_start = parse_coordinate(args.lower_corner) + bb_stop = parse_coordinate(args.upper_corner) + print("to pixel coordinates %s:%s" % (str(bb_start), + str(bb_stop))) + + print("Extracting raw data") + raw = make_cutout(args.scale_level, bb_start, bb_stop) + print("Saving data to tif-stack %s" % args.save_file) + save_tif(raw, args.save_file) diff --git a/scripts/export/extract_subvolume.py b/scripts/export/extract_subvolume.py new file mode 100644 index 0000000000000000000000000000000000000000..b9e26e5dba36ae634eec3f5ce2736662471dc1ec --- /dev/null +++ b/scripts/export/extract_subvolume.py @@ -0,0 +1,87 @@ +import argparse +import os +import h5py +import imageio + +from ..files.xml_utils import get_h5_path_from_xml + + +def get_res_level(level=None): + res0 = [.01, .01, .025] + res1 = [.02, .02, .025] + resolutions = [res0] + [[re * 2 ** i for re in res1] for i in range(8)] + if level is None: + return resolutions + else: + assert level <= 6,\ + "Scale level %i is not supported, only supporting up to level 8" % level + return resolutions[level] + + +def save_tif_stack(raw, save_file): + try: + imageio.volwrite(save_file, raw) + return True + except RuntimeError: + return False + + +def save_tif_slices(raw, save_file): + print("Could not save tif stack, saving slices to folder %s instead" + % save_file) + save_folder = os.path.splitext(save_file)[0] + os.makedirs(save_folder, exist_ok=True) + for z in range(raw.shape[0]): + out_file = os.path.join(save_folder, "z%05i.tif" % z) + io.imsave(out_file, raw[z]) + + +def save_tif(raw, save_file): + if imageio is None: + save_tif_slices(raw, save_file) + else: + if not save_tif_stack(raw, save_file): + save_tif_slices(raw, save_file) + + +def name_to_path(name): + name_dict = {'raw': 'images/', + 'cells': 'segmentations/', + 'nuclei': 'segmentations/', + 'cilia': 'segmentations/', + 'chromatin': 'segmentations/'} + assert name in name_dict, "Name must be one of %s, not %s" % (str(name_dict.keys()), + name) + return name_dict[name] + + +def make_cutout(tag, name, scale, bb_start, bb_stop): + assert all(sta < sto for sta, sto in zip(bb_start, bb_stop)) + + path = os.path.join('data', tag, name_to_path(name)) + data_resolution = read_resolution(path) + path = get_h5_path_from_xml(path, return_absolute_path=True) + resolution = get_res_level(scale) + + data_scale = match_scale(resolution, data_resolution) + + bb_start = [int(sta / re) for sta, re in zip(bb_start, resolution)][::-1] + bb_stop = [int(sto / re) for sto, re in zip(bb_stop, resolution)][::-1] + bb = tuple(slice(sta, sto) for sta, sto in zip(bb_start, bb_stop)) + + key = 't00000/s00/%i/cells' % data_scale + with h5py.File(path, 'r') as f: + ds = f[key] + data = ds[bb] + return data + + +def parse_coordinate(coord): + coord = coord.rstrip('\n') + pos_start = coord.find('(') + 1 + pos_stop = coord.find(')') + coord = coord[pos_start:pos_stop] + coord = coord.split(',') + coord = [float(co) for co in coord] + assert len(coord) == 3, "Coordinate conversion failed" + return coord