diff --git a/notebooks/SPECT/SPECT_OSEM.ipynb b/notebooks/SPECT/SPECT_OSEM.ipynb index 1306f1f3..902eef02 100644 --- a/notebooks/SPECT/SPECT_OSEM.ipynb +++ b/notebooks/SPECT/SPECT_OSEM.ipynb @@ -8,6 +8,21 @@ "A notebook to demonstrate the setup and basic OSEM reconstruction of a 2-dimensional dummy image using SIRF's SPECT projector" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Authors: Daniel Deidda and Sam Porter \n", + "\n", + "CCP SyneRBI Synergistic Image Reconstruction Framework (SIRF). \n", + "Copyright 2022 - 2024 National Physical Laboratory \n", + "Copyright 2022 - 2024, 2021 University College London\n", + "\n", + "This is software developed for the Collaborative Computational Project in Synergistic Reconstruction for Biomedical Imaging (http://www.ccpsynerbi.ac.uk/).\n", + "\n", + "SPDX-License-Identifier: Apache-2.0" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -183,7 +198,8 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Now we need to create the acquisition model by first creating an acquisition matrix (SPECTUBMatrix) object and apply this to a sirf Acquisition Model (AcquisitionModelUsingMatrix) object" + "Now we need to create the acquisition model by first creating an acquisition matrix (SPECTUBMatrix) object, add attenuation, PSF, and apply this to a sirf Acquisition Model \n", + "(AcquisitionModelUsingMatrix) object. Hint (help(spect.SPECTUBMatrix))" ] }, { @@ -192,7 +208,9 @@ "metadata": {}, "outputs": [], "source": [ - "### Acquisition Model code here ###" + "### Acquisition Model code here ###\n", + "# acq_model_matrix =?\n", + "help(spect.SPECTUBMatrix)" ] }, { @@ -239,12 +257,21 @@ "metadata": {}, "outputs": [], "source": [ - "### Objective Function code here ###\n", + "# create objective function\n", + "obj_fun = spect.make_Poisson_loglikelihood(noisy_data)\n", "\n", + "### and now we apply a different resolution mdoel ### and\n", + "acq_model_matrix.set_resolution_model(0.1,0.1,full_3D=False)\n", + "acq_model_2 = spect.AcquisitionModelUsingMatrix(acq_model_matrix)\n", + "obj_fun.set_acquisition_model(acq_model_2)\n", + "\n", + "# create OSEM reconstructor object\n", "num_subsets = 21 # number of subsets for OSEM reconstruction\n", "num_subiters = 42 #number of subiterations (i.e two full iterations)\n", - "\n", - "### OSEM reconstructor code here ###" + "OSEM_reconstructor = spect.OSMAPOSLReconstructor()\n", + "OSEM_reconstructor.set_objective_function(obj_fun)\n", + "OSEM_reconstructor.set_num_subsets(num_subsets)\n", + "OSEM_reconstructor.set_num_subiterations(num_subiters)" ] }, { @@ -326,12 +353,15 @@ }, "language_info": { "codemirror_mode": { - "name": "ipython" + "name": "ipython", + "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", - "nbconvert_exporter": "python" + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.12" } }, "nbformat": 4, diff --git a/notebooks/SPECT/SPECT_OSEM_measured_data.ipynb b/notebooks/SPECT/SPECT_OSEM_measured_data.ipynb index dda7d53c..2cf19833 100644 --- a/notebooks/SPECT/SPECT_OSEM_measured_data.ipynb +++ b/notebooks/SPECT/SPECT_OSEM_measured_data.ipynb @@ -165,12 +165,15 @@ }, "language_info": { "codemirror_mode": { - "name": "ipython" + "name": "ipython", + "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", - "nbconvert_exporter": "python" + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.12" } }, "nbformat": 4, diff --git a/notebooks/SPECT/SPECT_OSEM_measured_data_solution.ipynb b/notebooks/SPECT/SPECT_OSEM_measured_data_solution.ipynb index 0297f80d..2905ce79 100644 --- a/notebooks/SPECT/SPECT_OSEM_measured_data_solution.ipynb +++ b/notebooks/SPECT/SPECT_OSEM_measured_data_solution.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "962178b0", + "id": "0", "metadata": {}, "source": [ "Simple OSEM reconstruction demo for real data: \n", @@ -28,7 +28,7 @@ }, { "cell_type": "markdown", - "id": "d91eaa49", + "id": "1", "metadata": {}, "source": [ "In this exercise you are going to apply what you learned in the previous notebooks about SPECT reconstruction to reconstruct real data. \n" @@ -37,7 +37,7 @@ { "cell_type": "code", "execution_count": null, - "id": "f80fd2f7", + "id": "2", "metadata": {}, "outputs": [], "source": [ @@ -63,7 +63,7 @@ { "cell_type": "code", "execution_count": null, - "id": "9e8458d6", + "id": "3", "metadata": {}, "outputs": [], "source": [ @@ -75,7 +75,7 @@ }, { "cell_type": "markdown", - "id": "08524093", + "id": "4", "metadata": {}, "source": [ "The following is not needed if the data is already downloaded" @@ -84,7 +84,7 @@ { "cell_type": "code", "execution_count": null, - "id": "bbca69aa", + "id": "5", "metadata": {}, "outputs": [], "source": [ @@ -97,7 +97,7 @@ { "cell_type": "code", "execution_count": null, - "id": "931f651d", + "id": "6", "metadata": {}, "outputs": [], "source": [ @@ -110,7 +110,7 @@ { "cell_type": "code", "execution_count": null, - "id": "6ec9dcbc", + "id": "7", "metadata": {}, "outputs": [], "source": [ @@ -120,7 +120,7 @@ }, { "cell_type": "markdown", - "id": "54ce83d9", + "id": "8", "metadata": {}, "source": [ "# Exercise 1: Simple Reconstruction:\n", @@ -130,7 +130,7 @@ { "cell_type": "code", "execution_count": null, - "id": "58584f1f", + "id": "9", "metadata": {}, "outputs": [], "source": [ @@ -145,7 +145,7 @@ { "cell_type": "code", "execution_count": null, - "id": "850fd91f", + "id": "10", "metadata": {}, "outputs": [], "source": [ @@ -158,7 +158,7 @@ { "cell_type": "code", "execution_count": null, - "id": "d14521ed", + "id": "11", "metadata": {}, "outputs": [], "source": [ @@ -173,7 +173,7 @@ { "cell_type": "code", "execution_count": null, - "id": "4479d023", + "id": "12", "metadata": {}, "outputs": [], "source": [ @@ -188,7 +188,7 @@ { "cell_type": "code", "execution_count": null, - "id": "0663ab88", + "id": "13", "metadata": {}, "outputs": [], "source": [ @@ -200,7 +200,7 @@ { "cell_type": "code", "execution_count": null, - "id": "82b2fcec", + "id": "14", "metadata": {}, "outputs": [], "source": [ @@ -216,7 +216,7 @@ { "cell_type": "code", "execution_count": null, - "id": "dd1042e3", + "id": "15", "metadata": {}, "outputs": [], "source": [ @@ -231,7 +231,7 @@ { "cell_type": "code", "execution_count": null, - "id": "51264d96", + "id": "16", "metadata": {}, "outputs": [], "source": [ @@ -244,7 +244,7 @@ }, { "cell_type": "markdown", - "id": "57e0ffc1", + "id": "17", "metadata": {}, "source": [ "# Exercise 2: PSF Reconstruction:\n", @@ -260,7 +260,7 @@ { "cell_type": "code", "execution_count": null, - "id": "3f5fce65", + "id": "18", "metadata": {}, "outputs": [], "source": [ @@ -282,7 +282,7 @@ { "cell_type": "code", "execution_count": null, - "id": "0c56f23b", + "id": "19", "metadata": {}, "outputs": [], "source": [ diff --git a/notebooks/SPECT/SPECT_OSEM_solution.ipynb b/notebooks/SPECT/SPECT_OSEM_solution.ipynb index 84df3238..f145336a 100644 --- a/notebooks/SPECT/SPECT_OSEM_solution.ipynb +++ b/notebooks/SPECT/SPECT_OSEM_solution.ipynb @@ -8,6 +8,21 @@ "A notebook to demonstrate the setup and basic OSEM reconstruction of a 2-dimensional dummy image using SIRF's SPECT projector" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Authors: Daniel Deidda and Sam Porter \n", + "\n", + "CCP SyneRBI Synergistic Image Reconstruction Framework (SIRF). \n", + "Copyright 2022 - 2024 National Physical Laboratory \n", + "Copyright 2022 - 2024, 2021 University College London\n", + "\n", + "This is software developed for the Collaborative Computational Project in Synergistic Reconstruction for Biomedical Imaging (http://www.ccpsynerbi.ac.uk/).\n", + "\n", + "SPDX-License-Identifier: Apache-2.0" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -194,16 +209,20 @@ "metadata": {}, "outputs": [], "source": [ + "help(spect.SPECTUBMatrix.set_resolution_model)\n", "# select acquisition model that implements the geometric\n", - "# forward projection by a ray tracing matrix multiplication\n", + "\n", "acq_model_matrix = spect.SPECTUBMatrix();\n", + "def get_acquisition_model(uMap,templ_sino,resol_slope, resol_sigma0):\n", + " # forward projection by a ray tracing matrix multiplication\n", "\n", - "### here we add attenuation and set a PSF for the aquisition model ###\n", - "acq_model_matrix.set_attenuation_image(uMap) # add attenuation\n", - "acq_model_matrix.set_resolution_model(0.2,0.2,full_3D=False) #resolution modelling\n", + " ### here we add attenuation and set a PSF for the aquisition model ###\n", + " acq_model_matrix.set_attenuation_image(uMap) # add attenuation\n", + " acq_model_matrix.set_resolution_model(resol_slope,resol_sigma0,full_3D=False) #resolution modelling\n", "\n", - "acq_model = spect.AcquisitionModelUsingMatrix(acq_model_matrix)\n", - "acq_model.set_up(templ_sino, image)" + " acq_model = spect.AcquisitionModelUsingMatrix(acq_model_matrix)\n", + " acq_model.set_up(templ_sino, uMap)\n", + " return acq_model" ] }, { @@ -222,7 +241,7 @@ "print('projecting image...')\n", "# project the image to obtain simulated acquisition data\n", "# data from raw_data_file is used as a template\n", - "acq_model.set_up(templ_sino, image)\n", + "acq_model= get_acquisition_model(uMap,templ_sino,0.1,0.1)\n", "simulated_data = templ_sino.get_uniform_copy()\n", "acq_model.forward(image, 0, 1, simulated_data)\n", "\n", @@ -271,6 +290,13 @@ "OSEM_reconstructor.set_up(init_image)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And now for the reconstruction..." + ] + }, { "cell_type": "code", "execution_count": null, @@ -288,7 +314,103 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "What is different about this image?" + "How did we do?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Exercise: reconstruct with and without attenuation and resolution modelling\n", + "* Investigate the effect of resolution modelling and attenuation on the reconstructed image\n", + "* What happens if you have a different resolution model for the simulated data and the reconstruction?\n", + "\n", + "Hint: use the code created before for the acquisition model, a zero umap as the effect of no attenuation correction" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "acq_model_noatten = get_acquisition_model(uMap.clone().fill(np.zeros(uMap_array.shape)),templ_sino,0.1, 0.1)\n", + "\n", + "# create initialisation image and set up reconstructor\n", + "obj_fun = spect.make_Poisson_loglikelihood(noisy_data)\n", + "obj_fun.set_acquisition_model(acq_model_noatten)\n", + "\n", + "# create OSEM reconstructor object\n", + "num_subsets = 21 # number of subsets for OSEM reconstructionc\n", + "num_subiters = 42 #number of subiterations (i.e two full iterations)\n", + "OSEM_reconstructor = spect.OSMAPOSLReconstructor()\n", + "OSEM_reconstructor.set_objective_function(obj_fun)\n", + "OSEM_reconstructor.set_num_subsets(num_subsets)\n", + "OSEM_reconstructor.set_num_subiterations(num_subiters)\n", + "init_image = make_cylindrical_FOV(image.get_uniform_copy(1))\n", + "OSEM_reconstructor.set_up(init_image)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "OSEM_reconstructor.reconstruct(init_image)\n", + "out_image = OSEM_reconstructor.get_current_estimate()\n", + "out_image_array = out_image.as_array()\n", + "show_2D_array('Reconstructed image', out_image_array[0,:,:])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "now change let's change the sigma of the PSF and see what happens" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "acq_model_psf = get_acquisition_model(uMap,templ_sino,0.1, 0.001)\n", + "\n", + "# create initialisation image and set up reconstructor\n", + "obj_fun = spect.make_Poisson_loglikelihood(noisy_data)\n", + "obj_fun.set_acquisition_model(acq_model_psf)\n", + "\n", + "# create OSEM reconstructor object\n", + "num_subsets = 21 # number of subsets for OSEM reconstructionc\n", + "num_subiters = 42 #number of subiterations (i.e two full iterations)\n", + "OSEM_reconstructor = spect.OSMAPOSLReconstructor()\n", + "OSEM_reconstructor.set_objective_function(obj_fun)\n", + "OSEM_reconstructor.set_num_subsets(num_subsets)\n", + "OSEM_reconstructor.set_num_subiterations(num_subiters)\n", + "init_image = make_cylindrical_FOV(image.get_uniform_copy(1))\n", + "OSEM_reconstructor.set_up(init_image)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# OSEM_reconstructor.reconstruct(init_image)\n", + "# out_image = OSEM_reconstructor.get_current_estimate()\n", + "out_image_array = out_image.as_array()\n", + "show_2D_array('Reconstructed image', out_image_array[0,:,:])\n", + "show_2D_array('Reconstructed image', uMap_array[0,:,:])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "What is different about this images?" ] }, { @@ -315,12 +437,15 @@ }, "language_info": { "codemirror_mode": { - "name": "ipython" + "name": "ipython", + "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", - "nbconvert_exporter": "python" + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.12" } }, "nbformat": 4,