Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

adds an mlem notebook and 3.11 to CI (3.x is frequently broken) #184

Merged
merged 4 commits into from
Nov 22, 2023
Merged

adds an mlem notebook and 3.11 to CI (3.x is frequently broken) #184

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
5a0603b
Merge pull request #9 from FABLE-3DXRD/master
jonwright Nov 13, 2023
97fcedd
Thanks to Andrew Reader's youtube video
jonwright Nov 17, 2023
7b43cfe
Merge branch 'master' of https://github.com/jonwright/ImageD11
jonwright Nov 17, 2023
5204848
added 3.11 (3.x is busted due to numba)
jonwright Nov 22, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion .github/workflows/build_flake_pytest_ubuntu2004.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ jobs:
strategy:
fail-fast: false
matrix:
python-version: ["2.7", "3.8", "3.x"]
python-version: ["2.7", "3.8", "3.11", "3.x"]

steps:
- uses: actions/checkout@v3
Expand Down
281 changes: 281 additions & 0 deletions sandbox/mlem.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,281 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "d9f646df-ac47-4e8f-a65f-53922d9acaf1",
"metadata": {},
"source": [
"An \"MLEM\" algorithm from XRDUA was used in this paper:\n",
"\n",
"\"Impurity precipitation in atomized particles evidenced by nano x-ray diffraction computed tomography\"\n",
"Anne Bonnin; Jonathan P. Wright; Rémi Tucoulou; Hervé Palancher\n",
"Appl. Phys. Lett. 105, 084103 (2014)\n",
"https://doi.org/10.1063/1.4894009\n",
"\n",
"This python code implements something similar based on a youtube video (https://www.youtube.com/watch?v=IhETD4nSJec)\n",
"\n",
"There are lots of papers from mathematicians in the literature about MART (multiplicative ART). The conversion of latex algebra back and forth into computer code seems to be a bit of a problem for me (Jon Wright - 17 Nov 2023)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "cc820bc1-a9e9-4daa-b10f-fc23ef6a96d4",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"from ImageD11.sinograms.roi_iradon import iradon, radon # to have projection_shifts\n",
"from skimage.transform import iradon_sart\n",
"import numpy as np, pylab as pl\n",
"import skimage.data\n",
"\n",
"def backproject( sino, theta, projection_shifts = None ):\n",
" \"\"\" project the sinogram into the sample \"\"\"\n",
" return iradon( sino, theta, filter_name=None, projection_shifts = projection_shifts )\n",
"\n",
"def forwardproject( sample, theta, projection_shifts = None ):\n",
" \"\"\" project the sample into the experiment (sinogram) \"\"\"\n",
" return radon( sample, theta, projection_shifts = projection_shifts )\n",
"\n",
"def mlem( sino, theta, \n",
" startvalue = 1,\n",
" projection_shifts = None,\n",
" pad=0, niter=50, \n",
" divtol=1e-5, ):\n",
" # \n",
" # Also called \"MART\" for Multiplicative ART\n",
" # This keeps a positivity constraint for both the data and reconstruction\n",
" #\n",
" # This implementation was inspired from from:\n",
" # https://www.youtube.com/watch?v=IhETD4nSJec\n",
" # by Andrew Reader\n",
" #\n",
" # ToDo : implement a mask\n",
" # ToDo : check about the corners / circle=False aspects\n",
" #\n",
" # Number of pixels hitting each output in the sample:\n",
" sensitivity_image = backproject( np.ones_like(sino), theta,\n",
" projection_shifts = projection_shifts )\n",
" recip_sensitivity_image = 1./sensitivity_image\n",
" # The image reconstruction:\n",
" mlem_rec = np.empty( sensitivity_image.shape, np.float32)\n",
" mlem_rec[:] = startvalue\n",
" for i in range(niter):\n",
" calc_sino = forwardproject( mlem_rec, theta, projection_shifts = projection_shifts )\n",
" ratio = sino / (calc_sino + divtol )\n",
" correction = recip_sensitivity_image * backproject( ratio, theta, \n",
" projection_shifts = projection_shifts )\n",
" mlem_rec *= correction\n",
" return mlem_rec"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "78c6d5b3-a575-4286-a26a-da0cfbcbe588",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"im = skimage.data.shepp_logan_phantom().astype(np.float32)\n",
"pl.imshow(im)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "305bf8e6-64db-4c97-9d97-af7a85ebd639",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# generate some angles from a diffraction pattern\n",
"\n",
"import ImageD11.transform, ImageD11.unitcell\n",
"import scipy.spatial.transform\n",
"\n",
"# Angles to use are going to correspond to 5 rings of fcc\n",
"a = 3.06\n",
"cell = ImageD11.unitcell.unitcell( [a,a,a,90,90,90],'F')\n",
"cell.makerings(2)\n",
"hkls = []\n",
"for ds in cell.ringds[:5]:\n",
" hkls += cell.ringhkls[ds]\n",
"hkls = np.array( hkls )\n",
"hkls.shape"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "27e7634c-a59d-4134-9055-8eba0ef50c99",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"UB = scipy.spatial.transform.Rotation.random( random_state = 97 ).as_matrix()\n",
"gvecs = np.dot( UB/a, hkls.T )"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "38c231c6-bbfa-407c-9ae3-999c233883a4",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"tth, (eta0, eta1), (omega0, omega1) = ImageD11.transform.uncompute_g_vectors(gvecs, 12.39/50)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f745a86b-4a7b-4a22-8a74-e27a7e894678",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"measured = np.concatenate( [ omega0[ np.abs( np.sin( np.radians( eta0 ) ) ) > 0.1 ],\n",
" omega1[ np.abs( np.sin( np.radians( eta1 ) ) ) > 0.1 ] ] )\n",
"print(eta0.shape, measured.shape, measured.min(), measured.max())"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "dffc1087-70f5-4ca0-ad09-b67250dceaa1",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"order = np.argsort( measured )\n",
"theta = measured[ order ]\n",
"theta = theta[ theta > 0 ]"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3b17d387-bdc8-4633-9c48-46a1f700f435",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"sino = forwardproject( im, theta )"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6314be77-51a6-445a-814a-be1f305814de",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"pl.imshow(sino, aspect='auto', interpolation='nearest')"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "009c4844-6b8b-47a0-9ced-287e96fd4f54",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"def disp(im, title):\n",
" f, a = pl.subplots(1,2,figsize=(12,5))\n",
" f.colorbar( a[1].imshow(im), ax=a[0] )\n",
" a[1].set(title=title, xticks=(), yticks=())\n",
" a[0].plot( im[im.shape[0]//2] )\n",
" a[0].plot( im[:,im.shape[1]//2] )\n",
" pl.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "69a4c8d3-bb3c-4a4a-adcc-0ad024391234",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"fbp = iradon( sino, theta, output_size=im.shape[0])\n",
"disp(fbp,'Hann filter')"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8dd83a9a-41b5-44f8-a7a7-da920b6e649f",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"mlr = mlem( sino, theta, niter=1)\n",
"n = 1\n",
"disp( mlr, f'mlem {n} steps' )\n",
"for niter in (1,4,20,100):\n",
" mlr = mlem( sino, theta, startvalue = mlr, niter=niter)\n",
" n += niter\n",
" disp( mlr, f'mlem {n} steps' )"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fef2e35e-9330-475b-b8ff-56bcb2180247",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"rsart = skimage.transform.iradon_sart(sino, theta )\n",
"disp( rsart, 'iradon_sart')"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "62356c3e-9f7b-48a4-a0b4-6b29f01a75cf",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (main)",
"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.11.6"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Loading