<div dir="ltr"><div><div><div>Hi,<br></div>Again, I insist, use the mailing list.<br></div>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).<br></div><div>Thanks,<br></div><div>Simon<br></div><br></div><div class="gmail_extra"><br><div class="gmail_quote">On Tue, Jul 11, 2017 at 3:55 PM, Daniel Markel <span dir="ltr"><<a href="mailto:markeldaniel@gmail.com" target="_blank">markeldaniel@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Hi Simon here is the full script, minus some sub-functions. I apologize for the mess in advance.<div><br></div><div>-Dan<br><div><br></div><div><div>def main():</div><div> folderpath = 'D:\test_object'</div><div> reader = srtk.ImageSeriesReader() #initialize the reader</div><div> onlyfile = srtk.ImageSeriesReader_<wbr>GetGDCMSeriesFileNames(<wbr>folderpath)</div><div> Image = srtk.Image() #initialize the image object type</div><div> reader.SetFileNames(onlyfile) #input the file names to be read in the series</div><div> Image = reader.Execute() #import the dicom slices into a 3D volume</div><div><br></div><div> Image = Image+2000</div><div> Image.SetDirection([1,0,0,0,0,<wbr>1,0,1,0])</div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div> size2 = Image.GetSize()</div><div> spacing2 = [1,1,1]</div><div> numberOfProjections = 31</div><div> Disp = np.zeros((numberOfProjections-<wbr>1,2,20))#recorded residual geometric error</div><div> Normccs = np.zeros((numberOfProjections-<wbr>1,20))#recorded metric values post-registration</div><div> shiftsx = np.zeros((numberOfProjections-<wbr>1,20))#recorded image corrections</div><div> shiftsy = np.zeros((numberOfProjections-<wbr>1,20))</div><div> Inplaneshiftsx = np.zeros((numberOfProjections-<wbr>1,20)) #recorded perturbations to be corrected as projected into the plane of the kV panel for x and y</div><div> Inplaneshiftsy = np.zeros((numberOfProjections-<wbr>1,20))</div><div> shiftsx2 = np.zeros((numberOfProjections-<wbr>1,20))</div><div> shiftsy2 = np.zeros((numberOfProjections-<wbr>1,20))</div><div> shiftsz2 = np.zeros((numberOfProjections-<wbr>1,20))</div><div> #Acc = np.zeros(numberOfProjections-<wbr>1)</div><div> cast = srtk.CastImageFilter()</div><div> cast.SetOutputPixelType(7)</div><div> Image2 = cast.Execute(Image)</div><div> firstAngle = 0</div><div> angularArc = 360</div><div><br></div><div> sid = 10000 # source to isocenter distance in mm</div><div> sdd = 15360 # source to detector distance in mm</div><span class=""><div><br></div><div> isox = 0 # X coordinate on the projection image of isocenter</div></span><div> isoy = -90 # Y coordinate on the projection image of isocenter\</div><div><br></div><div> for k in range(0,1):</div><div> Image = Image2</div><div> geometry = srtk.<wbr>ThreeDCircularProjectionGeomet<wbr>ry()</div><div> angle = 180+k*18</div><div> geometry.AddProjection(sid,<wbr>sdd,angle,isox,isoy)</div><div><br></div><div> projInput = srtk.ConstantImageSource()</div><span class=""><div><br></div><div> </div><div><br></div><div> size = (425,425,1)</div><div> spacing = (1,1,1)</div><div> projInput.SetOrigin(origin)</div><div> projInput.SetSpacing(spacing)</div><div> projInput.SetSize(size)</div><div> projInput.SetConstant(int(1))</div><div> pixelID = Image.GetPixelIDValue()</div><div> projInput.SetOutputPixelType(<wbr>pixelID)</div><div> proj = projInput.Execute()</div></span><div> iters2 = numberOfProjections - 1</div><div><br></div><div> </div><div><br></div><div> </div><div><br></div><div> shifty2 = np.zeros(iters2+1)</div><div> #shiftx2 = np.zeros(iters2+1)</div><div> shiftz2 = np.zeros(iters2+1)</div><div><br></div><div> shiftx2[1:] = np.linspace(1,30,30)-15</div><div><br></div><div><br></div><div> Temp = srtk.GetArrayFromImage(Image)</div><div> Image_shifted_array = np.zeros((size[0],size[1],<wbr>iters2+1))</div><div><br></div><div> for i in range(0,iters2):</div><div><br></div><div> Temp_vol = np.roll(Temp,int(shiftx2[i]),<wbr>2)#Shift in the LR axis </div><div> Temp_vol = np.roll(Temp_vol,int(shifty2[<wbr>i]),1)#Shift in the AP axis</div><div> Temp_vol = np.roll(Temp_vol,int(shiftz2[<wbr>i]),0)#Shift in the I/S axis</div><div> shiftsx2[i,k] = shiftx2[i]</div><div> shiftsy2[i,k] = shifty2[i]</div><div> shiftsz2[i,k] = shiftz2[i]</div><div> Inplaneshiftsy[i,k] = shiftz2[i]*1.25</div><div> Inplaneshiftsx[i,k] = ((shiftx2[i]*0.39063)*np.cos(<wbr>angle*np.pi/180)+(shifty2[i]*<wbr>0.39063)*np.sin(angle*np.pi/<wbr>180))</div><span class=""><div> fwrd = srtk.<wbr>JosephForwardProjectionImageFi<wbr>lter()</div><div> projInput = srtk.ConstantImageSource()</div><div> geometry = srtk.<wbr>ThreeDCircularProjectionGeomet<wbr>ry()</div><div> </div><div><br></div><div> size = (425,425,1)</div><div> spacing = (1,1,1)</div><div> projInput.SetOrigin(origin)</div><div> projInput.SetSpacing(spacing)</div><div> projInput.SetSize(size)</div><div> projInput.SetConstant(int(1))</div><div> pixelID = Image.GetPixelIDValue()</div><div> projInput.SetOutputPixelType(<wbr>pixelID)</div><div> proj = projInput.Execute()</div><div> Image = srtk.GetImageFromArray(Temp_<wbr>vol)</div><div> Image.SetOrigin(origin)</div><div> Image.SetSpacing(spacing2)</div><div> Image.SetDirection([1,0,0,0,0,<wbr>1,0,1,0])</div></span><div> #geometry.AddProjection()</div><span class=""><div> geometry.AddProjection(sid,<wbr>sdd,angle,isox,isoy)</div><div> fwrd.SetGeometry(geometry)</div><div> proj2 = fwrd.Execute(proj,Image)</div></span><div> Image_shifted_array[:,:,i] = srtk.GetArrayFromImage(proj2[:<wbr>,:,0])</div><div><br></div><div> Image_unshifted = Image_shifted_array[:,:,0]</div><div> </div><div> for i in range(0,iters2-1):</div><div><br></div><div> Image_shifted = np.squeeze(Image_shifted_<wbr>array[:,:,i])</div><div><br></div><div> plt.imshow(Image_shifted_<wbr>array[:,:,i], cmap = cm.Greys_r)</div><div> plt.show()</div><div> plt.imshow(Image_unshifted, cmap = cm.Greys_r)</div><div> plt.show()</div><div><br></div><div> iters = 150</div><div> step_size = 30</div><div> randy = np.round((np.random.rand(<wbr>iters)-0.5)*step_size)</div><div> randx = np.round((np.random.rand(<wbr>iters)-0.5)*step_size)</div><div><br></div><div> normcc_last = NCC(Image_unshifted[40:(425-<wbr>40),40:(425-40)],Image_<wbr>shifted[40:(425-40),40:(425-<wbr>40)])</div><div><br></div><div> shiftx = 0</div><div> shifty = 0</div><div> time_start = time.clock()</div><div> for shift in range(iters):</div><div> Image_temp = np.roll(Image_shifted,int(<wbr>shiftx+randx[shift]),1)</div><div> Image_temp = np.roll(Image_temp,int(shifty+<wbr>randy[shift]),0)</div><div><br></div><div> normcc_new = NCC(Image_unshifted[40:(425-<wbr>40),40:(425-40)],Image_temp[<wbr>40:(425-40),40:(425-40)])</div><div><br></div><div> if normcc_new>normcc_last:</div><div> shiftx = shiftx+randx[shift]</div><div> shifty = shifty+randy[shift]</div><div> normcc_last = normcc_new</div><div> if normcc_new == 0.5:</div><div> </div><div> break</div><div><br></div><div> Disp[i,0,k] = np.sqrt((Inplaneshiftsx[i,k]*(<wbr>1536/1000))**2+(<wbr>Inplaneshiftsy[i,k]*(1536/<wbr>1000))**2)</div><div> Disp[i,1,k] = np.sqrt((shiftx*1.66 - (Inplaneshiftsx[i,k]*(1536/<wbr>1000)))**2+(shifty*1.66 - (Inplaneshiftsy[i,k]*(1536/<wbr>1000)))**2)</div><div> Normccs[i,k] = normcc_last</div><div> shiftsx[i,k] = shiftx</div><div> shiftsy[i,k] = shifty</div><div><br></div><div> ind = np.argsort(Disp,axis = 0)</div><div> Disp2 = Disp</div><div> Normccs2 = Normccs</div><div> Inplaneshiftsx_new = Inplaneshiftsx</div><div> Inplaneshiftsy_new = Inplaneshiftsy</div><div> shiftsx_new = shiftsx</div><div> shiftsy_new = shiftsy</div><div> shiftsx2_new = shiftsx2</div><div> shiftsy2_new = shiftsy2</div><div> shiftsz2_new = shiftsz2</div><div> for p in range(0,20):</div><div> Disp2[:,0,p] = Disp[ind[:,0,p],0,p]</div><div> Disp2[:,1,p] = Disp[ind[:,0,p],1,p]</div><div> Normccs2[:,p] = Normccs2[ind[:,0,p],p]</div><div> Inplaneshiftsy_new[:,p] = Inplaneshiftsy[ind[:,0,p],p]</div><div> Inplaneshiftsx_new[:,p] = Inplaneshiftsx[ind[:,0,p],p]</div><div> shiftsy_new[:,p] = shiftsy[ind[:,0,p],p]</div><div> shiftsx_new[:,p] = shiftsx[ind[:,0,p],p]</div><div> shiftsy2_new[:,p] = shiftsy2[ind[:,0,p],p]</div><div> shiftsx2_new[:,p] = shiftsx2[ind[:,0,p],p]</div><div> shiftsz2_new[:,p] = shiftsz2[ind[:,0,p],p]</div><div> with open('D:\\NCC_result_centered_<wbr>HD_systematic2.csv','w', newline = '') as csvfile:</div><div> spamwriter = csv.writer(csvfile, delimiter=' ', quotechar='|', quoting=csv.QUOTE_MINIMAL)</div><div> sz = np.shape(Disp2)</div><div> for t in range(0,sz[2]):</div><div> for l in range(0,sz[0]):</div><div> line = ''</div><div> t0 = str(180+t*18)</div><div> t1 = ' '</div><div> t2 = str(Disp2[l,0,t])</div><div> t3 = str(Disp2[l,1,t])</div><div> t4 = str(Inplaneshiftsx_new[l,t])</div><div> t5 = str(Inplaneshiftsy_new[l,t])</div><div> t6 = str(shiftsx2_new[l,t])</div><div> t7 = str(shiftsy2_new[l,t])</div><div> t8 = str(shiftsz2_new[l,t])</div><div> t9 = str(shiftsx_new[l,t])</div><div> t10 = str(shiftsy_new[l,t])</div><div> t11 = 'angle:'</div><div> t12 = str(Normccs2[l,t])</div><div> </div><div> </div><div> line = t11+t0+t1+t2+t1+t3+t1+t4+t1+<wbr>t5+t1+t6+t1+t7+t1+t8+t1+t9+t1+<wbr>t10+t1+t12</div><div><br></div><div> </div><div> </div><div> spamwriter.writerow(line)</div><div> </div><div> </div><div>main()</div></div></div></div><div class="HOEnZb"><div class="h5"><div class="gmail_extra"><br><div class="gmail_quote">On Tue, Jul 11, 2017 at 2:03 AM, Simon Rit <span dir="ltr"><<a href="mailto:simon.rit@creatis.insa-lyon.fr" target="_blank">simon.rit@creatis.insa-lyon.<wbr>fr</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div><div><div>Hi,<br></div>Please use the mailing list.<br></div>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?<span class="m_-1734077975825938577HOEnZb"><font color="#888888"><br></font></span></div><span class="m_-1734077975825938577HOEnZb"><font color="#888888">Simon<br></font></span></div><div class="m_-1734077975825938577HOEnZb"><div class="m_-1734077975825938577h5"><div class="gmail_extra"><br><div class="gmail_quote">On Mon, Jul 10, 2017 at 11:19 PM, Daniel Markel <span dir="ltr"><<a href="mailto:markeldaniel@gmail.com" target="_blank">markeldaniel@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><p>Hi Simon,</p><p><br></p><p> 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. </p><p><br></p><p>I'm getting projections using the following snippet of code</p><p><br></p><p>origin = [-200,-200,-180-(600-512)*0.39<wbr>063] </p><p>sid = 1000 # source to isocenter distance in mm<br>sdd = 1536 # source to detector distance in mm</p><p>isox = 0 # X coordinate on the projection image of isocenter<br>isoy = -90 # Y coordinate on the projection image of isocenter</p><p>angle = 180<br></p><p>Temp = srtk.GetArrayFromImage(Image)</p><p></p><div>Temp_vol = np.roll(Temp,int(shiftx2),2)#S<wbr>hift in the LR axis <br>Temp_vol = np.roll(Temp_vol,int(shifty2),<wbr>1)#Shift in the AP axis<br>Temp_vol = np.roll(Temp_vol,int(shiftz2),<wbr>0)#Shift in the I/S axis</div><p></p><p></p><div>fwrd = srtk.JosephForwardProjectionIm<wbr>ageFilter()<br>projInput = srtk.ConstantImageSource()<br> geometry = srtk.ThreeDCircularProjectionG<wbr>eometry()</div><div><div>size = (425,425,1)<br> spacing = (1,1,1)<br>projInput.SetOrigin(origin)<br>projInput.SetSpacing(spacing)<br>projInput.SetSize(size)<br>projInput.SetConstant(int(1))<br>pixelID = Image.GetPixelIDValue()<br>projInput.SetOutputPixelType(p<wbr>ixelID)<br>proj = projInput.Execute()</div><div>Image = srtk.GetImageFromArray(Temp_vo<wbr>l)<br>Image.SetOrigin(origin)<br>Image.SetSpacing(spacing2)<br>Image.SetDirection([1,0,0,0,0,<wbr>1,0,1,0])<br>geometry.AddProjection(sid,sdd<wbr>,angle,isox,isoy)<br>fwrd.SetGeometry(geometry)<br> proj2 = fwrd.Execute(proj,Image)</div></div><p></p><p>projim = srtk.GetArrayFromImage(proj2[:<wbr>,:,0])</p><p><br></p><p>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? </p><p><br></p><p>Any insight would be greatly appreciated.</p><span class="m_-1734077975825938577m_-5294988651364277486HOEnZb"><font color="#888888"><p><br></p><p>-Daniel</p></font></span></div>
</blockquote></div><br></div>
</div></div></blockquote></div><br></div>
</div></div></blockquote></div><br></div>