Comments (7)
from rtk.
Hi,
Thank you so much for your quick answer!
I didn't use SimpleRTK to generate my projections. I have attached them here as a tiff stack so you can look at them. I load them as a stack of files through another software that switches them into a ndarray.
I have 500 tif files of 256 x 256 and after the numpy conversion I have a ndarray (ndaData) of 500, 256, 256 (Z, Y, X convention). Maybe the axes convention causes problems?
I will still try your advice, to generate my projections with SimpleRTK methods and to compare the projections. I tried to reconstruct with the same code as I have now with another set of Cone Beam projections generated from a real CT machine and I still have unexpected results...
Thank you so much,
Emmanuelle
SheppLoganConeBeam.zip
from rtk.
Hi,
Just wanted to give an update: I generated projections with RTK (with the command line options in the built applications : http://wiki.openrtk.org/index.php/RTK/Scripts/FDK) and I obtain very similar results. And when I test the reconstruction on the projections made by the command line methods, I obtain basically the same reconstruction output with lots of weird artefacts.
But, when I input either negative values to my sad and sid distances (-20000 and -10000) or when I put my max angle negative (-360), I obtain pretty good results (see image below). I suppose there is a flaw in the way I create my geometry but I don't have any idea what it is!
Please tell me if you have an idea about what it is that's causing this!
Have a good day and thank you,
Emmanuelle
from rtk.
Thanks for the update. That's weird indeed. Are you writing it to disk in mhd format and reading with SimpleRTK? I don't understand how this is possible to be honest. Please give the command line and the Python code to read it.
from rtk.
Hi,
So I'm getting closer to an answer. I was able to create projections with these command line options:
C:\Users\admin\Desktop\RTK-BuildExamples>rtksimulatedgeometry -n 256 -o C:\Users\admin\Desktop\RTK-BuildExamples\ExternalData\testing\Data\Baseline\XRad\geometry.xml
C:\Users\admin\Desktop\RTK-BuildExamples>rtkprojectgeometricphantom -g C:\Users\admin\Desktop\RTK-BuildExamples\ExternalData\testing\Data\Baseline\XRad\geometry.xml -o projections.tif --spacing 2 --dimension 256 --phantomfile C:\Users\admin\Desktop\RTK-BuildExamples\ExternalData\testing\Data\Input\GeometricPhantom\SheppLogan.txt
and then reconstruct them with my simpleRTK code:
I am working with Dragonfly, I import the files in the software then with this code transform the channel containing my projections to a ndarray.
Needs a channel
channelGUIDList = (
WorkingContext.getEntitiesOfClass(self.getImplementation(), OrsSelectedObjects, Channel.getClassNameStatic()))
if len(channelGUIDList) == 0:
return
self.channel = orsObj(channelGUIDList[0]) # Take the first one
numberOfProjections = self.channel.getZSize()
self.sdd = 1536
self.sid = 1000
anglemin = 0 # first angle
anglemax = 360 # last angle
isox = 0
isoy = 0
outangle = 0
inangle = 0
srcx = 0
srcy = 0
for noProj in range(0, numberOfProjections):
angle = anglemin + noProj * anglemax / numberOfProjections
self.geometry.AddProjection(self.sid, self.sdd, angle, isox, isoy, outangle, inangle, srcx, srcy)
Convert from intensity to attenuation
ndaData = self.channel.getNDArray(0) # working with Dragonly, that’s how I a transform my projections images in ndarray
I0 = 2 ** 16 - 1
if ndaData.min() > 0:
intensityOffset = 0
else:
intensityOffset = 0.01*(ndaData.max()-ndaData.min())-ndaData.min() # offset to non negative numbers
ndaData = np.log(float(I0)/(ndaData+intensityOffset)) # Convert intensity to line integral
print("Attenuation range =", ndaData.min(), ndaData.max())
ndaData = ndaData.astype(np.float32)
must have the same dtype as the reconstructed image i.e. float32
Convert to a SimpleRTK image
srtkData = srtk.GetImageFromArray(ndaData)
uvSize = srtkData.GetSize()
convert to mm
spacing = [1000 * self.channel.getXSpacing(), 1000 * self.channel.getYSpacing(), 1000 * self.channel.getZSpacing()]
origin = [-0.5spacing[0](uvSize[0]-1), -0.5spacing[1](uvSize[1]-1), 0] # centered
srtkData.SetSpacing(spacing)
srtkData.SetOrigin(origin)
Create an empty output
constantImageRecons = srtk.ConstantImageSource()
sizeX = 256
sizeY = 256
sizeZ = 256
sizeOutput = [sizeX, sizeY, sizeZ]
sizeX = 1
sizeY = 1
sizeZ = 1
spacing = [sizeX, sizeY, sizeZ]
sizeX = -127.5
sizeY = -127.5
sizeZ = -127.5
origin = [sizeX, sizeY, sizeZ]
constantImageRecons.SetOrigin(origin)
constantImageRecons.SetSpacing(spacing)
constantImageRecons.SetSize(sizeOutput)
constantImageRecons.SetConstant(0.0)
imageRecons = constantImageRecons.Execute()
print("Performing reconstruction")
feldkamp = srtk.CudaFDKConeBeamReconstructionFilter()
feldkamp.SetGeometry(self.geometry)
feldkamp.SetTruncationCorrection(0.0)
feldkamp.SetHannCutFrequency(0.0)
srtkImage = feldkamp.Execute(imageRecons, srtkData)
Convert to ndarray
reconstructData = srtk.GetArrayFromImage(srtkImage)
Put it on a new channel
newChannel = Channel()
newChannel.resize(reconstructData.shape, dtype=reconstructData.dtype)
origin = Vector3() ; origin.setXYZ(float(self.ui.outputOriginX.text()), float(self.ui.outputOriginY.text()),
float(self.ui.outputOriginZ.text()))
newChannel.setOrigin(origin)
scale = srtkImage.GetSpacing()
newChannel.setXSpacing(scale[0] / 1000) # convert into meters
newChannel.setYSpacing(scale[1] / 1000)
newChannel.setZSpacing(scale[2] / 1000)
newChannel = createChannelFromNumpyArray(reconstructData, newChannel.getGUID())
save to files in directory
file = self.directory + r'\file.tif'
orsimagesaver.OrsImageSaver.exportDatasetToTiffFile(newChannel.getGUID(),
fileName=file,
useLZWCompression=False,
exportMultiFrameFile=False,
showProgress=False)
if self.ui.loadChannel.isChecked():
newChannel.setTitle('Reconstructions')
newChannel.publish()
print('Reconstruction complete.')
else:
print('Reconstruction files are now in the chosen directory.')
and obtain good results. When I was trying the same reconstruction code on the projections that I posted earlier in the tif stack, I was obtaining the weird artefacts on the reconstructed images. But when I inverted the signs of sad and sid, it worked. I figured out that when I inverted the Z axis in my projections and then applied the reconstruction code, it was all working out as well. I suppose there was a difference in the way I had created my projections. I will try to understand why I need to invert an axis to make my reconstruction code work with my own projections but I will close this issue since I found out the way to obtain good reconstructed images.
I'm sorry for not posting my code earlier but I wasn't sure if I was at liberty to do so.
Thank you for your help!
Emmanuelle
from rtk.
Hi,
I don't understand why you go through this path. I suggest that you do everything in Python with SimpleRTK following, e.g., the example here
http://wiki.openrtk.org/index.php/FanBeam
Then, you can try to understand what's wrong when you do your conversion via Dragon fly. My guess is that it is not using a C convention for arrays as ITK does but, maybe, a Fortran convention as numpy does. So, from ITK to numpy, the order [x,y,z] becomes [z,y,x] and it's possible that you can compensate that with some minus signs in your geometry.
Simon
from rtk.
Also, this is not an RTK issue but this is a user question, I would prefer to have that conversation on the mailing list [email protected] dedicated to this. I'm therefore closing the issue.
from rtk.
Related Issues (20)
- Read Raw Data HOT 1
- how to use "DisplacedDetectorForOffsetFieldOfViewImageFilter"? HOT 1
- how itk.CudaImage convert to itk.Image? HOT 1
- Python Wrapping fails for missing dimensions HOT 15
- Why use SetTruncationCorrection in RTK has no performance? HOT 1
- Why SetTruncationCorrection(0.4) is not useful? HOT 3
- "Please, is there any source code available for this function?" HOT 1
- How to output the image after truncation correction and ramp filtering? HOT 1
- reconstruction by cuda occur problem: ITK ERROR: CUDA ERROR: out of memory HOT 1
- “itkRTK-5.2.lib” is not created HOT 1
- Unable to successfully build the rtk library HOT 2
- Reconstructed slices are returning as zeros and NaN HOT 1
- Unable to successfully build the rtk library HOT 2
- Error C1001 An internal error has occurred in the compiler. HOT 3
- Cuda memory error when using Cuda forward projector HOT 11
- Unable to debug hxx files and h files HOT 8
- Unable to compile simpleRTK HOT 1
- Installation issue with pip HOT 1
- why rtk::CudaFDKConeBeamReconstructionFilter's GetRampFilter()->SetTruncationCorrection() failed!
- why rtk::CudaFDKConeBeamReconstructionFilter's GetRampFilter()->SetTruncationCorrection() failed! HOT 4
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from rtk.