Skip to content

Instantly share code, notes, and snippets.

@NicerNewerCar
Created August 24, 2023 14:42
Show Gist options
  • Select an option

  • Save NicerNewerCar/8aecc64913c3b6707002086d39b75da9 to your computer and use it in GitHub Desktop.

Select an option

Save NicerNewerCar/8aecc64913c3b6707002086d39b75da9 to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "a578193d-4ad4-42d7-9931-24cd43fdc469",
"metadata": {},
"outputs": [],
"source": [
"import itk\n",
"from itk import RTK as rtk\n",
"import numpy as np"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "f5b214ba-8e42-48ae-9640-8358c7911026",
"metadata": {},
"outputs": [],
"source": [
"volumeFName = \"myVolume.nrrd\"\n",
"\n",
"# Load CT Image\n",
"CT = itk.imread(volumeFName, pixel_type=itk.F)\n",
"\n",
"# Direction control\n",
"direction=np.array([[0.,1.,0.],[0.,0.,1.],[-1.,0.,0.]])\n",
"CT.SetDirection(itk.matrix_from_array(direction))\n",
"\n",
"# Set the origin to be around 0 - center of rotation\n",
"CT.SetOrigin([-0.5*(CT.GetLargestPossibleRegion().GetSize()[1]-1)*CT.GetSpacing()[1],\n",
" -0.5*(CT.GetLargestPossibleRegion().GetSize()[2]-1)*CT.GetSpacing()[2],\n",
" 0.5*(CT.GetLargestPossibleRegion().GetSize()[0]-1)*CT.GetSpacing()[0]])\n",
"\n",
"# Define the Image Type\n",
"ImageType = itk.Image[itk.F,3]"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "07b5b12e-05ad-4c1f-beef-12960ff7aee7",
"metadata": {},
"outputs": [],
"source": [
"# Empty Projection Stack\n",
"ConstantImageSourceType = rtk.ConstantImageSource[ImageType]\n",
"constantImageSource = ConstantImageSourceType.New()"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "3e4bc5c0-253a-4483-92f8-e3424cf1ccc1",
"metadata": {},
"outputs": [],
"source": [
"# Define the origin, size, and spacing of the projection stack \n",
"# Need to automatically get/calculate this info from the CT data\n",
"origin = [-200, -200, -500.]\n",
"nProjections = 360\n",
"size = [1000,1000,nProjections]\n",
"spacing = [0.39625,0.39625,0.625]\n",
"\n",
"constantImageSource.SetOrigin(origin)\n",
"constantImageSource.SetSpacing(spacing)\n",
"constantImageSource.SetSize(size)\n",
"constantImageSource.SetConstant(0.)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "9c76cd79-83d6-473e-a900-bdf3f060e478",
"metadata": {},
"outputs": [],
"source": [
"# RTK geometry\n",
"geo = rtk.ThreeDCircularProjectionGeometry.New()\n",
"firstAngle = 0.\n",
"angularArc = 360.\n",
"# Pull SID and SDD from Dicom Headers\n",
"sid = 541 # Source to Isocenter Distance\n",
"sdd = 949.075 # Source to Detector Distance\n",
"\n",
"for x in range(0,nProjections):\n",
" angle = firstAngle + x * angularArc / nProjections\n",
" geo.AddProjection(sid,sdd,angle)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "4d19f728-963d-48b1-85b6-1aa050ebfb81",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Writing the geometry to disk\n",
"xmlWriter = rtk.ThreeDCircularProjectionGeometryXMLFileWriter.New()\n",
"xmlWriter.SetFilename ( 'geo.xml' )\n",
"xmlWriter.SetObject ( geo )\n",
"xmlWriter.WriteFile()"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "d0a6c2c6-24cf-4af7-813e-d9b0bf31e960",
"metadata": {},
"outputs": [],
"source": [
"JFPType = rtk.JosephForwardProjectionImageFilter[ImageType, ImageType]\n",
"jfp = JFPType.New()\n",
"\n",
"jfp.SetGeometry(geo)\n",
"jfp.SetInput(0,constantImageSource.GetOutput())\n",
"jfp.SetInput(1,CT)\n",
"jfp.Update()"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "efc1c05b-5d8a-42f9-9e63-1e4aefb475c3",
"metadata": {},
"outputs": [],
"source": [
"itk.imwrite(jfp.GetOutput(), 'stack.tif')"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.4"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment