Custom Image Phantom & Sinograms

[singhtejus]

Hey! I was working on implementing the sinograms.jl package to create a QRM phantom. Since the original examples were given using the Shepp Logan phantom, some parameters had to be tweaked to make the CT Fan Scanner work for the new phantom. Below are the results and code.

The sinogram had many crossing lines through it, which was fixed by increasing the number of angular samples. This however introduced blurriness to the image.

I have a few questions on this package:

• How can I make this sinogram better? Are there parameters that would work better for the phantom I have?
• The Shepp Logan phantom was hard-coded with specific materials, but when building my own phantom, how can I replicate the material (ie change the density) of the different parts of the phantom to simulate lungs, calcium (of various densities), etc?
• The CT energy seems to be hardcoded into the package. Is there any way to change the tube voltages to have a variable energy CT scanner?

# ╔═╡ abf12001-5afe-4ff6-89c9-74aae9d38ff8
begin
	using ImagePhantoms
	using MIRTjim: jim
	using Unitful: mm,cm, unit, °, ustrip
	using ImageGeoms: ImageGeom, MaskCircle
	using CairoMakie
	using PlutoUI
	using Sinograms

end

# ╔═╡ 8aebcddd-2e52-4b03-a66e-8b145a9bc714
TableOfContents()

# ╔═╡ 387f7c9c-cf68-4ae5-9333-42434596a131
md"""
# QRM Phantom
"""

# ╔═╡ 680f7e0e-4f55-48e8-8630-9095f516144d
md"""
## Define Image Geometry & Phantom
"""

# ╔═╡ 5e386415-2234-46f9-9278-2e311f6ed391
begin
	#tissue
	center = (0mm, 0mm, 0mm)
	width = (150mm, 100mm, 15mm) # x radius, y radius, height

	#left lung
	center2 = (-82.5mm, 0mm, center[3])
	width2 = (27.5mm, 60mm, 15mm)
	
	#right lung
	center3 = (82.5mm, 0mm, center[3])
	width3 = (27.5mm, 60mm, 15mm)

	#center circle
	center4 = (0mm, 0mm, 0mm)
	width4 = (40mm, 40mm, 15mm)

	#1mm diameter calcium inserts
	center5 = (5mm,5mm,0mm)
	center5b = (-5mm,5mm,0mm)
	center5c = (0mm,-5mm,0mm)
	width5 = (1mm,1mm,15mm)

	#3mm diameter calcium inserts
	center6 = (10mm,10mm,0mm)
	center6b = (-10mm,10mm,0mm)
	center6c = (0mm,-10mm,0mm)
	width6 = (1.5mm,1.5mm,15mm)
	
	#5mm diameter calcium inserts
	center7 = (15mm,15mm,0mm)
	center7b = (-15mm,15mm,0mm)
	center7c = (0mm,-15mm,0mm)
	width7 = (2.5mm,2.5mm,15mm)

	#spine insert 
	center8 = (0mm, -57mm,0mm)
	width8 = (15mm, 15mm, 15mm)
	center9 = (0mm, -85mm,0mm)
	width9 = (25mm, 8mm, 15mm)
		
end

# ╔═╡ d73b27b1-5199-4b27-8465-c3bab0e5af26
begin
	angles = (0, 0, 0)
	angles2 = (0, 0, 0)
	angles3 = (0, 0, 0)
	angles4 = (0, 0, 0)
	angles5 = (0, 0, 0)
	angles6 = (3π/2,0,0)
end

# ╔═╡ 0a0a9a61-b993-447a-8b54-1a561e22b58b
ob1 = cylinder(center, width, angles, 10); # tissue

# ╔═╡ e4f97cda-4836-4e90-95a0-a750983f4e16
ob2 = cylinder(center2, width2, angles2, -500); # lung

# ╔═╡ 84c6b7ac-80d5-4737-b20a-b845b869ebe1
ob3 = cylinder(center3, width3, angles3, -500); # lung

# ╔═╡ 3b9286c0-6803-4bea-b216-2065851802a1
ob4 = cylinder(center4, width4, angles4, 40); # myocardium

# ╔═╡ 32893dc8-15dd-4311-9e6b-de12f57f3260
begin #1mm diameter calcium inserts
	ob5 = cylinder(center5, width5, angles5, -200)
	ob5b = cylinder(center5b, width5, angles5, -200)
	ob5c = cylinder(center5c, width5, angles5, -200)
end;

# ╔═╡ 572e9c2d-1631-4c0e-a563-a40c1253e267
begin #3mm diameter calcium inserts
	ob6 = cylinder(center6, width6, angles5, -200)
	ob6b = cylinder(center6b, width6, angles5, -200)
	ob6c = cylinder(center6c, width6, angles5, -200)
end;

# ╔═╡ 37ba8f65-9c40-4809-81e3-88ae4c5e7def
begin #5mm diameter calcium inserts
	ob7 = cylinder(center7, width7, angles5, -200)
	ob7b = cylinder(center7b, width7, angles5, -200)
	ob7c = cylinder(center7c, width7, angles5, -200)
end;

# ╔═╡ 33038548-64b0-4a6f-bb92-2e803251b9ab
begin #spine insert
	ob8 = cylinder(center8,width8,angles5,40)
	ob9 = cuboid(center9,width9,angles6,40)
end;

# ╔═╡ 6432b8e6-41f1-4e2e-b584-2e13178428cc
objects = [ob1,ob2,ob3,ob4,ob5,ob5b,ob5c,ob6,ob6b,ob6c,ob7,ob7b,ob7c,ob8,ob9];

# ╔═╡ ba994d45-f992-4d07-99c2-a3ed5efe1743
deltas = (1.0mm, 1.0mm, 1.0mm);

# ╔═╡ 1842125f-9dfc-42bc-b524-f2bad4a26021
dims = (400, 400, 30) # odd

# ╔═╡ e6f1406e-d281-4b9c-868a-35e9accc6110
image_of_phantom = phantom(axes(ig1)..., objects);

# ╔═╡ 9d5d03cd-b858-40df-b751-7fca88f853d8
ig1 = ImageGeom(MaskCircle(); dims = dims, deltas = deltas);

# ╔═╡ 1a0a57ab-f451-41fd-addb-c554de056919
md"""
## Define the CT Geometry
"""

# ╔═╡ 832e2c99-0a18-4839-b029-6612a9a6dcd0
begin
	p1 = (ns = 150, ds = 2mm, nt = 60, dt = 2mm, na = 200, dsd = 2000mm, dod = 400mm)
	rg1 = CtFanArc(; p1...)
end;

# ╔═╡ 365504bd-4234-48aa-84af-31e3b4862262
ct_geom_plot3(rg1, ig1)

# ╔═╡ a673ea08-3121-4455-91f0-c9c3801dc54d
proj_arc1 = radon(rays(rg1), objects);

# ╔═╡ 2fcd7aee-8512-47c7-b7eb-f69fb5598d8f
begin
	plan1 = plan_fbp(rg1, ig1; window = Window(Hamming(), 1.0))
	ct_scan = fdk(plan1, proj_arc1)
end;

# ╔═╡ 8b1f3ea9-084f-4194-9293-9b8472a593f7
@bind a1 PlutoUI.Slider(axes(ct_scan, 3), show_value=true, default=20)

# ╔═╡ e1070f91-acd2-4218-a97c-564b8bdc579c
let
	f = Figure(resolution=(1600, 800))

	ax = CairoMakie.Axis(
		f[1, 1],
		title="phantom"
	)
	CairoMakie.heatmap!(ustrip.(image_of_phantom[:, :, a1]); colormap=:grays)
	
	ax = CairoMakie.Axis(
		f[1, 2],
		title="sinogram"
	)
	CairoMakie.heatmap!(ustrip.(ct_scan[:, :, a1]); colormap=:grays)
	
	f
end

[JeffFessler]

• How can I make this sinogram better? Are there parameters that would work better for the phantom I have?

The right image above is the reconstructed image, it is not the sinogram! The sinogram is proj_arc1 in your code.
You could increase na to get finer angular sampling to reduce aliasing artifacts. You could also decrease the detector pixel size ds to have finer radial sampling (the increase ns to keep the detector FOV big enough).

change the density

You already did with -200 and 40 etc. as the last parameter.

• The CT energy seems to be hardcoded into the package

No it is not. The code just computes line integrals of the objects you provide.
You can set the object value to be whatever you want, e.g., the linear attenuation coefficient at the energy of interest to you.
You could even do a polyenergetic simulation if you know your X-ray spectrum by integrating numerically over the spectrum.

[JeffFessler]

The user gets to choose the units for the object values, and if you use -200 then you are probably using Hounsfield units. But you certainly can instead use linear attenuation coefficients corresponding to the materials and photon energy of interest to you.

When you are using ImagePhantoms you are doing simulations, so you can make up any energy spectrum you want.
I have a few spectra in this matlab toolbox: https://github.com/JeffFessler/mirt/tree/main/ct/xray-spectra
You'd have to search the literature to find out how to measure spectra.

[Dale-Black]

Gotcha! Is there an example or some source code that shows how to plug in a spectra to the simulation?

[Dale-Black]

If not, I am sure I can figure out what the numerical integration would look like. But just to be sure I am on the right track, the bulk of this would still be within the realm of the phantom design itself? Where the phantom materials are designed for specific kVs from whatever range we choose (e.g. 20:10:200 kV) and then sinograms for each one of those phantoms would be calculated and some sort of weighted sum of these sinograms would be the final product?

[JeffFessler]

There is a matlab example here that is pretty complicated:
https://github.com/JeffFessler/mirt/blob/main/example/xct_bh_example.m
I don't know any Julia code for it.
That code is based on the math this paper: https://doi.org/10.1117/12.467189
It would be great for someone to make Julia version!

[Dale-Black]

Hello Dr. Fessler,

I am finally getting close to setting up a "proper" polyenergetic spectrum for dual energy reconstruction. Does this look like it's on the right track? I am using the spectrum.x files from your mirt repo.

[singhtejus]

I've made some more progress in making my phantom.

Here's what I have so far:

Im trying to make the phantom as close to this:

Do have any suggestions on how I can get closer to that geometry (for example the lungs are currently simple cylinders) ?

[JeffFessler]

I doubt if anyone knows how to compute efficiently the analytical line integrals through those strange shapes. You could just make a voxelized phantom from the picture itself. You have to choose between making a phantom out of a few "simple" objects like ellipses that have simple Radon transforms, or out of numerous voxels (that also have simple Radon transforms but the forward project time would be much longer).