<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Volumetric Registration and Analysis
In this tutorial we align a partial mouse brain to the Allen Referece Atlas (ARA).  Thus the ARA is our reference image, it consists of image data stored in the default "average" channel and corresponding annotations stored in a separate "annotation" channel.  We begin by importing the relevant modules
Step1:  We define the server and our user token
Step2:  First we'll download the atlas image
Step3:  Next we'll visuaize the image.  To ensure that the visuization is has good contrast we'll only show intensity values below the 99th percentile.
Step4:  When given an image volume imgShow displays it as a set of evenly spaced z-slices (1st column), y-slices (2nd column) and x slices (3rd column).
Step5:  Now we download the corresponding annotations
Step6:  It's usually better to visuaize a set of annotations in color.  Let's display the atlas annotations using a random colormap.  We begin by creating a 1000x3 array of random values. Since we always want the backgrond (label 0) to be dark we make the 0th row is [0,0,0].
Step7:  Now we can display the annotations.
Step8:  Let's overlay these annotations on the atlas image using alpha transparency.  To do this we must set newFig to False so that matplotlib appends the annotation data to the current figure instead of creating new one.
Step9:  Downloading input image
ndstore stores the brains at multiple resolutions. Resolution 0 is always the highest resolution each subsequent resoltion takes up 1/4th (under slice scaling) or 1/8th (under isotropic scaling) as many bytes as the previous resoltion.  We can use ndio's get_metadata method to see a list of available resoltions.
Step10:  Clearly resolution 1 is a resonable resolution.  So we'll download the image at that resolution.  Depending on your internet connection downloading may take several minutes.
Step11:  We'll only use intensity values below the 95th percentile in the visualization.
Step12:  Reorienting input image
You may have noticed that the input brain is not oriented in the same way as the atlas.  Let's look at the atlas.
Step13:  By examining the the z slices in the first column its clear that the x-axis goes from Right to Left side of the brain. The y-axis varies from Superior to Inferior and the z-axis goes from Anterior to posterior.  Thus it is in RSA orientation.  Looking at the the input brain...
Step14:  ...we see that the x-axis goes from Inferior to Superior, the y-axis goes from Anterior to Posterior and the Z axis goes from Left to Right.  Thus it's in IAL orientation.  Therefore we reorient the input image from IAL to RSA
Step15:  Compare the above to the Atlas.  The slices should now correspond with the atlas.  Let's make a copy of the reorianted image size and spacing because we'll need it later.
Step16:  Downample images
Now we downsample the input and reference images to a more manageable size
Step17:  Create and apply mask of reference image
Our input image only consists of part of the brain.  Therefore we mask the reference image to match the input image.  First we define a region of interest starting at the following point with the following size in mm.
Step18:  We then convert these values to from mm to voxels
Step19:  We create an Region Of Interest of value 255
Step20:  We then paste this into an empty image to create a mask of the reference image
Step21:  Now we apply this mask to our downsampled reference image
Step22:  Create input image mask
Notice how the input image has extremely bright regions.
Step23:  This will likly interfere with the registration.  Therefore we create a registration mask which excludes those regions
Step24:  Affine Registration
We can finally begin the registration. We will initialize the affine registration using the reference image's region of interest
Step25:  Now we compute the affine transform.  Unlike in the basic registration tutiorial we'll use the imgAffineComposite instead of imgAffine.  imgAffine simply computes the affine transform between the input and reference images.  imgAffineComposite computes a translation then a rigid then an affine transformation.  It's output is the composition of those three transforms.  We use Mutual Information since the input and reference images have very differnt intensity profiles.  We also enable the verbose option so that each iteration is printed.
Step26:  Now we apply the affine transform to the input image and mask
Step27:  We can evaluate the affine registration by generating a checkerboard of the reference and input images.  In this method the input image is placed on the black squares of the metaphorical checkerboard while the reference image is placed on the red squares.  Idealy anatomical sturctures should line up across squares.
Step28:  We also apply the  affine to the input mask
Step29:  LDDMM registration
Now we run LDDMM registration. Here we use imgMetamorphosisComposite.  Unlike imgMetamorphosis introduced in the 2D registration notebook, this function runs LDDMM in multiple steps using the alpha values specified by in alphaList. The field and invField outputs are the composition of all steps.  Once agan we use a Mutual Information cost because the input brain and reference image have very differnt intensity profiles.
Step30:  Now we apply the displacement field
Step31:  Evaluation LDDMM registration
Now we evaluate the deformable registration using a checkerboard image
Step32:  We can also evaluate the registration by overlaying the atlas annotations over the deformed input image.
Step33:  Uploading Results
Uploading deformed input image
We can now upload the atlas-aligned input brain back into ndstore.
Step34:  Uploading deformed atlas annotations
Step35:  Before we can overlay the atlas annotations on the downloaded image we must transfom them to its space. Fortunatly this can be done since all spatial transforms in this tutorial are invertable.  First we construct an inverse displacement field which transforms the annotations from the reference space to the input space before registration.
Step36:  Were not done yet.  We still need to reorient these annotations to their original
Step37:  We can upload these annotations at the lowest possible resolution.
<ASSISTANT_TASK:>
Python Code:
from ndreg import *
import matplotlib
import ndio.remote.neurodata as neurodata
Explanation: Volumetric Registration and Analysis
In this tutorial we align a partial mouse brain to the Allen Referece Atlas (ARA).  Thus the ARA is our reference image, it consists of image data stored in the default "average" channel and corresponding annotations stored in a separate "annotation" channel.  We begin by importing the relevant modules
End of explanation
server = "dev.neurodata.io"
userToken = txtRead("userToken.pem").strip()
Explanation: We define the server and our user token
End of explanation
refToken = "ara3"
refImg = imgDownload(refToken, channel="average", server=server, userToken=userToken)
Explanation: First we'll download the atlas image
End of explanation
refThreshold = imgPercentile(refImg, 0.99)
print(refThreshold)
Explanation: Next we'll visuaize the image.  To ensure that the visuization is has good contrast we'll only show intensity values below the 99th percentile.
End of explanation
imgShow(refImg, vmax=refThreshold)
Explanation: When given an image volume imgShow displays it as a set of evenly spaced z-slices (1st column), y-slices (2nd column) and x slices (3rd column).
End of explanation
refAnnoImg = imgDownload(refToken, channel="annotation", server=server, userToken=userToken)
refAnnoImgOrig = refAnnoImg[:,:,:]
imgShow(refAnnoImg, vmax=1000)
Explanation: Now we download the corresponding annotations
End of explanation
randValues = np.random.rand(1000,3)
randValues = np.concatenate(([[0,0,0]],randValues))
Explanation: It's usually better to visuaize a set of annotations in color.  Let's display the atlas annotations using a random colormap.  We begin by creating a 1000x3 array of random values. Since we always want the backgrond (label 0) to be dark we make the 0th row is [0,0,0].
End of explanation
randCmap = matplotlib.colors.ListedColormap (randValues)
imgShow(refAnnoImg, vmax=1000, cmap=randCmap)
Explanation: Now we can display the annotations.
End of explanation
imgShow(refImg, vmax=refThreshold, newFig=False)
imgShow(refAnnoImg, vmax=1000, cmap=randCmap, alpha=0.3, newFig=False)
plt.show()
Explanation: Let's overlay these annotations on the atlas image using alpha transparency.  To do this we must set newFig to False so that matplotlib appends the annotation data to the current figure instead of creating new one.
End of explanation
inToken = "s275"
nd = neurodata(hostname="dev.neurodata.io", user_token=userToken)
print(nd.get_metadata(inToken)['dataset']['imagesize']) 
Explanation: Downloading input image
ndstore stores the brains at multiple resolutions. Resolution 0 is always the highest resolution each subsequent resoltion takes up 1/4th (under slice scaling) or 1/8th (under isotropic scaling) as many bytes as the previous resoltion.  We can use ndio's get_metadata method to see a list of available resoltions.
End of explanation
inImg = imgDownload(inToken, resolution=1, userToken=userToken, server=server)
inImg.SetSpacing(np.array(inImg.GetSpacing())*1000) ###
inImgOrig = inImg[:,:,:]
Explanation: Clearly resolution 1 is a resonable resolution.  So we'll download the image at that resolution.  Depending on your internet connection downloading may take several minutes.
End of explanation
inImg = inImgOrig[:,:,:]
inThreshold = imgPercentile(inImg, 0.95)
imgShow(inImg, vmax=inThreshold)
Explanation: We'll only use intensity values below the 95th percentile in the visualization.
End of explanation
imgShow(refImg, vmax=refThreshold)
Explanation: Reorienting input image
You may have noticed that the input brain is not oriented in the same way as the atlas.  Let's look at the atlas.
End of explanation
imgShow(inImg, vmax=inThreshold)
Explanation: By examining the the z slices in the first column its clear that the x-axis goes from Right to Left side of the brain. The y-axis varies from Superior to Inferior and the z-axis goes from Anterior to posterior.  Thus it is in RSA orientation.  Looking at the the input brain...
End of explanation
inOrient = "IAL"
refOrient = "RSA"
inImg = imgReorient(inImg, inOrient, refOrient)
imgShow(inImg, vmax=inThreshold)
Explanation: ...we see that the x-axis goes from Inferior to Superior, the y-axis goes from Anterior to Posterior and the Z axis goes from Left to Right.  Thus it's in IAL orientation.  Therefore we reorient the input image from IAL to RSA
End of explanation
inImgSize_reorient = inImg.GetSize()
inImgSpacing_reorient= inImg.GetSpacing()
Explanation: Compare the above to the Atlas.  The slices should now correspond with the atlas.  Let's make a copy of the reorianted image size and spacing because we'll need it later.
End of explanation
spacing = [0.1,0.1, 0.1]
inImg_ds = imgResample(inImg, spacing)
imgShow(inImg_ds, vmax=inThreshold)
refImg_ds = imgResample(refImg, spacing)
imgShow(refImg_ds, vmax=refThreshold)
imgWrite(inImg_ds,"/cis/project/clarity/data/ailey/s275_ch0_rsa_100um.img")
Explanation: Downample images
Now we downsample the input and reference images to a more manageable size
End of explanation
roiStart = [5.4, 1.2, 2.1]
roiSize = [4.5,6.5,7.5]
Explanation: Create and apply mask of reference image
Our input image only consists of part of the brain.  Therefore we mask the reference image to match the input image.  First we define a region of interest starting at the following point with the following size in mm.
End of explanation
roiStartVoxel = (roiStart / np.array(spacing)).astype('uint16').tolist()
print(roiStartVoxel)
roiSizeVoxel = (roiSize / np.array(spacing)).astype('uint16').tolist()
print(roiSizeVoxel)
Explanation: We then convert these values to from mm to voxels
End of explanation
roiImg = sitk.Image(roiSizeVoxel,sitk.sitkUInt8)
roiImg += 255
Explanation: We create an Region Of Interest of value 255
End of explanation
emptyImg = sitk.Image(refImg_ds.GetSize(),sitk.sitkUInt8) # Create an empty image
emptyImg.CopyInformation(refImg_ds) # Copy spacing, origin and direction from reference image
refMask = sitk.Paste(emptyImg, roiImg, roiSizeVoxel, [0,0,0], roiStartVoxel)
imgShow(refMask, vmax=255)
Explanation: We then paste this into an empty image to create a mask of the reference image
End of explanation
refImg_ds = sitk.Mask(refImg_ds, refMask)
imgShow(refImg_ds, vmax=refThreshold)
Explanation: Now we apply this mask to our downsampled reference image
End of explanation
imgShow(inImg_ds, vmax=inThreshold)
Explanation: Create input image mask
Notice how the input image has extremely bright regions.
End of explanation
threshold = imgPercentile(inImg_ds,0.95)
inMask_ds = sitk.BinaryThreshold(inImg_ds, 0, threshold, 255, 0)
imgShow(inMask_ds, vmax=255)
Explanation: This will likly interfere with the registration.  Therefore we create a registration mask which excludes those regions
End of explanation
translation = -np.array(roiStart)
inAffine = [1.2,0,0,0,1.2,0,0,0,1]+translation.tolist()
print(inAffine)
imgShow(imgApplyAffine(inImg_ds, inAffine, size=refImg_ds.GetSize()),vmax = inThreshold)
Explanation: Affine Registration
We can finally begin the registration. We will initialize the affine registration using the reference image's region of interest
End of explanation
affine = imgAffineComposite(inImg_ds, refImg_ds, inMask=inMask_ds, iterations=100, useMI=True, verbose=True, inAffine=inAffine)
Explanation: Now we compute the affine transform.  Unlike in the basic registration tutiorial we'll use the imgAffineComposite instead of imgAffine.  imgAffine simply computes the affine transform between the input and reference images.  imgAffineComposite computes a translation then a rigid then an affine transformation.  It's output is the composition of those three transforms.  We use Mutual Information since the input and reference images have very differnt intensity profiles.  We also enable the verbose option so that each iteration is printed.
End of explanation
inImg_affine = imgApplyAffine(inImg, affine, size=refImg.GetSize(), spacing=refImg.GetSpacing())
imgShow(inImg_affine, vmax=inThreshold)
Explanation: Now we apply the affine transform to the input image and mask
End of explanation
inImg_ds = imgResample(inImg_affine, spacing=spacing, size=refImg_ds.GetSize())
imgShow(imgChecker(inImg_ds, refImg_ds), vmax=refThreshold)
Explanation: We can evaluate the affine registration by generating a checkerboard of the reference and input images.  In this method the input image is placed on the black squares of the metaphorical checkerboard while the reference image is placed on the red squares.  Idealy anatomical sturctures should line up across squares.
End of explanation
inMask_ds = imgApplyAffine(inMask_ds, affine, useNearest=True, size=refImg_ds.GetSize())
imgShow(inMask_ds, vmax=255)
imgShow(inImg_ds, vmax=inThreshold)
Explanation: We also apply the  affine to the input mask
End of explanation
inImg_ds = imgResample(inImg_affine, spacing=spacing, size=refImg_ds.GetSize())
(field, invField) = imgMetamorphosisComposite(inImg_ds, refImg_ds, inMask=inMask_ds, alphaList=[0.1, 0.05,0.02],
                                              scaleList = [1.0, 1.0,1.0], useMI=True, iterations=100, verbose=True)
Explanation: LDDMM registration
Now we run LDDMM registration. Here we use imgMetamorphosisComposite.  Unlike imgMetamorphosis introduced in the 2D registration notebook, this function runs LDDMM in multiple steps using the alpha values specified by in alphaList. The field and invField outputs are the composition of all steps.  Once agan we use a Mutual Information cost because the input brain and reference image have very differnt intensity profiles.
End of explanation
inImg_lddmm = imgApplyField(inImg_affine, field, size=refImg.GetSize())
imgShow(inImg_lddmm, vmax=inThreshold)
Explanation: Now we apply the displacement field
End of explanation
inImg_ds = imgResample(inImg_lddmm, spacing=spacing, size=refImg_ds.GetSize())
imgShow(imgChecker(inImg_ds, refImg_ds), vmax=refThreshold)
Explanation: Evaluation LDDMM registration
Now we evaluate the deformable registration using a checkerboard image
End of explanation
imgShow(inImg_lddmm, vmax=inThreshold, newFig=False)
imgShow(refAnnoImg, vmax=1000, cmap=randCmap, alpha=0.2, newFig=False)
plt.show()
Explanation: We can also evaluate the registration by overlaying the atlas annotations over the deformed input image.
End of explanation
outToken = inToken + "_to_" + refToken
imgUpload(inImg_lddmm, outToken, server=server, userToken=userToken)
### imgWrite(inImg_lddmm, "/cis/project/clarity/data/ailey/"+outToken+"_new.img")
Explanation: Uploading Results
Uploading deformed input image
We can now upload the atlas-aligned input brain back into ndstore.
End of explanation
spacing_ds = invField.GetSpacing()
size_ds = np.ceil(np.array(refAnnoImg.GetSize())*np.array(refAnnoImg.GetSpacing())/np.array(spacing_ds))
size_ds = list(size_ds.astype(int))
Explanation: Uploading deformed atlas annotations
End of explanation
invAffine = affineInverse(affine)
invAffineField = affineToField(invAffine, size_ds, spacing_ds)
invField2 = fieldApplyField(invAffineField, invField)
inAnnoImg = imgApplyField(refAnnoImg, invField2,useNearest=True, size=inImgSize_reorient, spacing=inImgSpacing_reorient)
inAnnoThreshold = imgPercentile(inAnnoImg,0.99)
imgShow(inAnnoImg, vmax=inAnnoThreshold)
Explanation: Before we can overlay the atlas annotations on the downloaded image we must transfom them to its space. Fortunatly this can be done since all spatial transforms in this tutorial are invertable.  First we construct an inverse displacement field which transforms the annotations from the reference space to the input space before registration.
End of explanation
inAnnoImg = imgReorient(inAnnoImg, refOrient, inOrient)
imgShow(inAnnoImg, vmax=inAnnoThreshold)
Explanation: Were not done yet.  We still need to reorient these annotations to their original
End of explanation
outToken = "ara3_to_AutA"
outChannel = "annotation_draft"
imgUpload(inAnnoImg, outToken, outChannel, resolution=5)
Explanation: We can upload these annotations at the lowest possible resolution.
End of explanation
<END_TASK>