Skip to content

Commit

Permalink
Merge pull request #184 from jonwright/master
Browse files Browse the repository at this point in the history 
adds an mlem notebook and 3.11 to CI (3.x is frequently broken)
  • Loading branch information
@jonwright
jonwright authored Nov 22, 2023
2 parents 22e7759 + 5204848 commit 17ab635
Show file tree
Hide file tree
Showing 2 changed files with 282 additions and 1 deletion.
2 changes: 1 addition & 1 deletion .github/workflows/build_flake_pytest_ubuntu2004.yml
View file
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
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
}

0 comments on commit 17ab635

Please sign in to comment.