diff --git a/DriftCorrector.py b/DriftCorrector.py
new file mode 100644
index 0000000000000000000000000000000000000000..d5decd7ba38a4e87c720678fc849bb1741eef4ee
--- /dev/null
+++ b/DriftCorrector.py
@@ -0,0 +1,22 @@
+import src.UI as UI
+
+def create_socket_or_die():
+ '''Create a socket to prevent multiple instances of this program'''
+ import socket
+
+ HOST = 'localhost'
+ PORT = 44528
+ s = socket.socket()
+
+ try:
+ s.bind((HOST, PORT))
+ except OSError:
+ # Program is already running
+ exit()
+
+ return s
+
+if __name__ == '__main__':
+ s = create_socket_or_die()
+ app = UI.VesicleAnnotationInterface()
+ s.close()
diff --git a/DriftCorrector.spec b/DriftCorrector.spec
new file mode 100644
index 0000000000000000000000000000000000000000..a2e9f4f20d197e92988b6f4ca131924934593876
--- /dev/null
+++ b/DriftCorrector.spec
@@ -0,0 +1,36 @@
+# -*- mode: python -*-
+
+block_cipher = None
+
+
+a = Analysis(['DriftCorrector.py'],
+ pathex=['/home/dith/arngorf/code/drift_corrector'],
+ binaries=[],
+ datas=[],
+ hiddenimports=[],
+ hookspath=['.'],
+ runtime_hooks=[],
+ excludes=[],
+ win_no_prefer_redirects=False,
+ win_private_assemblies=False,
+ cipher=block_cipher,
+ noarchive=False)
+pyz = PYZ(a.pure, a.zipped_data,
+ cipher=block_cipher)
+exe = EXE(pyz,
+ a.scripts,
+ [],
+ exclude_binaries=True,
+ name='DriftCorrector',
+ debug=False,
+ bootloader_ignore_signals=False,
+ strip=False,
+ upx=True,
+ console=True )
+coll = COLLECT(exe,
+ a.binaries,
+ a.zipfiles,
+ a.datas,
+ strip=False,
+ upx=True,
+ name='DriftCorrector')
diff --git a/README.md b/README.md
index 13eb0620f216610a1a67925f636e66e5e365d7ff..778a571f1f0fdbfa6795c8d6baa1b0c6e0688596 100644
--- a/README.md
+++ b/README.md
@@ -14,3 +14,27 @@ This project contains a small program with user interface implementing the drift
publisher={Nature Publishing Group}
}
```
+
+# How to use
+
+1) $python DriftCorrector.py
+
+2) Click "Open dataset" to get startet
+
+3) navigate to the directory with your volume data as either a multipage tiff file or in the a number of consecutively named image files. I.e. my_image_0.png, my_image_1.png, my_image_2.png, ...
+
+4) Confirm that the tool is in "Add Points" mode in the left pane
+
+5) Annotate vesicles one at a time by using the following commands
+
+ Right-click to center at a vesicle
+ Zoom in/out using the buttons in the left pane
+ Left-click a vesicle boundary to annotate a point
+ Use the mouse scroll button or the "+" and "-" key to cycle through the image sections
+ Change viewing direction by clicking "a"
+ Undo point placement by clicking "ctrl+z"
+
+6) Click "Edit Points" to edit or delete previously placed points, do:
+
+ Right-click a previous annotation to make it active
+ Left-click a point to delete it
\ No newline at end of file
diff --git a/docs/images/add_points_01.png b/docs/images/add_points_01.png
new file mode 100644
index 0000000000000000000000000000000000000000..fe85ecd11235cbd22ab6390f50c89ce4414313f8
Binary files /dev/null and b/docs/images/add_points_01.png differ
diff --git a/docs/images/add_points_02.png b/docs/images/add_points_02.png
new file mode 100644
index 0000000000000000000000000000000000000000..8b80c62c0b7094926af0de845aa31d57ec043859
Binary files /dev/null and b/docs/images/add_points_02.png differ
diff --git a/docs/images/add_points_03.png b/docs/images/add_points_03.png
new file mode 100644
index 0000000000000000000000000000000000000000..3269c86b82efd1facd87dbd16bb35c57f1bde165
Binary files /dev/null and b/docs/images/add_points_03.png differ
diff --git a/docs/images/add_points_04.png b/docs/images/add_points_04.png
new file mode 100644
index 0000000000000000000000000000000000000000..91a187bcbf0ea4efa1cb53a1a1f55818dfa2ceac
Binary files /dev/null and b/docs/images/add_points_04.png differ
diff --git a/docs/images/add_points_05.png b/docs/images/add_points_05.png
new file mode 100644
index 0000000000000000000000000000000000000000..a18b8751551d8e6feac67d794961a029682436ae
Binary files /dev/null and b/docs/images/add_points_05.png differ
diff --git a/docs/images/add_points_06.png b/docs/images/add_points_06.png
new file mode 100644
index 0000000000000000000000000000000000000000..1bcdae0ad31100ccdf84cedcd9b179f5b3865a4b
Binary files /dev/null and b/docs/images/add_points_06.png differ
diff --git a/docs/images/add_points_07.png b/docs/images/add_points_07.png
new file mode 100644
index 0000000000000000000000000000000000000000..6eb50386781e370ce4a7a5753cf2cb43ca344edc
Binary files /dev/null and b/docs/images/add_points_07.png differ
diff --git a/docs/images/add_points_08.png b/docs/images/add_points_08.png
new file mode 100644
index 0000000000000000000000000000000000000000..49a2cb40c0ec3e95c1ec8ba5489afc4f9f30fc26
Binary files /dev/null and b/docs/images/add_points_08.png differ
diff --git a/docs/images/add_points_09.png b/docs/images/add_points_09.png
new file mode 100644
index 0000000000000000000000000000000000000000..86c0b087e0c5e81ab79d5281d3602a2d22378b9b
Binary files /dev/null and b/docs/images/add_points_09.png differ
diff --git a/docs/images/edit_points_01.png b/docs/images/edit_points_01.png
new file mode 100644
index 0000000000000000000000000000000000000000..7cf5628af9cf6a4c5191d8a98f9a0254b89a18eb
Binary files /dev/null and b/docs/images/edit_points_01.png differ
diff --git a/docs/images/edit_points_02.png b/docs/images/edit_points_02.png
new file mode 100644
index 0000000000000000000000000000000000000000..c7286461aef4fbd5e04310ef3e61eb6856a00149
Binary files /dev/null and b/docs/images/edit_points_02.png differ
diff --git a/docs/images/edit_points_03.png b/docs/images/edit_points_03.png
new file mode 100644
index 0000000000000000000000000000000000000000..f1ed0fcfefe1994c4cb1a18042690f66bd65b397
Binary files /dev/null and b/docs/images/edit_points_03.png differ
diff --git a/docs/images/edit_points_04.png b/docs/images/edit_points_04.png
new file mode 100644
index 0000000000000000000000000000000000000000..902b142a1cd89d064e831e4f8bc6bf287b436974
Binary files /dev/null and b/docs/images/edit_points_04.png differ
diff --git a/docs/images/edit_points_05.png b/docs/images/edit_points_05.png
new file mode 100644
index 0000000000000000000000000000000000000000..bf8e9a13e45cfb5d687ef57d65fb44e19e352125
Binary files /dev/null and b/docs/images/edit_points_05.png differ
diff --git a/docs/images/open_dataset_00.png b/docs/images/open_dataset_00.png
new file mode 100644
index 0000000000000000000000000000000000000000..91ce794684cf8882d529d60e0721f6c353738ca4
Binary files /dev/null and b/docs/images/open_dataset_00.png differ
diff --git a/docs/images/open_dataset_01.png b/docs/images/open_dataset_01.png
new file mode 100644
index 0000000000000000000000000000000000000000..0cee1e8628360e2b7184fc123b35b8523aa78e90
Binary files /dev/null and b/docs/images/open_dataset_01.png differ
diff --git a/docs/images/open_dataset_02.png b/docs/images/open_dataset_02.png
new file mode 100644
index 0000000000000000000000000000000000000000..5bd54a6fee9ce264939267f82b1a2bdc0a317e82
Binary files /dev/null and b/docs/images/open_dataset_02.png differ
diff --git a/docs/images/open_dataset_03.png b/docs/images/open_dataset_03.png
new file mode 100644
index 0000000000000000000000000000000000000000..2c81fb05f953050905d31dcad52bf13de94fb861
Binary files /dev/null and b/docs/images/open_dataset_03.png differ
diff --git a/docs/images/progress.png b/docs/images/progress.png
new file mode 100644
index 0000000000000000000000000000000000000000..48cd34f74df7326f808ec2e37f9342f63f3d8e73
Binary files /dev/null and b/docs/images/progress.png differ
diff --git a/docs/images/save_dataset_01.png b/docs/images/save_dataset_01.png
new file mode 100644
index 0000000000000000000000000000000000000000..61f04c6335e6637faaf689705f7df0ec80936722
Binary files /dev/null and b/docs/images/save_dataset_01.png differ
diff --git a/docs/images/save_dataset_02.png b/docs/images/save_dataset_02.png
new file mode 100644
index 0000000000000000000000000000000000000000..737bb87f84757b572ef2fe915110a261ee85c339
Binary files /dev/null and b/docs/images/save_dataset_02.png differ
diff --git a/docs/images/save_dataset_03.png b/docs/images/save_dataset_03.png
new file mode 100644
index 0000000000000000000000000000000000000000..e96ad8ac84b314edd33809683dae275a3f626a69
Binary files /dev/null and b/docs/images/save_dataset_03.png differ
diff --git a/docs/images/save_dataset_04.png b/docs/images/save_dataset_04.png
new file mode 100644
index 0000000000000000000000000000000000000000..710a73a7b8c2e8bd156dee1931f48b19aac42ee7
Binary files /dev/null and b/docs/images/save_dataset_04.png differ
diff --git a/docs/styles/default.css b/docs/styles/default.css
new file mode 100644
index 0000000000000000000000000000000000000000..30b1c2b4d111937a9d896a554f889727b09a29e9
--- /dev/null
+++ b/docs/styles/default.css
@@ -0,0 +1,23 @@
+body {
+ font-size: 14pt;
+ font-family: Verdana, Geneva, Arial, Helvetica, sans-serif;
+ color: black;
+ line-height: 14pt;
+ padding-left: 5pt;
+ padding-right: 5pt;
+ padding-top: 5pt;
+}
+
+li {
+ width: 600pt;
+}
+
+p {
+ width: 600pt;
+}
+
+p.author-subheading {
+ font-size: 12pt;
+ font-family: Verdana, Geneva, Arial, Helvetica, sans-serif;
+ width: 600pt;
+}
\ No newline at end of file
diff --git a/docs/user_guide.html b/docs/user_guide.html
new file mode 100644
index 0000000000000000000000000000000000000000..41451d09f6e05ec0da1470e74a532c6203936dd8
--- /dev/null
+++ b/docs/user_guide.html
@@ -0,0 +1,140 @@
+<!DOCTYPE html>
+<html>
+<head>
+ <title>Drift Corrector User Guide</title>
+ <meta charset="utf-8">
+ <link rel="stylesheet" type="text/css" href="styles/default.css">
+</head>
+<body>
+
+<h1>Drift Corrector User Guide</h1>
+
+<div>
+<p class="author-subheading">Written by Hans Jacob Teglbjærg Stephensen</p>
+</div>
+
+<p>This document is a user guide for the Drift Corrector toolbox for performing correctional translation of isotropic images with visible vesicles.</p>
+
+<p>The general process includes</p>
+<ol>
+ <li>Open a dataset to be corrected (we currently support a number image, including multipage tiff images or a numbered image files in various formats such as TIF, PNG, JPG, and BMP).</li>
+ <li>Annotate vesicles individually until the target error becomes low enough found for each image section.</li>
+ <li>Generate a corrected dataset.</li>
+</ol>
+
+<p>The process requires a few manual steps which are covered here in detail.</p>
+
+<h2>Open a dataset</h2>
+
+<ol>
+ <li>
+ <p>To open a dataset, make sure it is in a supported format. If stored as multiple image files, place these in a separate folder with names ending by a number (e.g., my_image_01.png, my_image_02.png, ...) starting at either 0 or 1, numbered consecutively. Additionally, multipage TIF files are also supported.</p>
+ <img src="images/open_dataset_00.png">
+ </li>
+ <li>
+ <p>To open the dataset, click <strong>Open Dataset</strong> in the left pane. When prompted, choose the dataset format.</p>
+ <img src="images/open_dataset_01.png">
+ </li>
+ <li>
+ <p>Navigate to the dataset.</p>
+ <img src="images/open_dataset_02.png">
+ </li>
+ <li>
+ <p>Large datasets may take a moment to load.</p>
+ <img src="images/open_dataset_03.png">
+ </li>
+</ol>
+
+<h2>Save Annotation Progress</h2>
+
+<p>Saving progress is done automatically in a automatically created subdirectory inside the datasets directory and does not require actions by the user.</p>
+<img src="images/progress.png">
+
+<h2>Annotate a vesicle</h2>
+
+<ol>
+ <li>
+ <p>Ensure <strong>Add Points</strong> mode is enabled in the left pane.</p>
+ <img src="images/add_points_01.png">
+ </li>
+ <li>
+ <p>Use mouse scroll to navigate the image sections if needed.</p>
+ <img src="images/add_points_02.png">
+ </li>
+ <li>
+ <p>Right-click a vesicle to pan to it.</p>
+ <img src="images/add_points_03.png">
+ </li>
+ <li>
+ <p>Zoom in a size where accurate annotation is possible.</p>
+ <img src="images/add_points_04.png">
+ </li>
+ <li>
+ <p>Left-click a vesicle boundary to annotate a point.</p>
+ <img src="images/add_points_05.png">
+ </li>
+ <li>
+ <p>Use mouse scroll to navigate the image sections to annotate the vesicle across multiple sections.</p>
+ <img src="images/add_points_06.png">
+ </li>
+
+ <li><p>Click the "a" button to change viewing and annotation axis. When not viewing in the standard image plane, the right pane shows statistics on the predicted accuracy of the drift estimate. In the standard view, this pane shows details about the current section.</p>
+ <img src="images/add_points_07.png">
+ </li>
+ <li><p>To delete a point, ensure the <strong>Edit Points</strong> mode is enabled. Left-click the point to delete it.</p>
+<img src="images/add_points_08.png">
+ </li>
+ <li><p>When done annotating the vesicle, click <strong>Next Vesicle</strong> in the left pane. Notice the changes in the right pane. The confidence interval (green line) and vesicle count (black line) changes to reflect a drift estimate with the new vesicle annotated.</p>
+<img src="images/add_points_09.png">
+ </li>
+</ol>
+
+<h2>Edit an annotated vesicle</h2>
+<ol>
+ <li>
+ <p>Ensure <strong>Edit Points</strong> is selected in the left pane.</p>
+ <img src="images/edit_points_01.png">
+ </li>
+ <li>
+ <p>Right-click a vesicle to activate it for editing points.</p>
+ <img src="images/edit_points_02.png">
+ </li>
+ <li>
+ <p>In <strong>Edit Points</strong> mode, left-click to delete points.</p>
+ <img src="images/edit_points_03.png">
+ </li>
+ <li>
+ <p>In <strong>Add Points</strong> mode, left-click to add points.</p>
+ <img src="images/edit_points_04.png">
+ </li>
+ <li>
+ <p>When done annotating the vesicle, click <strong>Next Vesicle</strong> in the left pane.</p>
+ <img src="images/edit_points_05.png">
+ </li>
+</ol>
+
+<h2>Save corrected dataset</h2>
+<ol>
+ <li>
+ <p>It is recommended first to create a folder to contain the corrected dataset.</p>
+ </li>
+ <li>
+ <p>Click <strong>Save Corrected Dataset</strong> in the left pane. Choose an output format and click <strong>submit</strong>.</p>
+ <img src="images/save_dataset_01.png">
+ </li>
+ <li>
+ <p>Navigate to a target directory for the output. <strong>Warning: Choosing a directory containing other files may overwrite the existing data.</strong>. When done, click <strong><u>O</u>K</strong>.</p>
+ <img src="images/save_dataset_02.png">
+ </li>
+ <li>
+ <p>Wait while the corrected dataset is produced.</p>
+ <img src="images/save_dataset_03.png">
+ </li>
+ <li>
+ <p>The output dataset is extended in size so that no parts of the dataset is transformed outside the view.</p>
+ <img src="images/save_dataset_04.png">
+ </li>
+</ol>
+
+</body>
+</html>
\ No newline at end of file
diff --git a/pil_hook.py b/pil_hook.py
new file mode 100644
index 0000000000000000000000000000000000000000..d0c357690515405d7b127e7b8aa5d314f4664bd5
--- /dev/null
+++ b/pil_hook.py
@@ -0,0 +1 @@
+hiddenimports=['PIL', 'PIL._imagingtk', 'PIL._tkinter_finder']
\ No newline at end of file
diff --git a/src/AskMultiDialog.py b/src/AskMultiDialog.py
new file mode 100644
index 0000000000000000000000000000000000000000..28b33f34ac5da20df932f00068c22dee269f34b7
--- /dev/null
+++ b/src/AskMultiDialog.py
@@ -0,0 +1,22 @@
+from tkinter import Tk, Label, Button, Radiobutton, IntVar, Toplevel
+
+def AskMultiDialog(root, prompt, options):
+
+ window = Toplevel(root)
+
+ def quit_and_destroy():
+ window.quit()
+ window.destroy()
+
+ if prompt:
+ Label(window, text=prompt).grid(row=0, column=0)
+
+ v = IntVar()
+
+ for i, option in enumerate(options):
+ Radiobutton(window, text=option, variable=v, value=i).grid(row=i+1, column=0)
+
+ Button(window, text="Submit", command=quit_and_destroy).grid(row=len(options)+1, column=0)
+ window.mainloop()
+
+ return options[v.get()]
diff --git a/src/DatasetManager.py b/src/DatasetManager.py
new file mode 100644
index 0000000000000000000000000000000000000000..2701bc52a00cb1ab1a23b48d18891e362d0da981
--- /dev/null
+++ b/src/DatasetManager.py
@@ -0,0 +1,181 @@
+import itertools
+from PIL import Image, ImageTk, TiffImagePlugin
+import os
+from collections import Counter
+import numpy as np
+
+class DatasetManager:
+
+ _accepted_formats = ['.png', '.jpg', '.jpeg', '.tif', '.tiff', '.bmp']
+
+ def __init__(self, dataset_directory, dataset_format):
+ self._dataset_directory = os.fsencode(dataset_directory)
+ self.dataset_format = dataset_format
+ self._volume, self.filenames = self._read_dataset()
+ self._dtype = None
+ self._axis = 0
+
+ def __len__(self):
+ return self._volume.shape[0]
+
+ def __getitem__(self, key):
+ return self._volume[key, :, :]
+
+ @property
+ def dtype(self):
+ return self._dtype
+
+ @property
+ def shape(self):
+ return self._volume.shape
+
+ @property
+ def axis(self):
+ return self._axis
+
+ @property
+ def axis_shape(self):
+ return tuple([dim for axis, dim in enumerate(self.shape) if axis != self.axis])
+
+ def set_axis(self, axis):
+ self._axis = axis
+
+ def _read_dataset(self):
+
+ if self.dataset_format == 'multi image':
+ vfilenames = self._get_multi_image_volume_and_filenames()
+ elif self.dataset_format == 'multipage file':
+ vfilenames = self._get_multipage_image_volume_and_filename()
+
+ return vfilenames
+
+ def _get_multi_image_volume_and_filenames(self):
+
+ filtered_string_list = []
+ file_formats_counter = Counter()
+
+ for file in os.listdir(self._dataset_directory):
+ filename = os.fsdecode(file)
+
+ formatted_string_list = [''.join(x) for _, x in itertools.groupby(filename, key=str.isdigit)]
+
+ file_format = formatted_string_list[-1]
+
+ if file_format in self._accepted_formats:
+ file_formats_counter[file_format] += 1
+ filtered_string_list.append(formatted_string_list)
+
+ if len(file_formats_counter) == 0:
+ raise IOError("No valid image files found")
+
+ chosen_ext = file_formats_counter.most_common()[0][0]
+
+ number_filename_list = []
+
+ for string_list in filtered_string_list:
+ ext = string_list[-1]
+ if chosen_ext == ext:
+ idx = int(string_list[-2])
+ filename = ''.join(string_list)
+ full_image_path = os.path.join(os.fsdecode(self._dataset_directory), filename)
+ number_filename_list.append([idx, full_image_path, filename])
+
+ number_filename_list = sorted(number_filename_list)
+ filenames = [filename for _, _, filename in number_filename_list]
+
+ fpairs = zip(number_filename_list[:-1], number_filename_list[1:])
+
+ # Check that sections are consecutive
+ for (idx_a, path_a, filename_a), (idx_b, path_b, filename_b) in fpairs:
+ if idx_b - idx_a != 1:
+ err_msg_info = filename_a + ' with number ' + str(idx_a) \
+ + ' and ' + filename_b + ' with number ' + str(idx_b)
+
+ raise IOError('Files must be numbered consecutively, got: ' + err_msg_info)
+
+ self.min_idx = number_filename_list[0][0]
+
+ volume = np.stack([np.array(Image.open(filepath)) for (N, filepath, _) in number_filename_list])
+ return volume, filenames
+
+ def _get_multipage_image_volume_and_filename(self):
+ img = Image.open(self._dataset_directory)
+ image_list = []
+
+ full_filename = os.fsdecode(self._dataset_directory).split('/')[-1]
+ filename_list = full_filename.split(".")
+ filename = ''.join(filename_list[:-1])
+ ext = filename_list[-1]
+
+ filenames = []
+
+ i = 0
+
+ while True:
+ try:
+ img.seek(i)
+ image_list.append(np.array(img))
+ filenames.append(filename + '_' + str(i) + '.' + ext)
+ i += 1
+ except EOFError:
+ # Not enough frames in img
+ break
+
+ volume = np.stack(image_list)
+
+ return volume, filenames
+
+ def get_image_patch(self, center, shape):
+
+ if self.axis == 0:
+ image_slice = self._volume[center[0], :, :]
+ i, j = center[1], center[2]
+ h, w = self._volume.shape[1], self._volume.shape[2]
+ if self.axis == 1:
+ image_slice = self._volume[:, center[1], :]
+ i, j = center[0], center[2]
+ h, w = self._volume.shape[0], self._volume.shape[2]
+ if self.axis == 2:
+ image_slice = self._volume[:, :, center[2]]
+ i, j = center[0], center[1]
+ h, w = self._volume.shape[0], self._volume.shape[1]
+
+ ph, pw = shape
+
+ if ph % 2 == 0 or pw % 2 == 0:
+ raise ValueError('Shape widths must be odd')
+
+ patch = np.zeros(shape, dtype=self.dtype)
+
+ for si in range(ph):
+ for sj in range(pw):
+ ii = i + si - ph//2
+ jj = j + sj - pw//2
+ if ii >= 0 and ii < h and jj >= 0 and jj < w:
+ patch[si, sj] = image_slice[ii, jj]
+
+ return patch, (j - pw//2, i - ph//2)
+
+
+ def initializeDatasetOutput(self, output_dir, output_shape, dataset_type):
+
+ self.output_dir = output_dir
+
+ #if not os.path.exists(output_dir):
+ # os.makedirs(output_dir)
+
+ if dataset_type == 'multipage file':
+ self.ouput_writer = TiffImagePlugin.AppendingTiffWriter(output_dir,True)
+
+ def saveTransformedImage(self, k, image, dataset_type):
+
+ if dataset_type == 'multi image':
+ filename = self.filenames[k]
+ output_filepath = os.path.join(self.output_dir, filename)
+ output_image = Image.fromarray(image).convert('L')
+ output_image.save(output_filepath)
+
+ elif dataset_type == 'multipage file':
+ im = Image.fromarray(image).convert('L')
+ im.save(self.ouput_writer)
+ self.ouput_writer.newFrame()
diff --git a/src/DriftEstimate.py b/src/DriftEstimate.py
new file mode 100644
index 0000000000000000000000000000000000000000..8612694c62019f63c83408ce25214a53eba85b93
--- /dev/null
+++ b/src/DriftEstimate.py
@@ -0,0 +1,46 @@
+import numpy as np
+
+def driftEstimate(ellipsoids, z_bounds, drift_spread=7, minimum_N=2):
+
+ drift_parameters = []
+
+ for E in ellipsoids:
+
+ v, radii, center = E
+
+ A, B, C, D, E, F = v
+
+ sx = (D*F - B*E) / (A*B - D**2)
+ sy = (D*E - A*F) / (A*B - D**2)
+
+ drift_parameters.append([sx, sy, center[0], center[1], center[2]])
+
+ z_lower, z_upper = z_bounds
+
+ drift_bins = [[] for i in range(z_lower, z_upper)]
+
+ drift_components = []
+
+ for param in drift_parameters:
+
+ drift_component = param[:2]
+ center_z = param[4]
+ center_k = int(round(center_z))
+ drift_components.append(drift_component)
+
+ for k in range(center_k - drift_spread, center_k + drift_spread + 1):
+ if k >= z_lower and k < z_upper:
+ k_array = k - z_lower
+ drift_bins[k_array].append(drift_component)
+
+ drift_bins_means = [np.mean(np.array(bin_group), axis=0) if len(bin_group) >= minimum_N else np.array([0, 0]) for bin_group in drift_bins]
+ drift_bins_std = [np.std(np.array(bin_group), axis=0) if len(bin_group) >= minimum_N else np.array([-1, -1]) for bin_group in drift_bins]
+ drift_bins_N = [len(bin_group) for bin_group in drift_bins]
+
+ # Interpolation, filling missing values etc., here
+
+ drift_estimate = np.array(drift_bins_means)
+ drift_std = np.array(drift_bins_std)
+ drift_bins_N = np.array(drift_bins_N)
+
+ return drift_estimate, drift_std, drift_bins_N, drift_components
\ No newline at end of file
diff --git a/src/DriftEstimate2.py b/src/DriftEstimate2.py
new file mode 100644
index 0000000000000000000000000000000000000000..b40b679432b25383e7162cb7973e331ca08260b3
--- /dev/null
+++ b/src/DriftEstimate2.py
@@ -0,0 +1,300 @@
+import numpy as np
+from EllipsoidFit import ellipsoidFit
+import pickle
+from scipy.interpolate import interp1d
+from scipy.optimize import minimize
+from tqdm import tqdm
+import matplotlib.pyplot as plt
+from time import sleep
+
+def cumulative_interp(k):
+ c = np.cumsum(k, axis=0)
+ c_interp_0 = interp1d(range(0, len(k)), c[:,0], kind='linear', fill_value='extrapolate')
+ c_interp_1 = interp1d(range(0, len(k)), c[:,1], kind='linear', fill_value='extrapolate')
+ return (c_interp_0, c_interp_1)
+
+def calculate_drift(E):
+ v, radii, center = E
+
+ A, B, C, D, E, F = v
+
+ sx = (D*F - B*E) / (A*B - D**2)
+ sy = (D*E - A*F) / (A*B - D**2)
+
+ return sx, sy
+
+def driftEstimate(ellipsoids, z_bounds, drift_spread=7, minimum_N=2):
+
+ drift_parameters = []
+
+ for E in ellipsoids:
+
+ v, radii, center = E
+
+ A, B, C, D, E, F = v
+
+ sx = (D*F - B*E) / (A*B - D**2)
+ sy = (D*E - A*F) / (A*B - D**2)
+
+ drift_parameters.append([sx, sy, center[0], center[1], center[2]])
+
+ z_lower, z_upper = z_bounds
+
+ drift_bins = [[] for i in range(z_lower, z_upper)]
+
+ drift_components = []
+
+ for param in drift_parameters:
+
+ drift_component = param[:2]
+ center_z = param[4]
+ center_k = int(round(center_z))
+ drift_components.append(drift_component)
+
+ for k in range(center_k - drift_spread, center_k + drift_spread + 1):
+ if k >= z_lower and k < z_upper:
+ k_array = k - z_lower
+ drift_bins[k_array].append(drift_component)
+
+ drift_bins_means = [np.mean(np.array(bin_group), axis=0) if len(bin_group) >= minimum_N else np.array([0, 0]) for bin_group in drift_bins]
+ drift_bins_std = [np.std(np.array(bin_group), axis=0) if len(bin_group) >= minimum_N else np.array([-1, -1]) for bin_group in drift_bins]
+ drift_bins_N = [len(bin_group) for bin_group in drift_bins]
+
+ # Interpolation, filling missing values etc., here
+
+ drift_estimate = np.array(drift_bins_means)
+ drift_std = np.array(drift_bins_std)
+ drift_bins_N = np.array(drift_bins_N)
+
+ return drift_estimate, drift_std, drift_bins_N, drift_components
+
+def drift_points(points, c_interp):
+
+ new_points = []
+
+ c_interp_0, c_interp_1 = c_interp
+
+ for X in points:
+ new_X = np.copy(X)
+ new_X[:,0] += c_interp_0(X[:,2])
+ new_X[:,1] += c_interp_1(X[:,2])
+ new_points.append(new_X)
+
+ return new_points
+
+def driftEstimate2(points, z_bounds):
+
+ drift_parameters = []
+
+ n = z_bounds[1] - z_bounds[0]
+
+ k = np.zeros((2*(n-1), 1), dtype=np.float32)
+
+ def callback(xk):
+ print(xk[:10:2])
+
+ def energy(x):
+ k = np.empty((n, 2), dtype=np.float32)
+ k[0,0] = 0
+ k[0,1] = 0
+ k[1:,:] = x.reshape((n-1, 2))
+ c_interp = cumulative_interp(k)
+ new_points = drift_points(points, c_interp)
+
+ energy = 0
+
+ for X in new_points:
+ if X.shape[0] >= 10:
+ E = ellipsoidFit(X)
+ sx, sy = calculate_drift(E)
+ print(sx, sy)
+ energy += abs(sx)**2 + abs(sy)**2
+
+ return energy
+
+ res = minimize(energy, k, callback=callback)
+
+ x = np.zeros(2*n, dtype=np.float32)
+ x[2:] = res.x
+
+ return x.reshape((n, 2))
+
+def driftEstimate3(points, z_bounds):
+
+ drift_parameters = []
+
+ n = z_bounds[1] - z_bounds[0]
+
+ k = np.zeros((n, 2), dtype=np.float32)
+
+ n_passes = 1
+
+ for i in tqdm(range(n_passes), desc='n passes'):
+
+ print('k before:', k)
+
+ for j in tqdm(range(1, n), desc='slices'):
+
+ def energy(x):
+
+ k_in = np.copy(k)
+ k_in[j,:] = x
+ c_interp = cumulative_interp(k_in)
+ new_points = drift_points(points, c_interp)
+
+ energy = 0
+
+ for X in new_points:
+ if X.shape[0] >= 10:
+ E = ellipsoidFit(X)
+ sx, sy = calculate_drift(E)
+ energy += sx**2 + sy**2
+
+ return energy
+
+ #plt.plot(np.linspace(-0.4, 0.4, 15), [energy([x,0]) for x in np.linspace(-0.4, 0.4, 15)])
+ #plt.show()
+
+ res = minimize(energy, k[j,:])
+ k[j,:] = res.x
+ print('k after',i, ':', k)
+
+ return k
+
+if __name__ == '__main__':
+
+ if False:
+ import matplotlib.pyplot as plt
+
+ data_path = '/home/dith/arngorf/code/fibsem_tools/data/image_transforms/synthetic_linear/'
+ transforms_path = data_path + 'transforms.txt'
+ points_path = data_path + 'vesicle_ellipsoid_manual_points_fixed.txt'
+
+ true_drift = []
+
+ with open(transforms_path, 'r') as f:
+ for line in f.readlines():
+ formated_line = list(map(float, line.rstrip().split(' ')))
+ true_drift.append([formated_line[0], formated_line[1]])
+
+ true_drift = np.array(true_drift)
+
+ vesicle_points = []
+
+ cur_vidx = -1
+ with open(points_path, 'r') as f:
+ for line in f.readlines():
+ fline = line.rstrip().split(' ')
+ vidx, x, y, z = int(fline[0]), float(fline[1]), float(fline[2]), float(fline[3])
+ if vidx != cur_vidx:
+ cur_vidx = vidx
+ if len(vesicle_points) > 0:
+ vesicle_points[-1] = np.array(vesicle_points[-1])
+ vesicle_points.append([])
+
+ vesicle_points[-1].append([x,y,z])
+
+ vesicle_points[-1] = np.array(vesicle_points[-1])
+
+ ellipsoids = []
+
+ for X in vesicle_points:
+ if X.shape[0] >= 10:
+ E = ellipsoidFit(X)
+ ellipsoids.append(E)
+
+ results = driftEstimate(ellipsoids, (0, 200))
+ results2 = driftEstimate3(vesicle_points, (0, 200))
+
+ plt.subplot(2,1,1)
+ plt.plot(true_drift[:,0], label='true', color='black', ls='--')
+ plt.plot(results[0][:,1], label='original')
+ plt.plot(results2[:,1], label='new')
+ plt.xlim(0,199)
+
+ plt.subplot(2,1,2)
+ plt.plot(true_drift[:,1], label='true', color='black', ls='--')
+ plt.plot(results[0][:,0], label='original')
+ plt.plot(results2[:,0], label='new')
+ plt.xlim(0,199)
+
+ plt.show()
+
+ sx = 0.25
+ sy = 0.0
+
+ a = 5
+ b = 5
+ c = 5
+
+ vesicle_test_points = []
+
+ for vidx in range(10):
+ vesicle_test_points.append([])
+
+ for u in np.linspace(0.1234, 2*np.pi, 9)[:-1]:
+ for v in np.linspace(0.1234, np.pi, 9)[:-1]:
+ x = a * np.cos(u) * np.sin(v) + np.random.uniform(-0.2, 0.2)
+ y = b * np.sin(u) * np.sin(v) + np.random.uniform(-0.2, 0.2)
+ z = c * np.cos(v) + np.random.uniform(-0.2, 0.2)
+
+ z += 3 + vidx
+
+ x += sx * z
+ y += sy * z
+
+ vesicle_test_points[-1].append([x,y,z])
+
+ vesicle_test_points[-1] = np.array(vesicle_test_points[-1])
+
+ X = np.round(vesicle_test_points[0],4)
+ print(X)
+ E = ellipsoidFit(X)
+ E_X_sx, E_X_sy = calculate_drift(E)
+
+ for i in range(16):
+
+ plt.subplot(4,4,i+1)
+
+ k = np.ones((10,2))
+ k[:,0] = -2 + i*0.25
+ k[:,1] = 0
+ #k[4,0] = 2
+ c_interp = cumulative_interp(k)
+ new_points = drift_points([X], c_interp)
+
+ plt.scatter(X[:,0], X[:,2], label='before, sx = ' + str(E_X_sx), marker='o')
+
+ Y = new_points[-1]
+ print(Y)
+ E = ellipsoidFit(Y)
+ E_Y_sx, E_Y_sy = calculate_drift(E)
+ plt.scatter(Y[:,0], Y[:,2], label='after, sx = ' + str(E_Y_sx), marker='x')
+
+ plt.legend()
+ plt.show()
+
+ sxs = np.linspace(-3, 3, 1000)
+ sxs_fit_sq = []
+
+ for sx in sxs:
+
+ k = np.ones((10,2))
+ k[:,0] = sx
+ k[:,1] = 0
+ #k[4,0] = 2
+ c_interp = cumulative_interp(k)
+ new_points = drift_points([X], c_interp)
+
+ Y = new_points[-1]
+ E = ellipsoidFit(Y)
+ E_Y_sx, E_Y_sy = calculate_drift(E)
+
+ sxs_fit_sq.append(E_Y_sx**2)
+
+ plt.plot(sxs, sxs_fit_sq)
+ plt.show()
+
+ r = driftEstimate3(vesicle_test_points, (0, 6))
+
+ print(np.round(r,2))
\ No newline at end of file
diff --git a/src/EllipsoidFit.py b/src/EllipsoidFit.py
new file mode 100644
index 0000000000000000000000000000000000000000..ac926adebd3b2e07032a3d35e487ebf29c2404db
--- /dev/null
+++ b/src/EllipsoidFit.py
@@ -0,0 +1,63 @@
+import numpy as np
+
+def ellipsoidFit(X, dtype=np.float32):
+
+ X = np.array(X)
+ N = X.shape[0]
+
+ mu = np.mean(X, axis=0)
+
+ for i in range(N):
+ X[i,:] = X[i,:] - mu
+
+ if N < 9:
+ raise ValueError('Need at least 9 points to fit an ellipsoid')
+
+ x = X[:,0]
+ y = X[:,1]
+ z = X[:,2]
+
+ A = np.empty((N,9), dtype=dtype)
+
+ A[:,0] = x**2 + y**2 - 2*z*z
+ A[:,1] = x**2 + z**2 - 2*y*y
+ A[:,2] = 2*x*y
+ A[:,3] = 2*x*z
+ A[:,4] = 2*y*z
+ A[:,5] = 2*x
+ A[:,6] = 2*y
+ A[:,7] = 2*z
+ A[:,8] = 1
+
+ b = x**2 + y**2 + z**2
+
+ u = np.linalg.solve(np.matmul(A.T, A), np.matmul(A.T, b))
+
+ v = np.zeros(10)
+
+ v[0] = u[0] + u[1] - 1
+ v[1] = u[0] - 2*u[1] - 1
+ v[2] = u[1] - 2*u[0] - 1
+ v[3] = u[2]
+ v[4] = u[3]
+ v[5] = u[4]
+ v[6] = u[5]
+ v[7] = u[6]
+ v[8] = u[7]
+ v[9] = u[8]
+
+ B = np.empty((4,4), dtype=dtype);
+ B = np.array([[v[0], v[3], v[4], v[6]],
+ [v[3], v[1], v[5], v[7]],
+ [v[4], v[5], v[2], v[8]],
+ [v[6], v[7], v[8], v[9]]])
+
+ center = np.linalg.solve(-B[:3,:3], B[:3,3])
+
+ center = center + mu
+
+ evals, evecs = np.linalg.eig(B[:3,:3]*(1./v[9]))
+
+ radii = np.sqrt(1./ abs(evals))
+
+ return v[:6]*(1./v[9]), radii, center
diff --git a/src/ImageWarp.py b/src/ImageWarp.py
new file mode 100644
index 0000000000000000000000000000000000000000..885ad68cbe0e913dc42d06447d1c3b5f6dfc4de6
--- /dev/null
+++ b/src/ImageWarp.py
@@ -0,0 +1,43 @@
+from skimage.transform import SimilarityTransform, warp
+import numpy as np
+import os
+
+def warpAndSave(driftEstimates, dataset, output_dir, dataset_type, interpolation_order=5, double_view_debug=False):
+
+ driftEstimates[0,:] = 0 # correcting wrt. first image, which is not drifted wrt. anything
+
+ cummulative_correction = np.cumsum(driftEstimates, axis=0)
+
+ # calculate new dataset bounds due to drifting of the sections
+ minx = -int(np.floor(np.min(-cummulative_correction[:,0])))
+ maxx = int(np.ceil( np.max(-cummulative_correction[:,0])))
+ miny = -int(np.floor(np.min(-cummulative_correction[:,1])))
+ maxy = int(np.ceil( np.max(-cummulative_correction[:,1])))
+
+ N = len(dataset)
+ shape = dataset.shape
+
+ h = shape[1] + (miny + maxy)
+ if double_view_debug:
+ w = (shape[2] + (minx + maxx))*2
+ else:
+ w = shape[2] + (minx + maxx)
+
+ output_shape = (h, w)
+
+ add_image = dataset.initializeDatasetOutput(output_dir, output_shape, dataset_type)
+
+ for k in range(N):
+
+ result = np.zeros((h, w), dtype=dataset.dtype)
+ result[miny:miny+shape[1], minx:minx+shape[2]] = dataset[k]
+
+ translation = cummulative_correction[k,:]
+
+ tform = SimilarityTransform(scale=1, rotation=0, translation=translation)
+ result = warp(result, tform, order=5)
+
+ if double_view_debug:
+ result[miny:miny+shape[1], (minx + w//2):(minx+shape[2] + w//2)] = dataset[k]
+
+ dataset.saveTransformedImage(k, result.astype(dataset.dtype), dataset_type)
diff --git a/src/VesiclePointManager.py b/src/VesiclePointManager.py
new file mode 100644
index 0000000000000000000000000000000000000000..fe9ca25ba1a41fd2da58f4a7268387c422ae1ef1
--- /dev/null
+++ b/src/VesiclePointManager.py
@@ -0,0 +1,103 @@
+from collections import defaultdict
+import numpy as np
+
+to_key = lambda val: int(round(val))
+
+def p_dist(p1, p2):
+ return np.linalg.norm(np.array(p1) - np.array(p2))
+
+class VesiclePointManager:
+
+ def __init__(self):
+ self._points = {0: defaultdict(list),
+ 1: defaultdict(list),
+ 2: defaultdict(list)}
+
+ self._points_vidx = defaultdict(list)
+
+ self._vidx = 0
+
+ @property
+ def vidx(self):
+ return self._vidx
+
+ def new_vesicle(self):
+ if len(self._points_vidx) > 0:
+ self._vidx = max(self._points_vidx.keys()) + 1
+ else:
+ self._vidx = 0
+
+ def add_point(self, x, y, z):
+
+ #for p in self._points_vidx[self._vidx]:
+ # if p_dist(p, [x,y,z]) > 40:
+ # print('Point too far from existing points', p, [x,y,z])
+ # return
+
+ self._points[0][to_key(z)].append((x,y,self._vidx,z))
+ self._points[1][to_key(y)].append((x,z,self._vidx,y))
+ self._points[2][to_key(x)].append((y,z,self._vidx,x))
+ self._points_vidx[self._vidx].append((x,y,z))
+
+ def delete_nearest_point(self, x, y, z):
+
+ min_dist = 3
+ min_p = None
+ min_vidx = None
+
+ for key in self._points_vidx.keys():
+ for p in self._points_vidx[key]:
+ dist = p_dist(p, [x,y,z])
+ if dist < min_dist:
+ min_dist = dist
+ min_p = p
+ min_vidx = key
+
+ if min_p == None:
+ print('no point was close enough')
+ else:
+ print('found:', min_p, min_vidx, min_dist)
+ x, y, z = min_p
+
+ try:
+ idx = self._points_vidx[min_vidx].index((x,y,z))
+ idx0 = self._points[0][to_key(z)].index((x,y,min_vidx,z))
+ idx1 = self._points[1][to_key(y)].index((x,z,min_vidx,y))
+ idx2 = self._points[2][to_key(x)].index((y,z,min_vidx,x))
+ del self._points_vidx[min_vidx][idx]
+ del self._points[0][to_key(z)][idx0]
+ del self._points[1][to_key(y)][idx1]
+ del self._points[2][to_key(x)][idx2]
+ except ValueError:
+ pass
+
+ def set_to_nearest_vesicle(self, x, y, z):
+ min_dist = 3
+ min_p = None
+ min_vidx = None
+
+ for key in self._points_vidx.keys():
+ for p in self._points_vidx[key]:
+ dist = p_dist(p, [x,y,z])
+ if dist < min_dist:
+ min_dist = dist
+ min_p = p
+ min_vidx = key
+
+ if min_vidx != None:
+ self._vidx = min_vidx
+
+ def get_vesicle_points(self):
+ return [np.array(self._points_vidx[key]) for key in self._points_vidx.keys()]
+
+ def get_vesicle_points_with_vidx(self):
+ points = []
+
+ for key in self._points_vidx.keys():
+ for x, y, z in self._points_vidx[key]:
+ points.append([key, x, y, z])
+
+ return points
+
+ def get_points_in_section(self, idx, axis):
+ return self._points[axis][idx]
\ No newline at end of file