Skip to content
GitLab
Commit d98b1a34 authored by Henri Payno's avatar Henri Payno
Browse files

try to generalize function for z stitching

parent e68195c6
Hide whitespace changes
Inline Side-by-side
doc/tutorials/notebooks/stitching_postprocessing.ipynb 0 → 100644
View file @ d98b1a34
{
"cells": [
{
"cell_type": "markdown",
"id": "6ac304ef-29eb-442a-8b37-1160f64f3f61",
"metadata": {
"tags": []
},
"source": [
"# stitching volume (post processing)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "381ece65-d7be-4674-9673-4999fd1df888",
"metadata": {},
"outputs": [],
"source": [
"from tomoscan.serie import Serie\n",
"from tomoscan.esrf.volume.hdf5volume import HDF5Volume"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "df9cefe7-d53b-4490-b298-906dc5ecb07f",
"metadata": {},
"outputs": [],
"source": [
"%pylab"
]
},
{
"cell_type": "markdown",
"id": "f5070d72-f5cf-49c2-89df-57270d0cf125",
"metadata": {},
"source": [
"### define volume and serie of volume"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b787eef8-608d-4c49-afab-aff6fd0f989d",
"metadata": {},
"outputs": [],
"source": [
"import os\n",
"data_folder = \"/tmp_14_days/payno/test_stitching/\"\n",
"volume_0 = HDF5Volume(\n",
" file_path=os.path.join(data_folder, \"SYS01S101_SYS01S101_3scans_0000_rec.hdf5\"),\n",
" data_path=\"/entry0000/reconstruction\",\n",
")\n",
"\n",
"volume_1 = HDF5Volume(\n",
" file_path=os.path.join(data_folder, \"SYS01S101_SYS01S101_3scans_0001_rec.hdf5\"),\n",
" data_path=\"/entry0001/reconstruction\",\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2ea13719-37fe-4bf2-8207-ed6230ad753c",
"metadata": {},
"outputs": [],
"source": [
"volume_z_serie = Serie([volume_0, volume_1])"
]
},
{
"cell_type": "markdown",
"id": "27f94805-8162-4c7d-ad7f-dfe39e8df5db",
"metadata": {},
"source": [
"## find volumes relative shift"
]
},
{
"cell_type": "markdown",
"id": "1f9dcea3-3196-4e03-a32c-454b549467be",
"metadata": {},
"source": [
"### z shift and z overlap"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b6514810-ebee-4a33-9184-08e6396ca010",
"metadata": {},
"outputs": [],
"source": [
"middle_XZ_slice_idx = 1059\n",
"vol_0_middle_XZ_slice = volume_0.get_slice(xz=middle_XZ_slice_idx)\n",
"vol_1_middle_XZ_slice = volume_1.get_slice(xz=middle_XZ_slice_idx)"
]
},
{
"cell_type": "markdown",
"id": "de777ae9-560f-47b7-91bc-bfe07ae0a56b",
"metadata": {
"tags": []
},
"source": [
"### plot raw volume 0 and volume 1 and plot composite image"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d4291ae4-dcf7-4feb-8b49-93389f310d16",
"metadata": {},
"outputs": [],
"source": [
"import numpy\n",
"def get_composite_image(img1, img2):\n",
" return img1.astype(numpy.float64) - img2.astype(numpy.float64)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "dd161bb6-41e0-49a0-ab09-e42c6f89f47e",
"metadata": {},
"outputs": [],
"source": [
"from matplotlib import pyplot as plt\n",
"fig = plt.figure()\n",
"plt.subplot(131)\n",
"plt.title(\"raw vol 0 middle XZ slice\")\n",
"plt.imshow(vol_0_middle_XZ_slice)\n",
"plt.subplot(132)\n",
"plt.title(\"raw vol 1 middle XZ slice\")\n",
"plt.imshow(vol_1_middle_XZ_slice)\n",
"plt.subplot(133)\n",
"plt.title(\"composite of middle XZ slice\")\n",
"plt.imshow(\n",
" get_composite_image(vol_0_middle_XZ_slice, vol_1_middle_XZ_slice)\n",
")\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "45079206-1af4-4243-b1a2-e6962c65dcad",
"metadata": {},
"source": [
"### find z overlap area"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "33506120-9507-4196-8b46-ee2d98160d89",
"metadata": {},
"outputs": [],
"source": [
"from skimage.registration import phase_cross_correlation\n",
"\n",
"shift, error, diffphase = phase_cross_correlation(\n",
" reference_image=vol_0_middle_XZ_slice,\n",
" moving_image=vol_1_middle_XZ_slice,\n",
" space='real',\n",
")\n",
"print(\"shift found is\", shift)"
]
},
{
"cell_type": "markdown",
"id": "676d1065-701a-48d8-b438-cc0347c64359",
"metadata": {},
"source": [
"### move image 2 to the found shift"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9f6500f7-5dfa-46ca-847d-3cf5a4eff777",
"metadata": {},
"outputs": [],
"source": [
"from scipy.ndimage import shift as shift_scipy\n",
"# here the y shift must be improved.\n",
"vol_1_middle_XZ_slice_shifted = shift_scipy(\n",
" vol_1_middle_XZ_slice,\n",
" mode=\"nearest\",\n",
" shift=-shift,\n",
")\n",
"vol_1_middle_XZ_slice_shifted = vol_1_middle_XZ_slice_shifted[-int(shift[0]):]\n",
"stitching_height_in_px = 300\n",
"# what to do from the x shift ? For now I think this is better to ignore it"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9e4c7c7c-9585-40cb-95cc-a5f88db9650b",
"metadata": {},
"outputs": [],
"source": [
"from matplotlib import pyplot as plt\n",
"fig = plt.figure()\n",
"overlap_vol_0 = vol_0_middle_XZ_slice[:y_stitching_in_px]\n",
"overlap_vol_1 = vol_1_middle_XZ_slice_shifted[-y_stitching_in_px:]\n",
"plt.subplot(131)\n",
"plt.title(\"overlap XZ slice vol 0\")\n",
"plt.imshow(overlap_vol_0)\n",
"plt.subplot(132)\n",
"plt.title(\"overlap XZ slice vol 1\")\n",
"plt.imshow(overlap_vol_1)\n",
"plt.subplot(133)\n",
"plt.title(\"composite of middle XZ slice with second shifted\")\n",
"plt.imshow(\n",
" get_composite_image(overlap_vol_0, overlap_vol_1)\n",
")"
]
},
{
"cell_type": "markdown",
"id": "f4ee5ec2-122b-41ad-8a2e-6a6a0bbec3f0",
"metadata": {},
"source": [
"## Check overlap volume\n",
"see result of the possible algorithm to compute the overlap area"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1d165c3a-3a95-4ac1-8035-14919306bcc8",
"metadata": {},
"outputs": [],
"source": [
"from tomwer.core.utils.stitching import z_stitch_frames"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c460073d-2005-4232-bbd9-913b5af1e0c9",
"metadata": {},
"outputs": [],
"source": [
"# TODO: rename stich projection to stitch_frames\n",
"stitched_projection_lin_weights, stitched_overlap_lin_weights = z_stitch_frames(\n",
" vol_0_middle_XZ_slice,\n",
" vol_1_middle_XZ_slice_shifted,\n",
" stitching_height_in_px=stitching_height_in_px,\n",
" strategy=\"linear weights\",\n",
")[0:2]"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a6097056-aa4a-4df0-9599-fe6057fd7f6d",
"metadata": {},
"outputs": [],
"source": [
"stitched_projection_cos_weights, stitched_overlap_cos_weights = z_stitch_frames(\n",
" vol_0_middle_XZ_slice,\n",
" vol_1_middle_XZ_slice_shifted,\n",
" stitching_height_in_px=stitching_height_in_px,\n",
" strategy=\"cosinus weights\",\n",
")[0:2]"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9249d6fd-e47b-4ff6-a3c9-a6ee457c2e05",
"metadata": {},
"outputs": [],
"source": [
"stitched_projection_mean, stitched_overlap_mean = z_stitch_frames(\n",
" vol_0_middle_XZ_slice,\n",
" vol_1_middle_XZ_slice_shifted,\n",
" stitching_height_in_px=stitching_height_in_px,\n",
" strategy=\"mean\",\n",
")[0:2]"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "57f72ff6-d469-442f-97cb-3e13a51dd8ee",
"metadata": {},
"outputs": [],
"source": [
"plt.figure()\n",
"plt.subplot(131)\n",
"plt.title(\"mean\")\n",
"plt.imshow(stitched_overlap_mean)\n",
"plt.subplot(132)\n",
"plt.title(\"linear weight\")\n",
"plt.imshow(stitched_overlap_lin_weights)\n",
"plt.subplot(133)\n",
"plt.title(\"cosinus weight\")\n",
"plt.imshow(stitched_overlap_cos_weights)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "999756a9-4ccc-49df-920b-f387d8e51de8",
"metadata": {},
"outputs": [],
"source": [
"plt.figure()\n",
"plt.subplot(131)\n",
"plt.title(\"mean\")\n",
"plt.imshow(stitched_projection_mean)\n",
"plt.subplot(132)\n",
"plt.title(\"linear weight\")\n",
"plt.imshow(stitched_projection_lin_weights)\n",
"plt.subplot(133)\n",
"plt.title(\"cosinus weight\")\n",
"plt.imshow(stitched_projection_cos_weights)"
]
},
{
"cell_type": "markdown",
"id": "bcdc8433-0018-47f7-b094-3cc5d0f7c9e8",
"metadata": {},
"source": [
"## compute final stitched volume"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "71a02b15-1b02-4aa4-8a81-4f36c54134a2",
"metadata": {},
"outputs": [],
"source": [
"# TODO"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.10"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
%% Cell type:markdown id:6ac304ef-29eb-442a-8b37-1160f64f3f61 tags:
# stitching volume (post processing)
%% Cell type:code id:381ece65-d7be-4674-9673-4999fd1df888 tags:
``` python
from tomoscan.serie import Serie
from tomoscan.esrf.volume.hdf5volume import HDF5Volume
```
%% Cell type:code id:df9cefe7-d53b-4490-b298-906dc5ecb07f tags:
``` python
%pylab
```
%% Cell type:markdown id:f5070d72-f5cf-49c2-89df-57270d0cf125 tags:
### define volume and serie of volume
%% Cell type:code id:b787eef8-608d-4c49-afab-aff6fd0f989d tags:
``` python
import os
data_folder = "/tmp_14_days/payno/test_stitching/"
volume_0 = HDF5Volume(
file_path=os.path.join(data_folder, "SYS01S101_SYS01S101_3scans_0000_rec.hdf5"),
data_path="/entry0000/reconstruction",
)
volume_1 = HDF5Volume(
file_path=os.path.join(data_folder, "SYS01S101_SYS01S101_3scans_0001_rec.hdf5"),
data_path="/entry0001/reconstruction",
)
```
%% Cell type:code id:2ea13719-37fe-4bf2-8207-ed6230ad753c tags:
``` python
volume_z_serie = Serie([volume_0, volume_1])
```
%% Cell type:markdown id:27f94805-8162-4c7d-ad7f-dfe39e8df5db tags:
## find volumes relative shift
%% Cell type:markdown id:1f9dcea3-3196-4e03-a32c-454b549467be tags:
### z shift and z overlap
%% Cell type:code id:b6514810-ebee-4a33-9184-08e6396ca010 tags:
``` python
middle_XZ_slice_idx = 1059
vol_0_middle_XZ_slice = volume_0.get_slice(xz=middle_XZ_slice_idx)
vol_1_middle_XZ_slice = volume_1.get_slice(xz=middle_XZ_slice_idx)
```
%% Cell type:markdown id:de777ae9-560f-47b7-91bc-bfe07ae0a56b tags:
### plot raw volume 0 and volume 1 and plot composite image
%% Cell type:code id:d4291ae4-dcf7-4feb-8b49-93389f310d16 tags:
``` python
import numpy
def get_composite_image(img1, img2):
return img1.astype(numpy.float64) - img2.astype(numpy.float64)
```
%% Cell type:code id:dd161bb6-41e0-49a0-ab09-e42c6f89f47e tags:
``` python
from matplotlib import pyplot as plt
fig = plt.figure()
plt.subplot(131)
plt.title("raw vol 0 middle XZ slice")
plt.imshow(vol_0_middle_XZ_slice)
plt.subplot(132)
plt.title("raw vol 1 middle XZ slice")
plt.imshow(vol_1_middle_XZ_slice)
plt.subplot(133)
plt.title("composite of middle XZ slice")
plt.imshow(
get_composite_image(vol_0_middle_XZ_slice, vol_1_middle_XZ_slice)
)
plt.show()
```
%% Cell type:markdown id:45079206-1af4-4243-b1a2-e6962c65dcad tags:
### find z overlap area
%% Cell type:code id:33506120-9507-4196-8b46-ee2d98160d89 tags:
``` python
from skimage.registration import phase_cross_correlation
shift, error, diffphase = phase_cross_correlation(
reference_image=vol_0_middle_XZ_slice,
moving_image=vol_1_middle_XZ_slice,
space='real',
)
print("shift found is", shift)
```
%% Cell type:markdown id:676d1065-701a-48d8-b438-cc0347c64359 tags:
### move image 2 to the found shift
%% Cell type:code id:9f6500f7-5dfa-46ca-847d-3cf5a4eff777 tags:
``` python
from scipy.ndimage import shift as shift_scipy
# here the y shift must be improved.
vol_1_middle_XZ_slice_shifted = shift_scipy(
vol_1_middle_XZ_slice,
mode="nearest",
shift=-shift,
)
vol_1_middle_XZ_slice_shifted = vol_1_middle_XZ_slice_shifted[-int(shift[0]):]
stitching_height_in_px = 300
# what to do from the x shift ? For now I think this is better to ignore it
```
%% Cell type:code id:9e4c7c7c-9585-40cb-95cc-a5f88db9650b tags:
``` python
from matplotlib import pyplot as plt
fig = plt.figure()
overlap_vol_0 = vol_0_middle_XZ_slice[:y_stitching_in_px]
overlap_vol_1 = vol_1_middle_XZ_slice_shifted[-y_stitching_in_px:]
plt.subplot(131)
plt.title("overlap XZ slice vol 0")
plt.imshow(overlap_vol_0)
plt.subplot(132)
plt.title("overlap XZ slice vol 1")
plt.imshow(overlap_vol_1)
plt.subplot(133)
plt.title("composite of middle XZ slice with second shifted")
plt.imshow(
get_composite_image(overlap_vol_0, overlap_vol_1)
)
```
%% Cell type:markdown id:f4ee5ec2-122b-41ad-8a2e-6a6a0bbec3f0 tags:
## Check overlap volume
see result of the possible algorithm to compute the overlap area
%% Cell type:code id:1d165c3a-3a95-4ac1-8035-14919306bcc8 tags:
``` python
from tomwer.core.utils.stitching import z_stitch_frames
```
%% Cell type:code id:c460073d-2005-4232-bbd9-913b5af1e0c9 tags:
``` python
# TODO: rename stich projection to stitch_frames
stitched_projection_lin_weights, stitched_overlap_lin_weights = z_stitch_frames(
vol_0_middle_XZ_slice,
vol_1_middle_XZ_slice_shifted,
stitching_height_in_px=stitching_height_in_px,
strategy="linear weights",
)[0:2]
```
%% Cell type:code id:a6097056-aa4a-4df0-9599-fe6057fd7f6d tags:
``` python
stitched_projection_cos_weights, stitched_overlap_cos_weights = z_stitch_frames(
vol_0_middle_XZ_slice,
vol_1_middle_XZ_slice_shifted,
stitching_height_in_px=stitching_height_in_px,
strategy="cosinus weights",
)[0:2]
```
%% Cell type:code id:9249d6fd-e47b-4ff6-a3c9-a6ee457c2e05 tags:
``` python
stitched_projection_mean, stitched_overlap_mean = z_stitch_frames(
vol_0_middle_XZ_slice,
vol_1_middle_XZ_slice_shifted,
stitching_height_in_px=stitching_height_in_px,
strategy="mean",
)[0:2]
```
%% Cell type:code id:57f72ff6-d469-442f-97cb-3e13a51dd8ee tags:
``` python
plt.figure()
plt.subplot(131)
plt.title("mean")
plt.imshow(stitched_overlap_mean)
plt.subplot(132)
plt.title("linear weight")
plt.imshow(stitched_overlap_lin_weights)
plt.subplot(133)
plt.title("cosinus weight")
plt.imshow(stitched_overlap_cos_weights)
```
%% Cell type:code id:999756a9-4ccc-49df-920b-f387d8e51de8 tags:
``` python
plt.figure()
plt.subplot(131)
plt.title("mean")
plt.imshow(stitched_projection_mean)
plt.subplot(132)
plt.title("linear weight")
plt.imshow(stitched_projection_lin_weights)
plt.subplot(133)
plt.title("cosinus weight")
plt.imshow(stitched_projection_cos_weights)
```
%% Cell type:markdown id:bcdc8433-0018-47f7-b094-3cc5d0f7c9e8 tags:
## compute final stitched volume
%% Cell type:code id:71a02b15-1b02-4aa4-8a81-4f36c54134a2 tags:
``` python
# TODO
```
doc/tutorials/notebooks/stitching_preprocessing.ipynb
View file @ d98b1a34
......@@ -3,16 +3,20 @@
{
"cell_type": "markdown",
"id": "ed525300",
"metadata": {},
"metadata": {
"tags": []
},
"source": [
"# stitching in pre processing"
"# stitching projections (pre processing)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a6f2555f",
"metadata": {},
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"from tomoscan.serie import Serie\n",
......@@ -859,7 +863,7 @@
"metadata": {},
"outputs": [],
"source": [
"from tomwer.core.utils.stitching import stich_projections\n",
"from tomwer.core.utils.stitching import z_stitch_frames\n",
"proj_0_scan_1_shifted = shift_scipy(\n",
" proj_0_scan_1,\n",
" shift=(final_y_shift, final_x_shift),\n",
......@@ -867,10 +871,10 @@
")\n",
"if final_y_shift > 0:\n",
" proj_0_scan_1_shifted = proj_0_scan_1_shifted[final_y_shift:]\n",
"stitched_projection_mean, stitched_overlap_mean = stich_projections(\n",
"stitched_projection_mean, stitched_overlap_mean = z_stitch_frames(\n",
" proj_0=proj_0_scan_0,\n",
" proj_1=proj_0_scan_1_shifted,\n",
" y_overlap_in_px=y_overlap_in_pixel,\n",
" stitching_height_in_px=y_overlap_in_pixel,\n",
" strategy=\"mean\",\n",
")[0:2]"
]
......@@ -906,10 +910,10 @@
"metadata": {},
"outputs": [],
"source": [
"stitched_projection_lin_weights, stitched_overlap_lin_weights = stich_projections(\n",
"stitched_projection_lin_weights, stitched_overlap_lin_weights = z_stitch_frames(\n",
" proj_0=proj_0_scan_0,\n",
" proj_1=proj_0_scan_1_shifted,\n",
" y_overlap_in_px=y_overlap_in_pixel,\n",
" stitching_height_in_px=y_overlap_in_pixel,\n",
" strategy=\"linear weights\",\n",
")[0:2]"
]
......@@ -941,10 +945,10 @@
"metadata": {},
"outputs": [],
"source": [
"stitched_projection_cos_weights, stitched_overlap_cos_weights = stich_projections(\n",
"stitched_projection_cos_weights, stitched_overlap_cos_weights = z_stitch_frames(\n",
" proj_0=proj_0_scan_0,\n",
" proj_1=proj_0_scan_1_shifted,\n",
" y_overlap_in_px=y_overlap_in_pixel,\n",
" stitching_height_in_px=y_overlap_in_pixel,\n",
" strategy=\"cosinus weights\",\n",
")[0:2]"
]
......@@ -1029,16 +1033,16 @@
"outputs": [],
"source": [
"from nxtomomill.utils import Progress\n",
"from tomwer.core.utils.stitching import dump_stiched_projections\n",
"from tomwer.core.utils.stitching import dump_z_stiched_projections\n",
"from silx.io.url import DataUrl\n",
"\n",
"progress = Progress(\"\")\n",
"stitched_proj_url = dump_stiched_projections(\n",
"stitched_proj_url = dump_z_stiched_projections(\n",
" scan_0,\n",
" scan_1,\n",
" file_path=\"/tmp_14_days/payno/test_stitching/stitched_projections.hdf5\",\n",