From 51cb03a6e3bf46aed2c3e9539af5461757766a2e Mon Sep 17 00:00:00 2001
From: Georgi <georgi.tushev@brain.mpg.de>
Date: Fri, 15 Jul 2022 12:55:42 +0200
Subject: [PATCH] use pyvips buffer
---
.gitignore | 1 +
README.md | 13 ++-
czistitcher.py | 181 +++++++++++++++++++++++++++++++++
pystitcher.py | 255 -----------------------------------------------
requirements.txt | 5 +
5 files changed, 199 insertions(+), 256 deletions(-)
create mode 100644 czistitcher.py
delete mode 100644 pystitcher.py
create mode 100644 requirements.txt
diff --git a/.gitignore b/.gitignore
index e43b0f9..ebf9dd9 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1 +1,2 @@
.DS_Store
+venv/
\ No newline at end of file
diff --git a/README.md b/README.md
index 5be0935..72d666c 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,4 @@
-# PyStitcher
+# CziStitcher
Reads a list of CZI files and stitches based on viewer coordinates
@@ -8,3 +8,14 @@ Reads a list of CZI files and stitches based on viewer coordinates
[aicspylibczi](https://pypi.org/project/aicspylibczi/)
[numpy](https://pypi.org/project/numpy/)
+
+## Installation
+
+```
+$ git clone https://gitlab.mpcdf.mpg.de/mpibr/scic/pystitcher.git
+$ virtualenv -p $(which python3) venv
+$ source venv/bin/activate
+$ pip install --upgrade pip
+$ pip install -r requirements.txt
+$ python czistitcher.py /<path_to_czi_files>/image_name.queue
+```
diff --git a/czistitcher.py b/czistitcher.py
new file mode 100644
index 0000000..50af256
--- /dev/null
+++ b/czistitcher.py
@@ -0,0 +1,181 @@
+import sys
+import os
+import numpy
+import aicspylibczi
+import pyvips
+import re
+
+
+def get_files_per_queue(file_queue: str) -> list:
+ """retrieve all czi files in the stitching queue"""
+ queue = []
+ file_path = os.path.dirname(file_queue)
+ fh = open(file_queue, "r")
+ for file_name in fh:
+ file_name = file_name.strip()
+ if file_name.endswith(".czi"):
+ file_czi = os.path.join(file_path, file_name)
+ queue.append(file_czi)
+ fh.close()
+ return queue
+
+
+def get_scene_grid(czi_queue: list) -> numpy.array:
+ """retrieve bounding box for each czi file in a common grid"""
+ scene_grid = []
+ for file_czi in czi_queue:
+ czi_obj = aicspylibczi.CziFile(file_czi)
+ bbox = czi_obj.get_scene_bounding_box()
+ dims = czi_obj.get_dims_shape()
+ scene_grid.append([bbox.x, bbox.y, bbox.w, bbox.h, dims[0]['Z'][1], dims[0]['C'][1]])
+ scene_grid = numpy.array(scene_grid)
+ scene_grid[:, 0] = scene_grid[:, 0] - scene_grid[:, 0].min()
+ scene_grid[:, 1] = scene_grid[:, 1] - scene_grid[:, 1].min()
+ return scene_grid
+
+
+def get_scene_dims(czi_scene: numpy.array, dim: chr) -> int:
+ if dim == 'X':
+ return czi_scene[:, 0].max() + czi_scene[:, 2].max()
+
+ elif dim == 'Y':
+ return czi_scene[:, 1].max() + czi_scene[:, 3].max()
+
+ elif dim == 'C':
+ return czi_scene[:, -1].max()
+
+ else:
+ return -1
+
+
+# def stitch_scene(czi_queue: list, czi_scene: numpy.array) -> numpy.ndarray:
+# """stitch czi images in a single stack, after max-Z-projection"""
+
+# stack_size_x = get_scene_dims(czi_scene, 'X')
+# stack_size_y = get_scene_dims(czi_scene, 'Y')
+# stack_size_c = get_scene_dims(czi_scene, 'C')
+
+# stack = numpy.zeros((stack_size_x, stack_size_y, stack_size_c), dtype=numpy.uint16)
+# for l in range(0, len(czi_queue)):
+
+# czi_obj = aicspylibczi.CziFile(czi_queue[l])
+# img, shp = czi_obj.read_image()
+# img = numpy.squeeze(img)
+# img = numpy.max(img, axis=1)
+# img = numpy.transpose(img)
+# scene_x = czi_scene[l, 0]
+# scene_y = czi_scene[l, 1]
+# scene_w = czi_scene[l, 2]
+# scene_h = czi_scene[l, 3]
+# print(czi_queue[l], scene_x, scene_y, scene_w, scene_h, img.shape)
+# stack[scene_x:(scene_x + scene_w), scene_y:(scene_y + scene_h), :] = img
+# #break
+
+# return stack
+
+# # print("convert to pyvips")
+# # images = []
+# # for c in range(0, stack_size_c):
+# # images.append(pyvips.Image.new_from_array(stack[:, :, c]))
+# # image = pyvips.Image.arrayjoin(images, across = 1)
+# # image.set_type(pyvips.GValue.gint_type, "page-height", stack_size_y)
+# # return image
+
+
+def stitch_scene(czi_queue: list, czi_scene: numpy.array) -> pyvips.Image:
+ """stitch czi images in a single stack, after max-Z-projection"""
+
+ stack_size_x = get_scene_dims(czi_scene, 'X')
+ stack_size_y = get_scene_dims(czi_scene, 'Y')
+ stack_size_c = get_scene_dims(czi_scene, 'C')
+
+ stack = []
+ for c in range(0, stack_size_c):
+ stack.append(pyvips.Image.black(stack_size_x, stack_size_y).cast(pyvips.enums.BandFormat.USHORT))
+
+ for l in range(0, len(czi_queue)):
+ print(czi_queue[l])
+ czi_obj = aicspylibczi.CziFile(czi_queue[l])
+ img, shp = czi_obj.read_image()
+ img = numpy.squeeze(img) #CZYX
+ img = numpy.max(img, axis=1) #CYX
+ scene_x = czi_scene[l, 0]
+ scene_y = czi_scene[l, 1]
+ for c in range(0, stack_size_c):
+ # pyvips will flip YX to XY with new_from_array()
+ frame = pyvips.Image.new_from_array(img[c, :, :])
+ stack[c] = stack[c].insert(frame, scene_x, scene_y)
+ #break
+
+ image = pyvips.Image.arrayjoin(stack, across = 1)
+ image.set_type(pyvips.GValue.gint_type, "page-height", stack_size_y)
+
+ return image
+
+
+def get_meta_info(file_czi: str) -> str:
+ czi_obj = aicspylibczi.CziFile(file_czi)
+ czi_meta = czi_obj.meta
+ return czi_obj.reader.read_meta()
+
+
+def replace_tag(text, tag, value):
+ pattern = "<" + tag + ">\d+</" + tag + ">"
+ replace = "<" + tag + ">" + str(value) + "</" + tag + ">"
+ return re.sub(pattern, replace, text)
+
+
+if __name__ == "__main__":
+
+ czi_queue = get_files_per_queue(sys.argv[1])
+ czi_scene = get_scene_grid(czi_queue)
+ image = stitch_scene(czi_queue, czi_scene)
+
+
+
+ # metadata = get_meta_info(czi_queue[0])
+ # metadata = replace_tag(metadata, "SizeX", get_scene_dims(czi_scene, 'X'))
+ # metadata = replace_tag(metadata, "SizeY", get_scene_dims(czi_scene, 'Y'))
+ # metadata = replace_tag(metadata, "SizeC", get_scene_dims(czi_scene, 'C'))
+ # metadata = replace_tag(metadata, "SizeZ", 1)
+ # metadata = re.sub("<\?xml version=\"1.0\"\?>\n","", metadata)
+ # metadata = re.sub("\ +<SizeM>\d+</SizeM>\n","", metadata)
+
+ print("write tif ...")
+ outPath = os.path.dirname(czi_queue[0])
+ outName = os.path.basename(czi_queue[0])
+ outName = os.path.splitext(outName)[0] + "_stitched.tif"
+ outFile = os.path.join(outPath, outName)
+
+ image_bands = get_scene_dims(czi_scene, 'C')
+ image_width = get_scene_dims(czi_scene, 'X')
+ image_height = get_scene_dims(czi_scene, 'Y')
+ image.set_type(pyvips.GValue.gstr_type, "image-description",
+ f"""<?xml version="1.0" encoding="UTF-8"?>
+ <OME xmlns="http://www.openmicroscopy.org/Schemas/OME/2016-06"
+ xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
+ xsi:schemaLocation="http://www.openmicroscopy.org/Schemas/OME/2016-06 http://www.openmicroscopy.org/Schemas/OME/2016-06/ome.xsd">
+ <Image ID="Image:0">
+ <!-- Minimum required fields about image dimensions -->
+ <Pixels
+ BigEndian="false"
+ DimensionOrder="XYCZT"
+ ID="Pixels:0"
+ Interleaved="true"
+ PhysicalSizeX="0.10378322801247437"
+ PhysicalSizeXUnit="µm"
+ PhysicalSizeY="0.10378322801247437"
+ PhysicalSizeYUnit="µm"
+ SignificantBits="16"
+ SizeC="{image_bands}"
+ SizeT="1"
+ SizeX="{image_width}"
+ SizeY="{image_height}"
+ SizeZ="1"
+ Type="uint16">
+ </Pixels>
+ </Image>
+ </OME>""")
+
+ #image.set_type(pyvips.GValue.gstr_type, "image-description", metadata)
+ image.write_to_file(outFile, pyramid=True, subifd=True, tile=True, compression="none")
\ No newline at end of file
diff --git a/pystitcher.py b/pystitcher.py
deleted file mode 100644
index 2ec3e12..0000000
--- a/pystitcher.py
+++ /dev/null
@@ -1,255 +0,0 @@
-import sys
-import os
-import numpy
-import aicspylibczi
-import pyvips
-import re
-
-
-def get_files_per_queue(file_queue: str) -> list:
- """retrieve all czi files in the stitching queue"""
- queue = []
- file_path = os.path.dirname(file_queue)
- fh = open(file_queue, "r")
- for file_name in fh:
- file_name = file_name.strip()
- if file_name.endswith(".czi"):
- file_czi = os.path.join(file_path, file_name)
- queue.append(file_czi)
- fh.close()
- return queue
-
-
-def get_scene_grid(czi_queue: list) -> numpy.array:
- """retrieve bounding box for each czi file in a common grid"""
- scene_grid = []
- for file_czi in czi_queue:
- czi_obj = aicspylibczi.CziFile(file_czi)
- bbox = czi_obj.get_scene_bounding_box()
- dims = czi_obj.get_dims_shape()
- scene_grid.append([bbox.x, bbox.y, bbox.w, bbox.h, dims[0]['Z'][1], dims[0]['C'][1]])
- scene_grid = numpy.array(scene_grid)
- scene_grid[:, 0] = scene_grid[:, 0] - scene_grid[:, 0].min()
- scene_grid[:, 1] = scene_grid[:, 1] - scene_grid[:, 1].min()
- return scene_grid
-
-
-# def stitch_scene(czi_queue: list, czi_scene: numpy.array) -> numpy.ndarray:
-# """stitch czi images in a single stack, after max-Z-projection"""
-# stack_size_x = czi_scene[:, 0].max()
-# stack_size_y = czi_scene[:, 1].max()
-# stack_size_c = czi_scene[:, -1].max()
-# stack = numpy.zeros((stack_size_x, stack_size_y, stack_size_c), dtype=numpy.uint16)
-# for l in range(0, len(czi_queue)):
-# print(czi_queue[l])
-# czi_obj = aicspylibczi.CziFile(czi_queue[l])
-# img, shp = czi_obj.read_image()
-# img = numpy.squeeze(img)
-# img = numpy.max(img, axis=1)
-# img = numpy.transpose(img)
-# scene_x = czi_scene[l, 0]
-# scene_y = czi_scene[l, 1]
-# scene_w = czi_scene[l, 2]
-# scene_h = czi_scene[l, 3]
-# stack[scene_x:(scene_x + scene_w), scene_y:(scene_y + scene_h), :] = img
-# break
-# return stack
-
-
-def get_scene_dims(czi_scene: numpy.array, dim: chr) -> int:
- if dim == 'X':
- return czi_scene[:, 0].max()
-
- elif dim == 'Y':
- return czi_scene[:, 1].max()
-
- elif dim == 'C':
- return czi_scene[:, -1].max()
-
- else:
- return -1
-
-
-def stitch_scene(czi_queue: list, czi_scene: numpy.array) -> pyvips.Image:
- """stitch czi images in a single stack, after max-Z-projection"""
-
- stack_size_x = get_scene_dims(czi_scene, 'X')
- stack_size_y = get_scene_dims(czi_scene, 'Y')
- stack_size_c = get_scene_dims(czi_scene, 'C')
-
- stack = []
- for c in range(0, stack_size_c):
- stack.append(pyvips.Image.black(stack_size_x, stack_size_y).cast(pyvips.enums.BandFormat.USHORT))
-
- for l in range(0, len(czi_queue)):
- print(czi_queue[l])
- czi_obj = aicspylibczi.CziFile(czi_queue[l])
- img, shp = czi_obj.read_image()
- img = numpy.squeeze(img)
- img = numpy.max(img, axis=1)
- img = numpy.transpose(img)
- scene_x = czi_scene[l, 0]
- scene_y = czi_scene[l, 1]
- for c in range(0, stack_size_c):
- frame = pyvips.Image.new_from_array(img[:, :, c])
- stack[0] = stack[0].insert(frame, scene_x, scene_y)
- break
-
- image = pyvips.Image.arrayjoin(stack, across = 1)
- image.set_type(pyvips.GValue.gint_type, "page-height", stack_size_y)
-
- return image
-
-
-def get_meta_info(file_czi: str) -> str:
- czi_obj = aicspylibczi.CziFile(file_czi)
- czi_meta = czi_obj.meta
- return czi_obj.reader.read_meta()
-
-
-def replace_tag(text, tag, value):
- pattern = "<" + tag + ">\d+</" + tag + ">"
- replace = "<" + tag + ">" + str(value) + "</" + tag + ">"
- return re.sub(pattern, replace, text)
-
-
-if __name__ == "__main__":
-
- czi_queue = get_files_per_queue(sys.argv[1])
- czi_scene = get_scene_grid(czi_queue)
- image = stitch_scene(czi_queue, czi_scene)
- metadata = get_meta_info(czi_queue[0])
- metadata = replace_tag(metadata, "SizeX", get_scene_dims(czi_scene, 'X'))
- metadata = replace_tag(metadata, "SizeY", get_scene_dims(czi_scene, 'Y'))
- metadata = replace_tag(metadata, "SizeC", get_scene_dims(czi_scene, 'C'))
- metadata = replace_tag(metadata, "SizeZ", 1)
-
-
- metadata = re.sub("<\?xml version=\"1.0\"\?>\n","", metadata)
- metadata = re.sub("\ +<SizeM>\d+</SizeM>\n","", metadata)
-
- image_bands = get_scene_dims(czi_scene, 'C')
- image_width = get_scene_dims(czi_scene, 'X')
- image_height = get_scene_dims(czi_scene, 'Y')
-
-
- image.set_type(pyvips.GValue.gstr_type, "image-description",
- f"""<?xml version="1.0" encoding="UTF-8"?>
- <OME xmlns="http://www.openmicroscopy.org/Schemas/OME/2016-06"
- xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
- xsi:schemaLocation="http://www.openmicroscopy.org/Schemas/OME/2016-06 http://www.openmicroscopy.org/Schemas/OME/2016-06/ome.xsd">
- <Image ID="Image:0">
- <!-- Minimum required fields about image dimensions -->
- <Pixels DimensionOrder="XYCZT"
- ID="Pixels:0"
- SizeC="{image_bands}"
- SizeT="1"
- SizeX="{image_width}"
- SizeY="{image_height}"
- SizeZ="1"
- Type="uint16">
- </Pixels>
- </Image>
- </OME>""")
-
- image.set_type(pyvips.GValue.gstr_type, "image-description", metadata)
- image.write_to_file("/Users/tushevg/Desktop/stitchme/stack.tif", pyramid=True, subifd=True, tile=True, compression="none")
-
-
-# def ExportNDImage(image, fileName, xmlstr):
-# image_width = image.shape[0]
-# image_height = image.shape[1]
-# image_bands = image.shape[2]
-# images = []
-# for c in range(0, image_bands):
-# images.append(pyvips.Image.new_from_array(image[:,:,c]))
-
-# out = pyvips.Image.arrayjoin(images, across=1)
-# out.set_type(pyvips.GValue.gint_type, "page-height", images[0].height)
-# out.set_type(pyvips.GValue.gstr_type, "image-description",
-# f"""<?xml version="1.0" encoding="UTF-8"?>
-# <OME xmlns="http://www.openmicroscopy.org/Schemas/OME/2016-06"
-# xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
-# xsi:schemaLocation="http://www.openmicroscopy.org/Schemas/OME/2016-06 http://www.openmicroscopy.org/Schemas/OME/2016-06/ome.xsd">
-# <Image ID="Image:0">
-# <!-- Minimum required fields about image dimensions -->
-# <Pixels DimensionOrder="XYCZT"
-# ID="Pixels:0"
-# SizeC="{image_bands}"
-# SizeT="1"
-# SizeX="{image_width}"
-# SizeY="{image_height}"
-# SizeZ="1"
-# Type="uint16">
-# </Pixels>
-# </Image>
-# </OME>""")
-# #out.set_type(pyvips.GValue.gstr_type, "image-description", xmlstr)
-# out.write_to_file(fileName, pyramid=True, subifd=True, tile=True, compression="none")
-
-
-# def ExportPyVipsImage(images, fileName, image_width, image_height, image_bands):
-
-# out = pyvips.Image.arrayjoin(images, across=1)
-# out.set_type(pyvips.GValue.gint_type, "page-height", images[0].height)
-# out.set_type(pyvips.GValue.gstr_type, "image-description",
-# f"""<?xml version="1.0" encoding="UTF-8"?>
-# <OME xmlns="http://www.openmicroscopy.org/Schemas/OME/2016-06"
-# xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
-# xsi:schemaLocation="http://www.openmicroscopy.org/Schemas/OME/2016-06 http://www.openmicroscopy.org/Schemas/OME/2016-06/ome.xsd">
-# <Image ID="Image:0">
-# <!-- Minimum required fields about image dimensions -->
-# <Pixels DimensionOrder="XYCZT"
-# ID="Pixels:0"
-# SizeC="{image_bands}"
-# SizeT="1"
-# SizeX="{image_width}"
-# SizeY="{image_height}"
-# SizeZ="1"
-# Type="uint16">
-# </Pixels>
-# </Image>
-# </OME>""")
-# #out.set_type(pyvips.GValue.gstr_type, "image-description", xmlstr)
-# out.write_to_file(fileName, pyramid=True, subifd=True, tile=True, compression="none")
-
-
-# if __name__ == "__main__":
-# fileQueue = "/Users/tushevg/Desktop/projects/InstinctSex/20220502_AFM_15431_I/20220502_AFM_15431_1d_I.queue"
-# queue = GetFilesInQueue(fileQueue)
-# grid = FindGridPoints(queue)
-# ndims = FindNDims(queue)
-
-# # calculate field of view
-# image_width = (grid[:,0] + grid[:,2]).max() + 1
-# image_height = (grid[:,1] + grid[:,3]).max() + 1
-# image_bands = ndims['C']
-
-# # allocate stitched image
-# images = AllocateImages(image_width, image_height, image_bands)
-# #ndimage = numpy.zeros((image_width, image_height, image_bands), dtype="uint16")
-# #print(ndimage.shape)
-
-# # fill
-# k = 0
-# for file in queue:
-# czi = CziFile(file)
-# data = czi.read_image()
-# img = numpy.squeeze(data[0])
-# img = numpy.max(img, axis = 1)
-# #img = numpy.moveaxis(img, [0,1,2], [2,1,0])
-# #xmlstr = ElementTree.tostring(czi.meta, encoding='unicode', method='xml')
-# print(file, img.shape)
-# xpos = grid[k,0]
-# ypos = grid[k,1]
-# #ndimage[xpos:(xpos+img.shape[0]), ypos:(ypos+img.shape[1]), :] = img
-# for c in range(0, ndims['C']):
-# slide = pyvips.Image.new_from_array(img[c,:,:])
-# images[c] = images[c].insert(slide, xpos, ypos)
-# #break
-# k = k + 1
-
-# # export
-# print("export tif")
-# #ExportNDImage(ndimage, "stitched.tif", xmlstr)
-# ExportPyVipsImage(images, "stack.tif", image_width, image_height, image_bands)
diff --git a/requirements.txt b/requirements.txt
new file mode 100644
index 0000000..718666e
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,5 @@
+aicspylibczi==3.0.5
+cffi==1.15.1
+numpy==1.23.1
+pycparser==2.21
+pyvips==2.2.1
--
GitLab