<html xmlns:v="urn:schemas-microsoft-com:vml" xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40"><head><meta http-equiv=Content-Type content="text/html; charset=iso-8859-1"><meta name=Generator content="Microsoft Word 14 (filtered medium)"><style><!--
/* Font Definitions */
@font-face
{font-family:Wingdings;
panose-1:5 0 0 0 0 0 0 0 0 0;}
@font-face
{font-family:Wingdings;
panose-1:5 0 0 0 0 0 0 0 0 0;}
@font-face
{font-family:Calibri;
panose-1:2 15 5 2 2 2 4 3 2 4;}
@font-face
{font-family:Tahoma;
panose-1:2 11 6 4 3 5 4 4 2 4;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0cm;
margin-bottom:.0001pt;
font-size:12.0pt;
font-family:"Times New Roman","serif";}
a:link, span.MsoHyperlink
{mso-style-priority:99;
color:blue;
text-decoration:underline;}
a:visited, span.MsoHyperlinkFollowed
{mso-style-priority:99;
color:purple;
text-decoration:underline;}
p.MsoAcetate, li.MsoAcetate, div.MsoAcetate
{mso-style-priority:99;
mso-style-link:"Sprechblasentext Zchn";
margin:0cm;
margin-bottom:.0001pt;
font-size:8.0pt;
font-family:"Tahoma","sans-serif";}
p.MsoListParagraph, li.MsoListParagraph, div.MsoListParagraph
{mso-style-priority:34;
margin-top:0cm;
margin-right:0cm;
margin-bottom:0cm;
margin-left:36.0pt;
margin-bottom:.0001pt;
font-size:12.0pt;
font-family:"Times New Roman","serif";}
span.E-MailFormatvorlage17
{mso-style-type:personal-compose;
font-family:"Calibri","sans-serif";
color:windowtext;}
span.SprechblasentextZchn
{mso-style-name:"Sprechblasentext Zchn";
mso-style-priority:99;
mso-style-link:Sprechblasentext;
font-family:"Tahoma","sans-serif";
mso-fareast-language:DE;}
.MsoChpDefault
{mso-style-type:export-only;
font-family:"Calibri","sans-serif";
mso-fareast-language:EN-US;}
@page WordSection1
{size:612.0pt 792.0pt;
margin:70.85pt 70.85pt 2.0cm 70.85pt;}
div.WordSection1
{page:WordSection1;}
/* List Definitions */
@list l0
{mso-list-id:662242275;
mso-list-type:hybrid;
mso-list-template-ids:-69031092 855254180 67567619 67567621 67567617 67567619 67567621 67567617 67567619 67567621;}
@list l0:level1
{mso-level-start-at:0;
mso-level-number-format:bullet;
mso-level-text:-;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-18.0pt;
font-family:"Calibri","sans-serif";
mso-fareast-font-family:Calibri;
mso-bidi-font-family:"Times New Roman";}
@list l0:level2
{mso-level-number-format:bullet;
mso-level-text:o;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-18.0pt;
font-family:"Courier New";}
@list l0:level3
{mso-level-number-format:bullet;
mso-level-text:\F0A7;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-18.0pt;
font-family:Wingdings;}
@list l0:level4
{mso-level-number-format:bullet;
mso-level-text:\F0B7;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-18.0pt;
font-family:Symbol;}
@list l0:level5
{mso-level-number-format:bullet;
mso-level-text:o;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-18.0pt;
font-family:"Courier New";}
@list l0:level6
{mso-level-number-format:bullet;
mso-level-text:\F0A7;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-18.0pt;
font-family:Wingdings;}
@list l0:level7
{mso-level-number-format:bullet;
mso-level-text:\F0B7;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-18.0pt;
font-family:Symbol;}
@list l0:level8
{mso-level-number-format:bullet;
mso-level-text:o;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-18.0pt;
font-family:"Courier New";}
@list l0:level9
{mso-level-number-format:bullet;
mso-level-text:\F0A7;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-18.0pt;
font-family:Wingdings;}
ol
{margin-bottom:0cm;}
ul
{margin-bottom:0cm;}
--></style><!--[if gte mso 9]><xml>
<o:shapedefaults v:ext="edit" spidmax="1026" />
</xml><![endif]--><!--[if gte mso 9]><xml>
<o:shapelayout v:ext="edit">
<o:idmap v:ext="edit" data="1" />
</o:shapelayout></xml><![endif]--></head><body lang=DE link=blue vlink=purple><div class=WordSection1><p class=MsoNormal><i><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'>Sorry, and again to the mailing list.<o:p></o:p></span></i></p><p class=MsoNormal><i><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'>I just signed up.<o:p></o:p></span></i></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'><o:p> </o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'><o:p> </o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'>Hello Simon,<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'>thank you for the fast reply.<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif"'>„<b><i>Joseph's method samples the ray with one pixel per slice in the main direction but it does not compute the intersection of the ray with each voxel. <a href="http://www.ncbi.nlm.nih.gov/pubmed/4000088">Siddon's method</a> does that. In fig 2 of <a href="http://www3.cs.stonybrook.edu/~mueller/papers/ISBI_06_quality_2.pdf">[Xu and Mueller, 2006]</a>, Joseph is referred to "slice interpolated" and Siddon to "box-line-integrated</i></b>".<br><span style='color:black'>OK, thanks for that hint. I think it will maybe also has the same problem with divergent rays and missing voxels during the reconstruction. Do you have a link to a paper or source to this algorithm. So the algorithm may work <o:p></o:p></span></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'>As following: <o:p></o:p></span></p><p class=MsoListParagraph style='text-indent:-18.0pt;mso-list:l0 level1 lfo1'><![if !supportLists]><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'><span style='mso-list:Ignore'>-<span style='font:7.0pt "Times New Roman"'> </span></span></span><![endif]><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'>calculate intersection with the physical bounding volume (entry exit points)<o:p></o:p></span></p><p class=MsoListParagraph style='text-indent:-18.0pt;mso-list:l0 level1 lfo1'><![if !supportLists]><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'><span style='mso-list:Ignore'>-<span style='font:7.0pt "Times New Roman"'> </span></span></span><![endif]><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'>from entrypoint, determine the 4 voxels that surround this entrypoint and bilinear interpolate the value at this position and sum it up<o:p></o:p></span></p><p class=MsoListParagraph style='text-indent:-18.0pt;mso-list:l0 level1 lfo1'><![if !supportLists]><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'><span style='mso-list:Ignore'>-<span style='font:7.0pt "Times New Roman"'> </span></span></span><![endif]><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'>go to next plane (plane that is most “perpendicular” to the current center ray (focus to detector center) ?<o:p></o:p></span></p><p class=MsoListParagraph style='text-indent:-18.0pt;mso-list:l0 level1 lfo1'><![if !supportLists]><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'><span style='mso-list:Ignore'>-<span style='font:7.0pt "Times New Roman"'> </span></span></span><![endif]><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'>at the end, the sum is normalized by the ray length ? (length(exitpoint – entrypoint))<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'>Is that right ?<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'>How can Joseph’s method be used for back projection ?<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif"'>“<b><i>This sounds very interesting, don't hesitate to share the code and/or the publication! BTW, what is DDA</i></b>”<span style='color:#1F497D'><o:p></o:p></span></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'>DDA stands for digital differential analyzer. I use this approach for a voxel-based forward projector. When I started my project I was looking for a fast and easy forward projector on voxel basis.<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'>I found an article about raytracing and acceleration structures where these people used this 3d dda. They actually needed it to traverse the bounding hierarchy to get to know what geometry<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'>the ray intersects with. I slightly modified it to get the intersecetion length of a ray within a voxel by substraction of the current and previous step width. I’v attached a zip file that contains this<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'>modification and the original source code. The original source code is from <a href="http://www.scratchapixel.com">www.scratchapixel.com</a> and there is also an article about this topic (</span><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:green;background:white;mso-highlight:white;mso-fareast-language:EN-US'><a href="http://www.scratchapixel.com/lessons/3d-basic-lessons/lesson-12-introduction-to-acceleration-structures/what-else/">http://www.scratchapixel.com/lessons/3d-basic-lessons/lesson-12-introduction-to-acceleration-structures/what-else/</a></span><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D;mso-fareast-language:EN-US'>).</span><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:green;mso-fareast-language:EN-US'><o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'>If you have a look at my code you’ll maybe miss some matrix calculation stuff etc. I use a scene based approach where I place the focus, object and detector in a so called record. Each record represents the scene geometry at a certain time, when a projection image was taken.<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'>Before reconstruction starts I calculate all these positions. I thought this could be a good approach to decouple the actuall reconstruction algorithm from the scene geometry.<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'><o:p> </o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'><o:p> </o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif"'>“<b><i>But generally we use matching resolution between pixels and voxels so the problem is minimal.</i></b>”<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif"'>Do you mean, that you use a volume resolution that fits to the current geometry setting and the detector’s pixel resolution ?<span style='color:black'><o:p></o:p></span></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'><o:p> </o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif"'>“<b><i>Thanks for the last trick, I am aware of it (<a href="http://dx.doi.org/10.1109/TMI.2006.876169">Riddell published it calling this "Rectification"</a>), I'm not sure that would change the computation time by a large factor but I should check. I think you then need an additional interpolation no to resample the "moved" object, no</i></b>”<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif"'>Not almost the same. What I mean seems to be simplier. Lets assume the detector has tilted around the x-axis by 0.25 degrees. The object’s center is at 0,0,0 and rotation axis is 0,1,0. In FDK you usually got a ray from source to the voxel center and then you calculate<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif"'>the intersection of this ray at the detector plane. To avoid ray plane intersection calculus, we can rotate the whole system by 0.25 degrees. Means that the detectors normal now is parallel to the z axis.Of course the rotation axis is not 0,1,0 any longer and the focus (xray source)<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif"'>also is a bit rotated. <span style='color:black'><o:p></o:p></span></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'><o:p> </o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif"'>“<b><i>For the backprojection, we typically used voxel-based backprojection using the center of the voxel which is faster than what you (seem to) use</i></b>”<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif"'>I think here we mixed something up. This type of back projection is used for FDK. All my questions were related to S-ART. I need to calculate the weights<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif"'>of a voxel for the back projection. To speed it up, I project its vertices on the detector plane, calc the MEB and the send rays from within this MEB through the<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif"'>voxel and caluclate the intersection length so those rays that will make my weights.<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif"'><o:p> </o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif"'><o:p> </o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif"'>Best regards,<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif"'>Robert<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'><o:p> </o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'><o:p> </o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'>P.S. Hello to all the other people here in the mailing list.<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:black'><o:p> </o:p></span></p><p class=MsoNormal><span lang=EN-US style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><b><span lang=EN-US style='font-size:10.0pt;font-family:"Tahoma","sans-serif"'>Von:</span></b><span lang=EN-US style='font-size:10.0pt;font-family:"Tahoma","sans-serif"'> <a href="mailto:simon.rit@gmail.com">simon.rit@gmail.com</a> [<a href="mailto:simon.rit@gmail.com">mailto:simon.rit@gmail.com</a>] <b>Im Auftrag von </b>Simon Rit<br><b>Gesendet:</b> Donnerstag, 29. Januar 2015 20:58<br><b>An:</b> Robert Calließ<br><b>Cc:</b> <a href="mailto:rtk-users@openrtk.org">rtk-users@openrtk.org</a><br><b>Betreff:</b> Re: RTK<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US><o:p> </o:p></span></p><p class=MsoNormal><span lang=EN-US>Hi,<br>Thank you for your interest in RTK. Please use the mailing list for questions that are of interest to anyone using RTK.<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US>There are many ways to model the direct problem (forward projection). Without going into too many details (available in the publications of each method) :<br>- "As far as I understand the goal of this approach is to calculate the intersection length of a ray through a voxel, right ?" False. Joseph's method samples the ray with one pixel per slice in the main direction but it does not compute the intersection of the ray with each voxel. <a href="http://www.ncbi.nlm.nih.gov/pubmed/4000088">Siddon's method</a> does that. In fig 2 of <a href="http://www3.cs.stonybrook.edu/~mueller/papers/ISBI_06_quality_2.pdf">[Xu and Mueller, 2006]</a>, Joseph is referred to "slice interpolated" and Siddon to "box-line-integrated".<br>- "I can calculate the intersection length of the ray within a voxel by a simple substraction, this runs very fast." This sounds very interesting, don't hesitate to share the code and/or the publication! BTW, what is DDA?<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US>- Small voxels / pixels, "Did you find a way to handle this ? " We don't handle this in RTK except if you consider that spatial regularisation (e.g., total variation) will overcome this problem in a way. But generally we use matching resolution between pixels and voxels so the problem is minimal. For the backprojection, we typically used voxel-based backprojection using the center of the voxel which is faster than what you (seem to) use. I think that if these things are a problem for you, there is a nice solution called <a href="http://stacks.iop.org/PMB/49/2463">distance driven (back-)projection</a> (by De Man and Basu). I think it will do exactly what you want. I haven't implemented it in RTK (yet).<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US>Thanks for the last trick, I am aware of it (<a href="http://dx.doi.org/10.1109/TMI.2006.876169">Riddell published it calling this "Rectification"</a>), I'm not sure that would change the computation time by a large factor but I should check. I think you then need an additional interpolation no to resample the "moved" object, no?<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US>I hope this helps. Let me know if sg is not clear in my answer!<br>Cheers,<br>Simon<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif"'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif"'><o:p> </o:p></span></p></div></body></html>