<jupyter_start><jupyter_text>MAPEM de Pierro algorithm for the Bowsher priorOne of the more popular methods for guiding a reconstruction based on a high quality image was suggested by Bowsher. This notebook explores this prior.We highly recommend you look at the [PET/MAPEM](../PET/MAPEM.ipynb) notebook first. This example extends upon the quadratic prior used in that notebook to use an anatomical prior. Authors: <NAME>, <NAME>, <NAME>, <NAME>  First version: 22nd of October 2019  Second version: 27th of October 2019  Third version: June 2021CCP SyneRBI Synergistic Image Reconstruction Framework (SIRF)  Copyright 2019,20201  University College London  Copyright 2019  King's College London  This is software developed for the Collaborative ComputationalProject in Synergistic Reconstruction for Biomedical Imaging. (http://www.synerbi.ac.uk/).Licensed under the Apache License, Version 2.0 (the "License");you may not use this file except in compliance with the License.You may obtain a copy of the License at    http://www.apache.org/licenses/LICENSE-2.0Unless required by applicable law or agreed to in writing, softwaredistributed under the License is distributed on an "AS IS" BASIS,WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.See the License for the specific language governing permissions andlimitations under the License.  Brief description of the Bowsher priorThe "usual" quadratic prior penalises differences between neighbouring voxels (using the square of the difference). This tends to oversmooth parts of the image where you know there should be an edge. To overcome this, it is natural to not penalise the difference between those "edge" voxels. This can be done after segmentation of the anatomical image for instance.Bowsher suggested a segmentation-free approach to use an anatomical (or any "side" image) as follows:- compute edge information on the anatomical image.- for each voxel, consider only the $N_B$ neighbours which have the lowest difference in the anatomical image.The paper is  Bowsher, J. E., <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, et al. ‘Utilizing MRI Information to Estimate F18-FDG Distributions in Rat Flank Tumors’. In IEEE Symposium Conference Record Nuclear Science 2004., 4:2488-2492 Vol. 4, 2004. https://doi.org/10.1109/NSSMIC.2004.1462760.  All the normal imports and handy functions<jupyter_code>%matplotlib notebook

# Setup the working directory for the notebook
import notebook_setup
from sirf_exercises import cd_to_working_dir
cd_to_working_dir('Synergistic', 'MAPEM_Bowsher')

#%% Initial imports etc
import numpy
import matplotlib.pyplot as plt
import os
import sys
import shutil
from tqdm.auto import tqdm, trange
import time
from scipy.ndimage.filters import gaussian_filter
import sirf.STIR as pet
from numba import jit
from sirf_exercises import exercises_data_path

brainweb_sim_data_path = exercises_data_path('working_folder', 'Synergistic', 'BrainWeb')
# set-up redirection of STIR messages to files
msg_red = pet.MessageRedirector('info.txt', 'warnings.txt', 'errors.txt')
# plotting settings
plt.ion() # interactive 'on' such that plots appear during loops

#%% some handy function definitions
def imshow(image, limits=None, title=''):
    """Usage: imshow(image, [min,max], title)"""
    plt.title(title)
    bitmap = plt.imshow(image)
    if limits is None:
        limits = [image.min(), image.max()]

    plt.clim(limits[0], limits[1])
    plt.colorbar(shrink=.6)
    plt.axis('off')
    return bitmap

def make_cylindrical_FOV(image):
    """truncate to cylindrical FOV"""
    filt = pet.TruncateToCylinderProcessor()
    filt.apply(image)

#%% define a function for plotting images and the updates
# This is the same function as in `ML_reconstruction`
def plot_progress(all_images, title, subiterations, cmax):
    if len(subiterations)==0:
        num_subiters = all_images[0].shape[0]-1
        subiterations = range(1, num_subiters+1)
    num_rows = len(all_images);
    slice_show = 60
    for it in subiterations:
        plt.figure()
        for r in range(num_rows):
            plt.subplot(num_rows,2,2*r+1)
            imshow(all_images[r][it,slice_show,:,:], [0,cmax], '%s at %d' % (title[r],  it))
            plt.subplot(num_rows,2,2*r+2)
            imshow(all_images[r][it,slice_show,:,:]-all_images[r][it-1,slice_show,:,:],[-cmax*.1,cmax*.1], 'update')
        plt.show();  

def subplot_(idx,vol,title,clims=None,cmap="viridis"):
    plt.subplot(*idx)
    plt.imshow(vol,cmap=cmap)
    if not clims is None:
        plt.clim(clims)
    plt.colorbar()
    plt.title(title)
    plt.axis("off")<jupyter_output><empty_output><jupyter_text>Load the dataTo generate the data needed for this notebook, run the [BrainWeb](./BrainWeb.ipynb) notebook first.<jupyter_code>full_acquired_data = pet.AcquisitionData(os.path.join(brainweb_sim_data_path, 'FDG_sino_noisy.hs'))
atten = pet.ImageData(os.path.join(brainweb_sim_data_path, 'uMap_small.hv'))

# Anatomical image
anatomical = pet.ImageData(os.path.join(brainweb_sim_data_path, 'T1_small.hv')) # could be T2_small.hv
anatomical_arr = anatomical.as_array()

# create initial image
init_image=atten.get_uniform_copy(atten.as_array().max()*.1)
make_cylindrical_FOV(init_image)

plt.figure()
imshow(anatomical.as_array()[64, :, :])
plt.show()
plt.figure()
imshow(full_acquired_data.as_array()[0, 64, :, :])
plt.show()<jupyter_output><empty_output><jupyter_text>Code from first MAPEM notebookThe following chunk of code is copied and pasted more-or-less directly from the other notebook as a starting point. First, run the code chunk to get the objective functions etc  construction of Likelihood objective functions and OSEM<jupyter_code>def get_obj_fun(acquired_data, atten):
    print('\n------------- Setting up objective function')
    #     #%% create objective function
    #%% create acquisition model
    am = pet.AcquisitionModelUsingRayTracingMatrix()
    am.set_num_tangential_LORs(5)

    # Set up sensitivity due to attenuation
    asm_attn = pet.AcquisitionSensitivityModel(atten, am)
    asm_attn.set_up(acquired_data)
    bin_eff = pet.AcquisitionData(acquired_data)
    bin_eff.fill(1.0)
    asm_attn.unnormalise(bin_eff)
    asm_attn = pet.AcquisitionSensitivityModel(bin_eff)

    # Set sensitivity of the model and set up
    am.set_acquisition_sensitivity(asm_attn)
    am.set_up(acquired_data,atten);

    #%% create objective function
    obj_fun = pet.make_Poisson_loglikelihood(acquired_data)
    obj_fun.set_acquisition_model(am)

    print('\n------------- Finished setting up objective function')
    return obj_fun

def get_reconstructor(num_subsets, num_subiters, obj_fun, init_image):
    print('\n------------- Setting up reconstructor') 

    #%% create OSEM reconstructor
    OSEM_reconstructor = pet.OSMAPOSLReconstructor()
    OSEM_reconstructor.set_objective_function(obj_fun)
    OSEM_reconstructor.set_num_subsets(num_subsets)
    OSEM_reconstructor.set_num_subiterations(num_subiters)

    #%% initialise
    OSEM_reconstructor.set_up(init_image)

    print('\n------------- Finished setting up reconstructor')
    return OSEM_reconstructor
# Use rebin to create a smaller sinogram to speed up calculations
acquired_data = full_acquired_data.clone()
acquired_data = acquired_data.rebin(3)
# Get the objective function
obj_fun = get_obj_fun(acquired_data, atten)<jupyter_output><empty_output><jupyter_text>Implement de Pierro MAP-EM for a quadratic prior with arbitrary weightsThe following code is almost a copy-paste of the implementation by <NAME> and <NAME> [contributed during one of our hackathons](https://github.com/SyneRBI/SIRF-Contribs/tree/master/src/Python/sirf/contrib/kcl). It is copied here for you to have an easier look.Note that the code avoids the `for` loops in our simplistic version above and hence should be faster (however, the construction of the neighbourhood is still slow, but you should have to do this only once). Also, this is a Python reimplementation of MATLAB code (hence the use of "Fortran order" below).<jupyter_code>def dePierroReg(image,weights,nhoodIndVec):
    """Get the de Pierro regularisation image (xreg)"""
    imSize = image.shape

    # vectorise image for indexing 
    imageVec = image.reshape(-1,order='F')

    # retrieve voxel intensities for neighbourhoods 
    resultVec = imageVec[nhoodIndVec]
    result = resultVec.reshape(weights.shape,order='F')

    # compute xreg
    imageReg = 0.5*numpy.sum(weights*(result + image.reshape(-1,1,order='F')),axis=1)/numpy.sum(weights,axis=1)
    imageReg = imageReg.reshape(imSize,order='F')

    return imageReg

def compute_nhoodIndVec(imageSize,weightsSize):
    """Get the neigbourhoods of each voxel"""
    w = int(round(weightsSize[1]**(1.0/3))) # side length of neighbourhood
    nhoodInd    = neighbourExtract(imageSize,w)
    return nhoodInd.reshape(-1,order='F')

def neighbourExtract(imageSize,w):
    """Adapted from kcl.Prior class"""
    n = imageSize[0]
    m = imageSize[1]
    h = imageSize[2]
    wlen = 2*numpy.floor(w/2)
    widx = xidx = yidx = numpy.arange(-wlen/2,wlen/2+1)

    if h==1:
        zidx = [0]
        nN = w*w
    else:
        zidx = widx
        nN = w*w*w

    Y,X,Z = numpy.meshgrid(numpy.arange(0,m), numpy.arange(0,n), numpy.arange(0,h))                
    N = numpy.zeros([n*m*h, nN],dtype='int32')
    l = 0
    for x in xidx:
        Xnew = setBoundary(X + x,n)
        for y in yidx:
            Ynew = setBoundary(Y + y,m)
            for z in zidx:
                Znew = setBoundary(Z + z,h)
                N[:,l] = ((Xnew + (Ynew)*n + (Znew)*n*m)).reshape(-1,1).flatten('F')
                l += 1
    return N

def setBoundary(X,n):
    """Boundary conditions for neighbourExtract.
    Adapted from kcl.Prior class"""
    idx = X<0
    X[idx] = X[idx] + n
    idx = X>n-1
    X[idx] = X[idx] - n
    return X.flatten('F')

@jit
def dePierroUpdate(xEM, imageReg, beta):
    """Update the image based on the de Pierro regularisation image"""
    return (2*xEM)/(((1 - beta*imageReg)**2 + 4*beta*xEM)**0.5 + (1 - beta*imageReg) + 0.00001)
def MAPEM_iteration(OSEM_reconstructor,current_image,weights,nhoodIndVec,beta):
    image_reg = dePierroReg(current_image.as_array(),weights,nhoodIndVec) # compute xreg
    OSEM_reconstructor.update(current_image); # compute EM update
    image_EM=current_image.as_array() # get xEM as a numpy array
    updated = dePierroUpdate(image_EM, image_reg, beta) # compute new uxpdate
    current_image.fill(updated) # store for next iteration
    return current_image<jupyter_output><empty_output><jupyter_text>Create uniform and Bowsher weightsWe will use the `kcl.Prior` class here to construct the Bowsher weights given an anatomical image. The `kcl.Prior` class (and the above code) assumes that the `weights` are returned an $N_v \times N_n$ array, with $N_v$ the number of voxels and $N_n$ the number of neighbours (here 27 as the implementation is in 3D).<jupyter_code>import sirf.contrib.kcl.Prior as pr
def update_bowsher_weights(prior,side_image,num_bowsher_neighbours):
    return prior.BowshserWeights\
        (side_image.as_array(),num_bowsher_neighbours)<jupyter_output><empty_output><jupyter_text>For illustration, we will keep only a few neighbours in the Bowsher prior. This makes the contrast with "uniform" weights higher of course.<jupyter_code>num_bowsher_neighbours = 3
myPrior = pr.Prior(anatomical_arr.shape)
BowsherWeights = update_bowsher_weights(myPrior,anatomical,num_bowsher_neighbours)<jupyter_output><empty_output><jupyter_text>Ignore the warning about `divide by zero`, it is actually handled in the `kcl.Prior` class.<jupyter_code># compute indices of the neighbourhood for each voxel
nhoodIndVec=compute_nhoodIndVec(anatomical_arr.shape,BowsherWeights.shape)
# illustrate that only a few of the weights in the neighbourhood are kept
# (taking an arbitrary voxel)
print(BowsherWeights[500,:])<jupyter_output><empty_output><jupyter_text>You could try to understand the neighbourhood structure using the following, but it is quite complicated due to the Fortran order and linear indices.```toLinearIndices=nhoodIndVec.reshape(BowsherWeights.shape,order='F')print(toLinearIndices[500,:])``` We will also use uniform weights where every neighbour is counted the same (often people will use 1/distance between voxels as weighting, but this isn't implemented here).<jupyter_code>uniformWeights=BowsherWeights.copy()
uniformWeights[:,:]=1
# set "self-weight" of the voxel to zero
uniformWeights[:,27//2]=0
print(uniformWeights[500,:])<jupyter_output><empty_output><jupyter_text>Run some experiments<jupyter_code>num_subsets = 21
num_subiters = 42<jupyter_output><empty_output><jupyter_text>Do a normal OSEM (for comparison and initialisation)<jupyter_code># Do initial OSEM recon
OSEM_reconstructor = get_reconstructor(num_subsets, num_subiters, obj_fun, init_image)
osem_image = init_image.clone()
OSEM_reconstructor.reconstruct(osem_image)

plt.figure()
imshow(osem_image.as_array()[60,:,:])
plt.show();<jupyter_output><empty_output><jupyter_text>Run MAP-EM with the 2 different sets of weightsTo save some time, we will initialise the algorithms with the OSEM image. This makes sense of course as in the initial iterations, the penalty will just slow everything down (as it smooths an already too smooth image even more!).<jupyter_code># arbitrary value for the weight of the penalty. You might have to tune it
beta=1<jupyter_output><empty_output><jupyter_text>Compute with Bowsher penalty<jupyter_code>current_image=osem_image.clone()

for it in trange(1, num_subiters+1):
    current_image = MAPEM_iteration(OSEM_reconstructor,current_image,BowsherWeights,nhoodIndVec,beta)
Bowsher=current_image.clone()<jupyter_output><empty_output><jupyter_text>Compute with uniform weights (we'll call the result UQP for "uniform quadratic penalty")<jupyter_code>current_image=osem_image.clone()

for it in trange(1, num_subiters+1):
    current_image = MAPEM_iteration(OSEM_reconstructor,current_image,uniformWeights,nhoodIndVec,beta)

UQP=current_image.clone()
# Plot the anatomical, OSEM, and two MAPEM images
plt.figure()
cmax=osem_image.max()*.6
clim=[0,cmax]
subplot_([1,2,1],anatomical.as_array()[60,:,:],"anatomical")
subplot_([1,2,2],osem_image.as_array()[60,:,:],"OSEM",clim)
plt.figure()
subplot_([1,2,1],UQP.as_array()[60,:,:],"Uniform Quadratic prior",clim)
subplot_([1,2,2],Bowsher.as_array()[60,:,:],"Bowsher Quadratic prior",clim)
plt.figure()
y_idx=osem_image.dimensions()[1]//2
plt.plot(osem_image.as_array()[60,y_idx,:],label="OSEM")
plt.plot(UQP.as_array()[60,y_idx,:],label="Uniform Quadratic prior")
plt.plot(Bowsher.as_array()[60,y_idx,:],label="Bowsher Quadratic prior")
plt.legend()<jupyter_output><empty_output><jupyter_text>You will probably see that the MAP-EM are quite smooth, and that there is very little difference between the "uniform" and "Bowsher" weights after this number of updates. The difference will get larger with higher number of updates (try it!).Also, with the Bowsher weights you should be able to increase `beta` more than for the uniform weights without oversmoothing the image too much.  Misalignment between anatomical and emission imagesWhat happens if you want to use an anatomical prior but the image isn't aligned with the image you're trying to reconstruct?  You'll have to register them of course! Have a look at the [registration notebook](../Reg/sirf_registration.ipynb) if you haven't already.  The idea here would be to run an initial reconstruction (say, OSEM), and then register the anatomical image to the resulting reconstruction...Once we've got the anatomical image in the correct space, we can calculate the Bowsher weights.<jupyter_code>import sirf.Reg as Reg

registration = Reg.NiftyAladinSym()
registration.set_reference_image
registration.set_reference_image(osem_image)
registration.set_floating_image(anatomical)
registration.set_parameter('SetPerformRigid','1')
registration.set_parameter('SetPerformAffine','0')
registration.process()
anatomical_in_emission_space = registration.get_output()

Bweights = update_bowsher_weights(myPrior,anatomical_in_emission_space,num_bowsher_neighbours)<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup