[Rtk-users] Obtaining Cone Beam projections using sRTK
Simon Rit
simon.rit at creatis.insa-lyon.fr
Wed Jul 12 05:08:22 EDT 2017
Hi,
Again, I insist, use the mailing list.
Can you make it a self contained code? You can replace your image with a
simulated Shepp Logan phantom using DrawSheppLoganFilter. And please make
it more minimal (just one shift).
Thanks,
Simon
On Tue, Jul 11, 2017 at 3:55 PM, Daniel Markel <markeldaniel at gmail.com>
wrote:
> Hi Simon here is the full script, minus some sub-functions. I apologize
> for the mess in advance.
>
> -Dan
>
> def main():
> folderpath = 'D:\test_object'
> reader = srtk.ImageSeriesReader() #initialize the reader
> onlyfile = srtk.ImageSeriesReader_GetGDCMSeriesFileNames(folderpath)
> Image = srtk.Image() #initialize the image object type
> reader.SetFileNames(onlyfile) #input the file names to be read in the
> series
> Image = reader.Execute() #import the dicom slices into a 3D volume
>
> Image = Image+2000
> Image.SetDirection([1,0,0,0,0,1,0,1,0])
>
>
>
>
>
> size2 = Image.GetSize()
> spacing2 = [1,1,1]
> numberOfProjections = 31
> Disp = np.zeros((numberOfProjections-1,2,20))#recorded residual
> geometric error
> Normccs = np.zeros((numberOfProjections-1,20))#recorded metric values
> post-registration
> shiftsx = np.zeros((numberOfProjections-1,20))#recorded image
> corrections
> shiftsy = np.zeros((numberOfProjections-1,20))
> Inplaneshiftsx = np.zeros((numberOfProjections-1,20)) #recorded
> perturbations to be corrected as projected into the plane of the kV panel
> for x and y
> Inplaneshiftsy = np.zeros((numberOfProjections-1,20))
> shiftsx2 = np.zeros((numberOfProjections-1,20))
> shiftsy2 = np.zeros((numberOfProjections-1,20))
> shiftsz2 = np.zeros((numberOfProjections-1,20))
> #Acc = np.zeros(numberOfProjections-1)
> cast = srtk.CastImageFilter()
> cast.SetOutputPixelType(7)
> Image2 = cast.Execute(Image)
> firstAngle = 0
> angularArc = 360
>
> sid = 10000 # source to isocenter distance in mm
> sdd = 15360 # source to detector distance in mm
>
> isox = 0 # X coordinate on the projection image of isocenter
> isoy = -90 # Y coordinate on the projection image of isocenter\
>
> for k in range(0,1):
> Image = Image2
> geometry = srtk.ThreeDCircularProjectionGeometry()
> angle = 180+k*18
> geometry.AddProjection(sid,sdd,angle,isox,isoy)
>
> projInput = srtk.ConstantImageSource()
>
>
>
> size = (425,425,1)
> spacing = (1,1,1)
> projInput.SetOrigin(origin)
> projInput.SetSpacing(spacing)
> projInput.SetSize(size)
> projInput.SetConstant(int(1))
> pixelID = Image.GetPixelIDValue()
> projInput.SetOutputPixelType(pixelID)
> proj = projInput.Execute()
> iters2 = numberOfProjections - 1
>
>
>
>
>
> shifty2 = np.zeros(iters2+1)
> #shiftx2 = np.zeros(iters2+1)
> shiftz2 = np.zeros(iters2+1)
>
> shiftx2[1:] = np.linspace(1,30,30)-15
>
>
> Temp = srtk.GetArrayFromImage(Image)
> Image_shifted_array = np.zeros((size[0],size[1],iters2+1))
>
> for i in range(0,iters2):
>
> Temp_vol = np.roll(Temp,int(shiftx2[i]),2)#Shift in the LR
> axis
> Temp_vol = np.roll(Temp_vol,int(shifty2[i]),1)#Shift in the
> AP axis
> Temp_vol = np.roll(Temp_vol,int(shiftz2[i]),0)#Shift in the
> I/S axis
> shiftsx2[i,k] = shiftx2[i]
> shiftsy2[i,k] = shifty2[i]
> shiftsz2[i,k] = shiftz2[i]
> Inplaneshiftsy[i,k] = shiftz2[i]*1.25
> Inplaneshiftsx[i,k] = ((shiftx2[i]*0.39063)*np.cos(
> angle*np.pi/180)+(shifty2[i]*0.39063)*np.sin(angle*np.pi/180))
> fwrd = srtk.JosephForwardProjectionImageFilter()
> projInput = srtk.ConstantImageSource()
> geometry = srtk.ThreeDCircularProjectionGeometry()
>
>
> size = (425,425,1)
> spacing = (1,1,1)
> projInput.SetOrigin(origin)
> projInput.SetSpacing(spacing)
> projInput.SetSize(size)
> projInput.SetConstant(int(1))
> pixelID = Image.GetPixelIDValue()
> projInput.SetOutputPixelType(pixelID)
> proj = projInput.Execute()
> Image = srtk.GetImageFromArray(Temp_vol)
> Image.SetOrigin(origin)
> Image.SetSpacing(spacing2)
> Image.SetDirection([1,0,0,0,0,1,0,1,0])
> #geometry.AddProjection()
> geometry.AddProjection(sid,sdd,angle,isox,isoy)
> fwrd.SetGeometry(geometry)
> proj2 = fwrd.Execute(proj,Image)
> Image_shifted_array[:,:,i] = srtk.GetArrayFromImage(proj2[:
> ,:,0])
>
> Image_unshifted = Image_shifted_array[:,:,0]
>
> for i in range(0,iters2-1):
>
> Image_shifted = np.squeeze(Image_shifted_array[:,:,i])
>
> plt.imshow(Image_shifted_array[:,:,i], cmap = cm.Greys_r)
> plt.show()
> plt.imshow(Image_unshifted, cmap = cm.Greys_r)
> plt.show()
>
> iters = 150
> step_size = 30
> randy = np.round((np.random.rand(iters)-0.5)*step_size)
> randx = np.round((np.random.rand(iters)-0.5)*step_size)
>
> normcc_last = NCC(Image_unshifted[40:(425-
> 40),40:(425-40)],Image_shifted[40:(425-40),40:(425-40)])
>
> shiftx = 0
> shifty = 0
> time_start = time.clock()
> for shift in range(iters):
> Image_temp = np.roll(Image_shifted,int(
> shiftx+randx[shift]),1)
> Image_temp = np.roll(Image_temp,int(shifty+
> randy[shift]),0)
>
> normcc_new = NCC(Image_unshifted[40:(425-
> 40),40:(425-40)],Image_temp[40:(425-40),40:(425-40)])
>
> if normcc_new>normcc_last:
> shiftx = shiftx+randx[shift]
> shifty = shifty+randy[shift]
> normcc_last = normcc_new
> if normcc_new == 0.5:
>
> break
>
> Disp[i,0,k] = np.sqrt((Inplaneshiftsx[i,k]*(1536/1000))**2+(
> Inplaneshiftsy[i,k]*(1536/1000))**2)
> Disp[i,1,k] = np.sqrt((shiftx*1.66 -
> (Inplaneshiftsx[i,k]*(1536/1000)))**2+(shifty*1.66 -
> (Inplaneshiftsy[i,k]*(1536/1000)))**2)
> Normccs[i,k] = normcc_last
> shiftsx[i,k] = shiftx
> shiftsy[i,k] = shifty
>
> ind = np.argsort(Disp,axis = 0)
> Disp2 = Disp
> Normccs2 = Normccs
> Inplaneshiftsx_new = Inplaneshiftsx
> Inplaneshiftsy_new = Inplaneshiftsy
> shiftsx_new = shiftsx
> shiftsy_new = shiftsy
> shiftsx2_new = shiftsx2
> shiftsy2_new = shiftsy2
> shiftsz2_new = shiftsz2
> for p in range(0,20):
> Disp2[:,0,p] = Disp[ind[:,0,p],0,p]
> Disp2[:,1,p] = Disp[ind[:,0,p],1,p]
> Normccs2[:,p] = Normccs2[ind[:,0,p],p]
> Inplaneshiftsy_new[:,p] = Inplaneshiftsy[ind[:,0,p],p]
> Inplaneshiftsx_new[:,p] = Inplaneshiftsx[ind[:,0,p],p]
> shiftsy_new[:,p] = shiftsy[ind[:,0,p],p]
> shiftsx_new[:,p] = shiftsx[ind[:,0,p],p]
> shiftsy2_new[:,p] = shiftsy2[ind[:,0,p],p]
> shiftsx2_new[:,p] = shiftsx2[ind[:,0,p],p]
> shiftsz2_new[:,p] = shiftsz2[ind[:,0,p],p]
> with open('D:\\NCC_result_centered_HD_systematic2.csv','w', newline =
> '') as csvfile:
> spamwriter = csv.writer(csvfile, delimiter=' ', quotechar='|',
> quoting=csv.QUOTE_MINIMAL)
> sz = np.shape(Disp2)
> for t in range(0,sz[2]):
> for l in range(0,sz[0]):
> line = ''
> t0 = str(180+t*18)
> t1 = ' '
> t2 = str(Disp2[l,0,t])
> t3 = str(Disp2[l,1,t])
> t4 = str(Inplaneshiftsx_new[l,t])
> t5 = str(Inplaneshiftsy_new[l,t])
> t6 = str(shiftsx2_new[l,t])
> t7 = str(shiftsy2_new[l,t])
> t8 = str(shiftsz2_new[l,t])
> t9 = str(shiftsx_new[l,t])
> t10 = str(shiftsy_new[l,t])
> t11 = 'angle:'
> t12 = str(Normccs2[l,t])
>
>
> line = t11+t0+t1+t2+t1+t3+t1+t4+t1+
> t5+t1+t6+t1+t7+t1+t8+t1+t9+t1+t10+t1+t12
>
>
>
> spamwriter.writerow(line)
>
>
> main()
>
> On Tue, Jul 11, 2017 at 2:03 AM, Simon Rit <simon.rit at creatis.insa-lyon.fr
> > wrote:
>
>> Hi,
>> Please use the mailing list.
>> A 1 mm shift of an object at the isocenter should indeed translate in a
>> 1.536 mm shift on the projection with the defined geometry. I don't see
>> anything wrong in your code but you have only provided a part of it. Could
>> you provide a full self-contained script?
>> Simon
>>
>> On Mon, Jul 10, 2017 at 11:19 PM, Daniel Markel <markeldaniel at gmail.com>
>> wrote:
>>
>>> Hi Simon,
>>>
>>>
>>> Thank you for helping me earlier in getting SRTK up and
>>> running. I've had some success running it to attain some projections using
>>> a 3D volume but I have run into a problem and cannot find the solution
>>> anywhere online. Perhaps you can help.
>>>
>>>
>>> I'm getting projections using the following snippet of code
>>>
>>>
>>> origin = [-200,-200,-180-(600-512)*0.39063]
>>>
>>> sid = 1000 # source to isocenter distance in mm
>>> sdd = 1536 # source to detector distance in mm
>>>
>>> isox = 0 # X coordinate on the projection image of isocenter
>>> isoy = -90 # Y coordinate on the projection image of isocenter
>>>
>>> angle = 180
>>>
>>> Temp = srtk.GetArrayFromImage(Image)
>>>
>>>
>>> Temp_vol = np.roll(Temp,int(shiftx2),2)#Shift in the LR axis
>>> Temp_vol = np.roll(Temp_vol,int(shifty2),1)#Shift in the AP axis
>>> Temp_vol = np.roll(Temp_vol,int(shiftz2),0)#Shift in the I/S axis
>>>
>>>
>>>
>>>
>>> fwrd = srtk.JosephForwardProjectionImageFilter()
>>> projInput = srtk.ConstantImageSource()
>>> geometry = srtk.ThreeDCircularProjectionGeometry()
>>> size = (425,425,1)
>>> spacing = (1,1,1)
>>> projInput.SetOrigin(origin)
>>> projInput.SetSpacing(spacing)
>>> projInput.SetSize(size)
>>> projInput.SetConstant(int(1))
>>> pixelID = Image.GetPixelIDValue()
>>> projInput.SetOutputPixelType(pixelID)
>>> proj = projInput.Execute()
>>> Image = srtk.GetImageFromArray(Temp_vol)
>>> Image.SetOrigin(origin)
>>> Image.SetSpacing(spacing2)
>>> Image.SetDirection([1,0,0,0,0,1,0,1,0])
>>> geometry.AddProjection(sid,sdd,angle,isox,isoy)
>>> fwrd.SetGeometry(geometry)
>>> proj2 = fwrd.Execute(proj,Image)
>>>
>>>
>>>
>>> projim = srtk.GetArrayFromImage(proj2[:,:,0])
>>>
>>>
>>> What I wanted to do was shift an object in the 3D volume 'temp' by a
>>> known amount recreate the image object using the shifted volume and then
>>> create a kV projection through the shifted volume. I'm doing this to test
>>> the ability of detecting shifts in the volume using the kV projected
>>> images. The problem is that I was calculating the anticipated shifts using
>>> a cone beam geometry where the shift in the 3D volume is magnified by an
>>> amount related to the distance of the detector from the source and the
>>> position of the 3D object in relation to the source. For the geometry I
>>> have above it's 1000 mm to the center of my 3D object and 1536 mm from the
>>> source to the detector so the factor of magnification of any shift in my 3D
>>> object should be 1.536 on my projection but when I look at my results the
>>> magnification is always 1 as if I'm not using a cone beam geometry at all
>>> but simply attaining a projection straight on using parallel x-rays. Am I
>>> missing something in my code to get a cone geometry?
>>>
>>>
>>> Any insight would be greatly appreciated.
>>>
>>>
>>> -Daniel
>>>
>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/rtk-users/attachments/20170712/897ba48a/attachment-0001.html>
More information about the Rtk-users
mailing list