diff --git a/docs/build/doctrees/environment.pickle b/docs/build/doctrees/environment.pickle
index 901c5f54cb9fcb8ef376c25a195f992e0557a0fd..22575387852092bebd18f307e320b0c446a4b69b 100644
Binary files a/docs/build/doctrees/environment.pickle and b/docs/build/doctrees/environment.pickle differ
diff --git a/docs/build/doctrees/matlab.doctree b/docs/build/doctrees/matlab.doctree
index 4eca9da13b4350823bc5e801a9f63c2dd5e43d66..3d8afac69bf3baff835a0bc02a373bc804ff66fb 100644
Binary files a/docs/build/doctrees/matlab.doctree and b/docs/build/doctrees/matlab.doctree differ
diff --git a/docs/build/html/_sources/matlab.rst.txt b/docs/build/html/_sources/matlab.rst.txt
index d9789a905d136a4fa80f409439dc31598e35eb6e..06f25958c307168a624a0d0c7cbc9b14f1d153d0 100644
--- a/docs/build/html/_sources/matlab.rst.txt
+++ b/docs/build/html/_sources/matlab.rst.txt
@@ -7,3 +7,13 @@ Phase Retrieval
.. automodule:: functions.phaseretrieval
.. autofunction:: phaserec_ctf
+.. autofunction:: phaserec_mba
+.. autofunction:: phaserec_smo
+
+##########
+Tomography
+##########
+
+.. automodule:: functions.tomography
+
+.. autofunction:: ringremove
diff --git a/docs/build/html/genindex.html b/docs/build/html/genindex.html
index 05c2c874512f6b1fed3bb1b82bb414c9ae4f35d7..af569caceb000a3ef0a2b25a3615b491c17fdebe 100644
--- a/docs/build/html/genindex.html
+++ b/docs/build/html/genindex.html
@@ -39,12 +39,17 @@
<div class="genindex-jumpbox">
<a href="#F"><strong>F</strong></a>
| <a href="#P"><strong>P</strong></a>
+ | <a href="#R"><strong>R</strong></a>
</div>
<h2 id="F">F</h2>
<table style="width: 100%" class="indextable genindextable"><tr>
<td style="width: 33%; vertical-align: top;"><ul>
<li><a href="matlab.html#module-functions.phaseretrieval">functions.phaseretrieval (module)</a>
+</li>
+ </ul></td>
+ <td style="width: 33%; vertical-align: top;"><ul>
+ <li><a href="matlab.html#module-functions.tomography">functions.tomography (module)</a>
</li>
</ul></td>
</tr></table>
@@ -53,6 +58,20 @@
<table style="width: 100%" class="indextable genindextable"><tr>
<td style="width: 33%; vertical-align: top;"><ul>
<li><a href="matlab.html#functions.phaseretrieval.phaserec_ctf">phaserec_ctf() (in module functions.phaseretrieval)</a>
+</li>
+ </ul></td>
+ <td style="width: 33%; vertical-align: top;"><ul>
+ <li><a href="matlab.html#functions.phaseretrieval.phaserec_mba">phaserec_mba() (in module functions.phaseretrieval)</a>
+</li>
+ <li><a href="matlab.html#functions.phaseretrieval.phaserec_smo">phaserec_smo() (in module functions.phaseretrieval)</a>
+</li>
+ </ul></td>
+</tr></table>
+
+<h2 id="R">R</h2>
+<table style="width: 100%" class="indextable genindextable"><tr>
+ <td style="width: 33%; vertical-align: top;"><ul>
+ <li><a href="matlab.html#functions.tomography.ringremove">ringremove() (in module functions.tomography)</a>
</li>
</ul></td>
</tr></table>
@@ -78,6 +97,7 @@
<p class="caption"><span class="caption-text">Dies ist ein Test</span></p>
<ul>
<li class="toctree-l1"><a class="reference internal" href="matlab.html">Phase Retrieval</a></li>
+<li class="toctree-l1"><a class="reference internal" href="matlab.html#module-functions.tomography">Tomography</a></li>
</ul>
<div class="relations">
diff --git a/docs/build/html/index.html b/docs/build/html/index.html
index 6d4df74211c1a779d0e90b48cc325614961cd483..b75a2321eae9db11c62cc2ef381cba3d501fccb5 100644
--- a/docs/build/html/index.html
+++ b/docs/build/html/index.html
@@ -41,6 +41,7 @@
<p class="caption"><span class="caption-text">Dies ist ein Test</span></p>
<ul>
<li class="toctree-l1"><a class="reference internal" href="matlab.html">Phase Retrieval</a></li>
+<li class="toctree-l1"><a class="reference internal" href="matlab.html#module-functions.tomography">Tomography</a></li>
</ul>
</div>
</div>
@@ -73,6 +74,7 @@
<p class="caption"><span class="caption-text">Dies ist ein Test</span></p>
<ul>
<li class="toctree-l1"><a class="reference internal" href="matlab.html">Phase Retrieval</a></li>
+<li class="toctree-l1"><a class="reference internal" href="matlab.html#module-functions.tomography">Tomography</a></li>
</ul>
<div class="relations">
diff --git a/docs/build/html/mat-modindex.html b/docs/build/html/mat-modindex.html
index 2bebd56112b522662e9c3215687bc4dddcfa577d..69fc9ff1e1ee268a7bed6f3b99b8162e2dbddf4b 100644
--- a/docs/build/html/mat-modindex.html
+++ b/docs/build/html/mat-modindex.html
@@ -57,6 +57,11 @@
<td>   
<a href="matlab.html#module-functions.phaseretrieval"><code class="xref">functions.phaseretrieval</code></a></td><td>
<em></em></td></tr>
+ <tr class="cg-1">
+ <td></td>
+ <td>   
+ <a href="matlab.html#module-functions.tomography"><code class="xref">functions.tomography</code></a></td><td>
+ <em></em></td></tr>
</table>
@@ -79,6 +84,7 @@
<p class="caption"><span class="caption-text">Dies ist ein Test</span></p>
<ul>
<li class="toctree-l1"><a class="reference internal" href="matlab.html">Phase Retrieval</a></li>
+<li class="toctree-l1"><a class="reference internal" href="matlab.html#module-functions.tomography">Tomography</a></li>
</ul>
<div class="relations">
diff --git a/docs/build/html/matlab.html b/docs/build/html/matlab.html
index a17e6ea0dd0dd9ebc4f1d3bca1292e3a4aee7db2..c17efd266d190cba0a6b848c77d80f64989b102a 100644
--- a/docs/build/html/matlab.html
+++ b/docs/build/html/matlab.html
@@ -80,7 +80,7 @@ y-dimension before reconstruction</li>
</table>
<div class="admonition seealso">
<p class="first admonition-title">See also</p>
-<p class="last"><code class="xref mat mat-func docutils literal notranslate"><span class="pre">phaserec_mba()</span></code></p>
+<p class="last"><a class="reference internal" href="#functions.phaseretrieval.phaserec_mba" title="functions.phaseretrieval.phaserec_mba"><code class="xref mat mat-func docutils literal notranslate"><span class="pre">phaserec_mba()</span></code></a></p>
</div>
<p class="rubric">Notes</p>
<p>Example 1:</p>
@@ -116,7 +116,32 @@ y-dimension before reconstruction</li>
<span class="n">figure</span><span class="p">(</span><span class="s">'name'</span><span class="p">,</span> <span class="s">'Reconstructed phase-image'</span><span class="p">);</span> <span class="n">imagesc</span><span class="p">(</span><span class="n">phaseRec</span><span class="p">);</span> <span class="n">colorbar</span><span class="p">;</span> <span class="n">colormap</span> <span class="n">bone</span><span class="p">;</span>
</pre></div>
</div>
-<p>See also phaserec_mba</p>
+<p>See also PHASEREC_MBA</p>
+</dd></dl>
+
+<dl class="function">
+<dt id="functions.phaseretrieval.phaserec_mba">
+<code class="descclassname">functions.phaseretrieval.</code><code class="descname">phaserec_mba</code><span class="sig-paren">(</span><em>im1</em>, <em>fresnelNumber</em>, <em>settings</em><span class="sig-paren">)</span><a class="headerlink" href="#functions.phaseretrieval.phaserec_mba" title="Permalink to this definition">¶</a></dt>
+<dd><p>UNTITLED4 Summary of this function goes here
+Detailed explanation goes here</p>
+</dd></dl>
+
+<dl class="function">
+<dt id="functions.phaseretrieval.phaserec_smo">
+<code class="descclassname">functions.phaseretrieval.</code><code class="descname">phaserec_smo</code><span class="sig-paren">(</span><em>im1</em>, <em>fresnelNumber</em>, <em>settings</em><span class="sig-paren">)</span><a class="headerlink" href="#functions.phaseretrieval.phaserec_smo" title="Permalink to this definition">¶</a></dt>
+<dd><p>UNTITLED Summary of this function goes here
+Detailed explanation goes here</p>
+</dd></dl>
+
+</div>
+<div class="section" id="module-functions.tomography">
+<span id="tomography"></span><h1>Tomography<a class="headerlink" href="#module-functions.tomography" title="Permalink to this headline">¶</a></h1>
+<dl class="function">
+<dt id="functions.tomography.ringremove">
+<code class="descclassname">functions.tomography.</code><code class="descname">ringremove</code><span class="sig-paren">(</span><em>sino</em>, <em>settings</em><span class="sig-paren">)</span><a class="headerlink" href="#functions.tomography.ringremove" title="Permalink to this definition">¶</a></dt>
+<dd><p>this function can be used to remove ring artifacts in tomographic
+reconstructions by reducing the amount of horizontal lines in the given
+sinogram</p>
</dd></dl>
</div>
@@ -141,6 +166,7 @@ y-dimension before reconstruction</li>
<p class="caption"><span class="caption-text">Dies ist ein Test</span></p>
<ul class="current">
<li class="toctree-l1 current"><a class="current reference internal" href="#">Phase Retrieval</a></li>
+<li class="toctree-l1"><a class="reference internal" href="#module-functions.tomography">Tomography</a></li>
</ul>
<div class="relations">
diff --git a/docs/build/html/objects.inv b/docs/build/html/objects.inv
index 0d5fcde6d5ffebac066dd6ea3a1a949a24e3335b..0f89c5119b074248f56849c96b5a5bd0dbee6f74 100644
--- a/docs/build/html/objects.inv
+++ b/docs/build/html/objects.inv
@@ -2,5 +2,4 @@
# Project: HoloTomoToolbox
# Version:
# The remainder of this file is compressed using zlib.
-xÚM
-Â0…÷=ňn+vÛèB\JšL É”&•ºó^Ï“Ø4D[qÈÌ|oò^òVs[‘6˺d´M…&A1›*DH\!Y¶,’óÄQþ‡j5QùfüK J~æ6 6ò‹ƒu¥v`¬Hû JˆÝƒ»FoLwЈ8¡ä¤,Áš$ Iõ‡dFÝãv7Ð Z…Ú2gbyãU#W°sîa Ҍ闿iŠ‰/™êk þ¡ÿ_ŸÎ kx9Uøž‡~¾cFOSg°Ý
\ No newline at end of file
+xÚ¥MNÃ0…÷9Å ºMÔn{‚²@ªh¥.‘cO~$ÛÙ“ªÝq
®ÇIˆk™$¢‚ vñÌû^Þ¼ª·’[²¾èáÑ!»ÏBƒ¼5¤z° -Ê¢a£ WgÕj=£â0ÿ NOù*¹º9%å Œ&¦ÿ7ñ†–›0ªèšëÒöî šÁµ¶vh茿%ѶVá<«í°A
ùÒô¦ƒ§ð™2E2ˆ&ŠjI v¤é8$9é’. oï 7hY„YÌ0µš¤‚}è^RïSõW¾ùßRwnê®yZÄ =Çfãu …“Íœˆ³(>Äý^Ô˜}ço!1
\ No newline at end of file
diff --git a/docs/build/html/search.html b/docs/build/html/search.html
index 33fa615d880486f07387519db2783f2704938872..5f37cee7bc83288315ce26aca62afdfd1f8a903f 100644
--- a/docs/build/html/search.html
+++ b/docs/build/html/search.html
@@ -83,6 +83,7 @@
<p class="caption"><span class="caption-text">Dies ist ein Test</span></p>
<ul>
<li class="toctree-l1"><a class="reference internal" href="matlab.html">Phase Retrieval</a></li>
+<li class="toctree-l1"><a class="reference internal" href="matlab.html#module-functions.tomography">Tomography</a></li>
</ul>
<div class="relations">
diff --git a/docs/build/html/searchindex.js b/docs/build/html/searchindex.js
index 853063554ceeb991fe065e5cc0bc9eaa0e1f923e..e6aa983de458867b20cd118b2d06afaa2bc475a3 100644
--- a/docs/build/html/searchindex.js
+++ b/docs/build/html/searchindex.js
@@ -1 +1 @@
-Search.setIndex({docnames:["index","matlab"],envversion:53,filenames:["index.rst","matlab.rst"],objects:{"functions.phaseretrieval":{phaserec_ctf:[1,1,1,""]},functions:{phaseretrieval:[1,0,0,"-"]}},objnames:{"0":["mat","module","MATLAB module"],"1":["mat","function","MATLAB function"]},objtypes:{"0":"mat:module","1":"mat:function"},terms:{"default":1,"function":1,"return":1,"true":1,Dies:0,The:1,absorpt:1,acquir:1,addit:1,along:1,also:1,amount:1,approach:1,arrai:1,asdf:[],asdfasdfasdf:[],assign:1,assum:1,astigmat:1,autofunct:[],automodul:[],base:1,befor:1,betadeltaratio:1,between:1,bla:[],block:[],blubb:[],bone:1,cell:1,cloeten:1,code:[],colorbar:1,colormap:1,contain:1,correspond:1,ctf:1,deviat:1,dicti:1,dicty_sketch_2048x2048:1,dimens:1,dummi:[],ein:0,empti:1,exampl:1,fals:1,figur:1,fix:1,four:1,frequenc:1,fresn:1,fresnel:1,from:1,high:1,holo:1,hologram:1,homogen:1,imag:1,imagesc:1,index:0,input:1,invers:1,ist:0,lim1:1,lim2:1,load:1,low:1,mat:1,matlab:[],modul:0,name:1,note:1,number:1,numer:1,object:1,one:1,option:1,other:1,pad:1,padi:1,padx:1,page:0,paramet:1,peter:1,phantom:1,phase:0,phase_reconstruct:[],phaserec:1,phaserec_ctf:1,phaserec_mba:1,phaseretriev:1,phasetru:1,pixel:1,possibl:1,pure:1,rand:1,ratio:1,reconstruct:1,regularis:1,replic:1,result:1,retriev:0,same:1,search:0,see:1,set:1,sever:1,shift:1,simset:1,simulatehologram:1,simulatelinear:1,size:1,structur:1,symmetr:1,test:0,type:1,valu:1,vector:1,were:1,which:1},titles:["Welcome to HoloTomoToolbox\u2019s documentation!","Phase Retrieval"],titleterms:{"function":[],api:[],document:0,holotomotoolbox:0,indic:0,matlab:[],phase:1,retriev:1,tabl:0,welcom:0}})
\ No newline at end of file
+Search.setIndex({docnames:["index","matlab"],envversion:53,filenames:["index.rst","matlab.rst"],objects:{"functions.phaseretrieval":{phaserec_ctf:[1,1,1,""],phaserec_mba:[1,1,1,""],phaserec_smo:[1,1,1,""]},"functions.tomography":{ringremove:[1,1,1,""]},functions:{phaseretrieval:[1,0,0,"-"],tomography:[1,0,0,"-"]}},objnames:{"0":["mat","module","MATLAB module"],"1":["mat","function","MATLAB function"]},objtypes:{"0":"mat:module","1":"mat:function"},terms:{"default":1,"function":1,"return":1,"true":1,Dies:0,The:1,absorpt:1,acquir:1,addit:1,along:1,also:1,amount:1,approach:1,arrai:1,artifact:1,asdf:[],asdfasdfasdf:[],assign:1,assum:1,astigmat:1,autofunct:[],automodul:[],base:1,befor:1,betadeltaratio:1,between:1,bla:[],block:[],blubb:[],bone:1,can:1,cell:1,cloeten:1,code:[],colorbar:1,colormap:1,contain:1,correspond:1,ctf:1,detail:1,deviat:1,dicti:1,dicty_sketch_2048x2048:1,dimens:1,dummi:[],ein:0,empti:1,exampl:1,explan:1,fals:1,figur:1,fix:1,four:1,frequenc:1,fresn:1,fresnel:1,fresnelnumb:1,from:1,given:1,goe:1,here:1,high:1,holo:1,hologram:1,homogen:1,horizont:1,im1:1,imag:1,imagesc:1,index:0,input:1,invers:1,ist:0,lim1:1,lim2:1,line:1,load:1,low:1,mat:1,matlab:[],modul:0,name:1,note:1,number:1,numer:1,object:1,one:1,option:1,other:1,pad:1,padi:1,padx:1,page:0,paramet:1,peter:1,phantom:1,phase:0,phase_reconstruct:[],phaserec:1,phaserec_ctf:1,phaserec_mba:1,phaserec_smo:1,phaseretriev:1,phasetru:1,pixel:1,possibl:1,pure:1,rand:1,ratio:1,reconstruct:1,reduc:1,regularis:1,remov:1,replic:1,result:1,retriev:0,ring:1,ringremov:1,same:1,search:0,see:1,set:1,sever:1,shift:1,simset:1,simulatehologram:1,simulatelinear:1,sino:1,sinogram:1,size:1,structur:1,summari:1,symmetr:1,test:0,thi:1,tomograph:1,tomographi:0,type:1,untitl:1,untitled4:1,used:1,valu:1,vector:1,were:1,which:1},titles:["Welcome to HoloTomoToolbox\u2019s documentation!","Phase Retrieval"],titleterms:{"function":[],api:[],document:0,holotomotoolbox:0,indic:0,matlab:[],phase:1,retriev:1,tabl:0,tomographi:1,welcom:0}})
\ No newline at end of file
diff --git a/docs/source/matlab.rst b/docs/source/matlab.rst
index d9789a905d136a4fa80f409439dc31598e35eb6e..06f25958c307168a624a0d0c7cbc9b14f1d153d0 100644
--- a/docs/source/matlab.rst
+++ b/docs/source/matlab.rst
@@ -7,3 +7,13 @@ Phase Retrieval
.. automodule:: functions.phaseretrieval
.. autofunction:: phaserec_ctf
+.. autofunction:: phaserec_mba
+.. autofunction:: phaserec_smo
+
+##########
+Tomography
+##########
+
+.. automodule:: functions.tomography
+
+.. autofunction:: ringremove
diff --git a/functions/tomography/alignProfiles.m b/functions/tomography/alignProfiles.m
index a4850e8a7b59e8556ec44b4c2b27098d751f32c1..c8fd058dcc279a930f65fee743c3bc03568a592e 100755
--- a/functions/tomography/alignProfiles.m
+++ b/functions/tomography/alignProfiles.m
@@ -1,73 +1,73 @@
-function shiftval = alignProfiles(profToAlign,profReference,settings)
-% Aligns two 1D signals and returns the shift between them
-
-% profToAlign will be shifted with stepsize (subpixel) between -range and
-% +range. The squared difference or standard deviation is calculated for each shift and the position of
-% minimal error is returned.
-%
-% Calculation is performed recursively with decreasing stepsize and
-% range until the stepsize is smaller than the accuracy.
-% If the inital range is to small it will be enlarged.
-
-
-if nargin<3
- settings=struct;
-end
-
-defaults.method = 'diff';
-
-defaults.accuracy=0.01;
-defaults.range=20;
-defaults.stepsize=0;
-defaults.offset=0;
-
-if (nargin == 0)
- shiftval=defaults;
- return
-end
-
-settings=completeStruct(settings, defaults);
-
-% dependent settings
-if settings.stepsize == 0
- settings.stepsize = settings.range/20;
-end
-
-
-% shifts to try out
-shifts = -settings.range+settings.offset:settings.stepsize:settings.range+settings.offset;
-
-err = zeros(numel(shifts),1);
-for iShift=1:numel(shifts)
- profShifted = circshiftSubPixel(profToAlign,[shifts(iShift) 0]);
- switch settings.method
- case 'diff'
- err(iShift) = sum((profReference(ceil(settings.range)+1:end-ceil(settings.range))-profShifted(ceil(settings.range)+1:end-ceil(settings.range))).^2);
- case 'std'
- err(iShift) = sum(std([profReference(ceil(settings.range)+1:end-ceil(settings.range)) profShifted(ceil(settings.range)+1:end-ceil(settings.range))],0,2));
- end
-end
-
-% if multiple minima exist, take the mean value of them
-shiftvals = shifts(err==min(err));
-shiftval = mean(shiftvals);
-
-if sum(abs(shiftvals) == abs(settings.range))>=1
- % if minimum value is on the edge of the range, enlarge the range by a factor of 2 to
- % ensure that it is a real minimum
- disp (['enlarging range to ', num2str(settings.range*2)])
- settings.range=settings.range*2;
- shiftval = alignProfiles(profToAlign,profReference,settings);
-else
- % value was not on the edge, increase accuracy by decreasing
- % stepsize by a factor of 5; use current guess as offset
- if settings.stepsize>settings.accuracy
- settings.range=settings.range/5;
- settings.stepsize=settings.stepsize/5;
- settings.offset=shiftval;
- shiftval = alignProfiles(profToAlign,profReference,settings);
- end
-end
-
-end
-
+function shiftval = alignProfiles(profToAlign,profReference,settings)
+% Aligns two 1D signals and returns the shift between them
+
+% profToAlign will be shifted with stepsize (subpixel) between -range and
+% +range. The squared difference or standard deviation is calculated for each shift and the position of
+% minimal error is returned.
+%
+% Calculation is performed recursively with decreasing stepsize and
+% range until the stepsize is smaller than the accuracy.
+% If the inital range is to small it will be enlarged.
+
+
+if nargin<3
+ settings=struct;
+end
+
+defaults.method = 'diff';
+
+defaults.accuracy=0.01;
+defaults.range=20;
+defaults.stepsize=0;
+defaults.offset=0;
+
+if (nargin == 0)
+ shiftval=defaults;
+ return
+end
+
+settings=completeStruct(settings, defaults);
+
+% dependent settings
+if settings.stepsize == 0
+ settings.stepsize = settings.range/20;
+end
+
+
+% shifts to try out
+shifts = -settings.range+settings.offset:settings.stepsize:settings.range+settings.offset;
+
+err = zeros(numel(shifts),1);
+for iShift=1:numel(shifts)
+ profShifted = circshiftSubPixel(profToAlign,[shifts(iShift) 0]);
+ switch settings.method
+ case 'diff'
+ err(iShift) = sum((profReference(ceil(settings.range)+1:end-ceil(settings.range))-profShifted(ceil(settings.range)+1:end-ceil(settings.range))).^2);
+ case 'std'
+ err(iShift) = sum(std([profReference(ceil(settings.range)+1:end-ceil(settings.range)) profShifted(ceil(settings.range)+1:end-ceil(settings.range))],0,2));
+ end
+end
+
+% if multiple minima exist, take the mean value of them
+shiftvals = shifts(err==min(err));
+shiftval = mean(shiftvals);
+
+if sum(abs(shiftvals) == abs(settings.range))>=1
+ % if minimum value is on the edge of the range, enlarge the range by a factor of 2 to
+ % ensure that it is a real minimum
+ disp (['enlarging range to ', num2str(settings.range*2)])
+ settings.range=settings.range*2;
+ shiftval = alignProfiles(profToAlign,profReference,settings);
+else
+ % value was not on the edge, increase accuracy by decreasing
+ % stepsize by a factor of 5; use current guess as offset
+ if settings.stepsize>settings.accuracy
+ settings.range=settings.range/5;
+ settings.stepsize=settings.stepsize/5;
+ settings.offset=shiftval;
+ shiftval = alignProfiles(profToAlign,profReference,settings);
+ end
+end
+
+end
+
diff --git a/functions/tomography/astraConeProjection.m b/functions/tomography/astraConeProjection.m
index 3cdcbe3e32a7d13c14e97c144fde2f59f9db5ae7..7c567d9937fc2e2ca30c4a38a71a2aca6b3f90bf 100755
--- a/functions/tomography/astraConeProjection.m
+++ b/functions/tomography/astraConeProjection.m
@@ -1,65 +1,65 @@
-function projs = astraConeProjection(vol,tomoAngles,z01,z02,dx,dx_eff,settings)
-% uses the Astra Toolbox to create projections of a given volume in
-% cone-beam geometry
-%
-% Important: z01, z02, dx and dx_eff must have the same dimensions (e.g.
-% all mm or all m)
-
-%
-% Optional settings to be passed as 7th argument (and defaults)
-%
-% outputSize - size of the projections (default: size of projections of Matlab's radon)
-% radonOrientation - should orientation of Matlab's radon and ASTRA
-% coincide? (1=true)
-%
-% author: Martin Krenkel and Mareike Toepperwien
-% credits go to programmers of the Astra Toolbox
-
- if (nargin < 7)
- settings = struct;
- end
-
- % defaults
- defaults.outputSize = 0; % size of the projections, the default of 0 results in the same size as the output of Matlab's radon transform
- defaults.radonOrientation = 1; % returns projections with the same orientation as Matlab's radon transform
-
- % return default settings
- if (nargin == 0)
- projs=defaults;
- return
- end
-
- settings = completeStruct(settings, defaults);
-
- % dependent settings
- if settings.outputSize == 0
- % as in radon
- sizeProjs = 2*ceil(norm([size(vol,1) size(vol,2)]-floor(([size(vol,1) size(vol,2)]-1)/2)-1))+3;
- settings.outputSize = [size(vol,1) sizeProjs];
- end
-
- M = z02/z01;
-
- if (nargin < 6)
- dx_eff = dx/M;
- end
-
-%% here it starts
-if settings.radonOrientation==1
- vol = permute(vol,[2 1 3]);
-end
-% volume geometry
-volGeom = astra_create_vol_geom(size(vol,2),size(vol,1),size(vol,3));
-
-% Create a data object for the projection
-projGeom = astra_create_proj_geom('cone', M, M, settings.outputSize(1), settings.outputSize(2), -tomoAngles/180*pi, z01/dx_eff, (z02-z01)/dx_eff);
-
-[projID, projs] = astra_create_sino3d_cuda(vol, projGeom, volGeom);
-
-% Clean up
-astra_mex_data3d('delete', projID);
-
-projs = permute(projs,[3 1 2]);
-
-end
-
+function projs = astraConeProjection(vol,tomoAngles,z01,z02,dx,dx_eff,settings)
+% uses the Astra Toolbox to create projections of a given volume in
+% cone-beam geometry
+%
+% Important: z01, z02, dx and dx_eff must have the same dimensions (e.g.
+% all mm or all m)
+
+%
+% Optional settings to be passed as 7th argument (and defaults)
+%
+% outputSize - size of the projections (default: size of projections of Matlab's radon)
+% radonOrientation - should orientation of Matlab's radon and ASTRA
+% coincide? (1=true)
+%
+% author: Martin Krenkel and Mareike Toepperwien
+% credits go to programmers of the Astra Toolbox
+
+ if (nargin < 7)
+ settings = struct;
+ end
+
+ % defaults
+ defaults.outputSize = 0; % size of the projections, the default of 0 results in the same size as the output of Matlab's radon transform
+ defaults.radonOrientation = 1; % returns projections with the same orientation as Matlab's radon transform
+
+ % return default settings
+ if (nargin == 0)
+ projs=defaults;
+ return
+ end
+
+ settings = completeStruct(settings, defaults);
+
+ % dependent settings
+ if settings.outputSize == 0
+ % as in radon
+ sizeProjs = 2*ceil(norm([size(vol,1) size(vol,2)]-floor(([size(vol,1) size(vol,2)]-1)/2)-1))+3;
+ settings.outputSize = [size(vol,1) sizeProjs];
+ end
+
+ M = z02/z01;
+
+ if (nargin < 6)
+ dx_eff = dx/M;
+ end
+
+%% here it starts
+if settings.radonOrientation==1
+ vol = permute(vol,[2 1 3]);
+end
+% volume geometry
+volGeom = astra_create_vol_geom(size(vol,2),size(vol,1),size(vol,3));
+
+% Create a data object for the projection
+projGeom = astra_create_proj_geom('cone', M, M, settings.outputSize(1), settings.outputSize(2), -tomoAngles/180*pi, z01/dx_eff, (z02-z01)/dx_eff);
+
+[projID, projs] = astra_create_sino3d_cuda(vol, projGeom, volGeom);
+
+% Clean up
+astra_mex_data3d('delete', projID);
+
+projs = permute(projs,[3 1 2]);
+
+end
+
diff --git a/functions/tomography/astraFDK.m b/functions/tomography/astraFDK.m
index d5ee43703f3bb2f413ac414c9b9a7ac5d905a611..793c1a5332531411716c471a2f175a08f8010e5d 100755
--- a/functions/tomography/astraFDK.m
+++ b/functions/tomography/astraFDK.m
@@ -1,95 +1,95 @@
-function [vol] = astraFDK( projs,tomoAngles,z01,z02,dx,dx_eff,settings)
-% uses the Astra Toolbox to perform FDK reconstruction
-%
-% usage: vol = astraFDK( projs,tomoAngles,z01,z02,dx,dx_eff,settings)
-%
-% Important: z01, z02, dx and dx_eff must have the same dimensions (e.g.
-% all mm or all m)
-%
-% If the reconstruction anyway looks weird try to change the sign of
-% tomoAngles.
-%
-% Optional settings to be passed as 7th argument (and defaults)
-%
-% outputSize - size of reconstructed slices (0 = automatic)
-% shortScan - scan was not from 0 to 360° (0=false)
-% numSlices - number of slices to be reconstructed (1)
-% offset - offset to central slice of detector (0)
-% iradonOrientation - should orientation of Matlab's iradon and ASTRA
-% coincide? (1=true)
-%
-% authors: Martin Krenkel and Mareike Töpperwien
-% credits go to programmers of the Astra Toolbox
-
- if (nargin < 7)
- settings = struct;
- end
-
- % defaults
- defaults.outputSize = 0;
- defaults.shortScan = 1;
- defaults.numSlices = 1;
- defaults.offset=0;
- defaults.iradonOrientation = 1;
-
- % return default settings
- if (nargin == 0)
- vol=defaults;
- return
- end
-
- settings = completeStruct(settings, defaults);
-
-
- % dependent settings
- if settings.outputSize == 0
- % as in iradon
- settings.outputSize = 2*floor(numel(tomoAngles)/(2*sqrt(2)));
- end
-
- %% check usage
- if numel(tomoAngles) ~= size(projs,3)
- error('Number of tomographic angles does not match number of projections');
- end
-
-%% here it starts
-
-% volume geometry
-% sizeY, sizeX, sizeZ, minX, maxX, minY, maxY, minZ, maxZ
-
-volGeom = astra_create_vol_geom(settings.outputSize, settings.outputSize, settings.numSlices, -settings.outputSize/2, settings.outputSize/2, -settings.outputSize/2, settings.outputSize/2, settings.offset-floor(settings.numSlices/2), settings.offset+ceil(settings.numSlices/2));
-
-% create a projection geometry and load data into memory
-M = dx/dx_eff;
-
-% keyboard;
-projGeom = astra_create_proj_geom('cone', M, M, size(projs,1), size(projs,2), -tomoAngles/180*pi, z01/dx_eff, (z02-z01)/dx_eff);
-projID = astra_mex_data3d('create', '-proj3d', projGeom, permute(projs,[2 3 1]));
-% Create a data object for the reconstruction
-recID = astra_mex_data3d('create', '-vol', volGeom);
-
-% Set up the parameters for a reconstruction algorithm using the GPU
-cfg = astra_struct('FDK_CUDA');
-cfg.ReconstructionDataId = recID;
-cfg.ProjectionDataId = projID;
-cfg.option.ShortScan = settings.shortScan;
-
-% Create the algorithm object from the configuration structure
-algID = astra_mex_algorithm('create', cfg);
-% astra_mex_algorithm('run', algID, 1);
-astra_mex_algorithm('iterate', algID, 1);
-% Get the result
-vol = astra_mex_data3d('get', recID);
-
-% same orientation as iradon
-if settings.iradonOrientation==1
- vol=permute(vol,[2 1 3]);
-end
-
-% Clean up
-astra_mex_algorithm('delete', algID);
-astra_mex_data3d('delete', recID);
-astra_mex_data3d('delete', projID);
-
-end
-
+function [vol] = astraFDK( projs,tomoAngles,z01,z02,dx,dx_eff,settings)
+% uses the Astra Toolbox to perform FDK reconstruction
+%
+% usage: vol = astraFDK( projs,tomoAngles,z01,z02,dx,dx_eff,settings)
+%
+% Important: z01, z02, dx and dx_eff must have the same dimensions (e.g.
+% all mm or all m)
+%
+% If the reconstruction anyway looks weird try to change the sign of
+% tomoAngles.
+%
+% Optional settings to be passed as 7th argument (and defaults)
+%
+% outputSize - size of reconstructed slices (0 = automatic)
+% shortScan - scan was not from 0 to 360° (0=false)
+% numSlices - number of slices to be reconstructed (1)
+% offset - offset to central slice of detector (0)
+% iradonOrientation - should orientation of Matlab's iradon and ASTRA
+% coincide? (1=true)
+%
+% authors: Martin Krenkel and Mareike Töpperwien
+% credits go to programmers of the Astra Toolbox
+
+ if (nargin < 7)
+ settings = struct;
+ end
+
+ % defaults
+ defaults.outputSize = 0;
+ defaults.shortScan = 1;
+ defaults.numSlices = 1;
+ defaults.offset=0;
+ defaults.iradonOrientation = 1;
+
+ % return default settings
+ if (nargin == 0)
+ vol=defaults;
+ return
+ end
+
+ settings = completeStruct(settings, defaults);
+
+
+ % dependent settings
+ if settings.outputSize == 0
+ % as in iradon
+ settings.outputSize = 2*floor(numel(tomoAngles)/(2*sqrt(2)));
+ end
+
+ %% check usage
+ if numel(tomoAngles) ~= size(projs,3)
+ error('Number of tomographic angles does not match number of projections');
+ end
+
+%% here it starts
+
+% volume geometry
+% sizeY, sizeX, sizeZ, minX, maxX, minY, maxY, minZ, maxZ
+
+volGeom = astra_create_vol_geom(settings.outputSize, settings.outputSize, settings.numSlices, -settings.outputSize/2, settings.outputSize/2, -settings.outputSize/2, settings.outputSize/2, settings.offset-floor(settings.numSlices/2), settings.offset+ceil(settings.numSlices/2));
+
+% create a projection geometry and load data into memory
+M = dx/dx_eff;
+
+% keyboard;
+projGeom = astra_create_proj_geom('cone', M, M, size(projs,1), size(projs,2), -tomoAngles/180*pi, z01/dx_eff, (z02-z01)/dx_eff);
+projID = astra_mex_data3d('create', '-proj3d', projGeom, permute(projs,[2 3 1]));
+% Create a data object for the reconstruction
+recID = astra_mex_data3d('create', '-vol', volGeom);
+
+% Set up the parameters for a reconstruction algorithm using the GPU
+cfg = astra_struct('FDK_CUDA');
+cfg.ReconstructionDataId = recID;
+cfg.ProjectionDataId = projID;
+cfg.option.ShortScan = settings.shortScan;
+
+% Create the algorithm object from the configuration structure
+algID = astra_mex_algorithm('create', cfg);
+% astra_mex_algorithm('run', algID, 1);
+astra_mex_algorithm('iterate', algID, 1);
+% Get the result
+vol = astra_mex_data3d('get', recID);
+
+% same orientation as iradon
+if settings.iradonOrientation==1
+ vol=permute(vol,[2 1 3]);
+end
+
+% Clean up
+astra_mex_algorithm('delete', algID);
+astra_mex_data3d('delete', recID);
+astra_mex_data3d('delete', projID);
+
+end
+
diff --git a/functions/tomography/checkRotaxis.m b/functions/tomography/checkRotaxis.m
index 91dfbba4b997a3a3984680e3c71d0772a911c5f8..8fe0ab4e1409ef57ec58c13a05345cd66459a386 100755
--- a/functions/tomography/checkRotaxis.m
+++ b/functions/tomography/checkRotaxis.m
@@ -1,55 +1,56 @@
-function result = checkRotaxis(proj0deg,proj180deg,settings)
-
-% Compares two images which were recorded at an angular difference of 180°; one of the images is
-% flipped around the vertical axis; if the rotation axis is in the center
-% of the detector, both images should cancel each other out; in case of a
-% misaligned rotation axis, the images have to be shifted against each other
-% by the distance of the actual rotation axis to the center of the detector
-
-% Optional settings to be passed as 3rd argument
-%
-% shiftx - horizontal shift needed so that the images cancel each
-% other out
-% shifty - vertical shift needed so that the images cancel each
-% other out (this should be zero as otherwise the sample moved in vertical direction during the acquisition)
-% tiltAngle - tilt angle of the detector with respect to the rotation
-% axis
-% cutoff - cuts away the specified number of pixels from the edges of
-% the images
-
-if nargin<3
- settings=struct;
-end
-
-defaults.shiftx = 0;
-defaults.shifty = 0;
-defaults.tiltAngle = 0;
-defaults.cutoff = 0;
-
-if (nargin == 0)
- result = defaults;
- return
-end
-
-settings=completeStruct(settings, defaults);
-
-% rotate the projections
-if settings.tiltAngle~=0
- proj0deg = imrotate(proj0deg,settings.tiltAngle,'bilinear','crop');
- proj180deg = imrotate(proj180deg,settings.tiltAngle,'bilinear','crop');
-end
-
-% shift the projections
-if mod(settings.shiftx,1)==0 || mod(settings.shifty,1)==0
- proj0deg = circshiftSubPixel(proj0deg,[0 -settings.shiftx]);
- proj180deg = circshiftSubPixel(proj180deg,[settings.shifty -settings.shiftx]);
-else
- proj0deg = circshiftSubPixel(proj0deg,[0 -settings.shiftx]);
- proj180deg = circshiftSubPixel(proj180deg,[settings.shifty -settings.shiftx]);
-end
-
-% show the difference between the first projection and the second
-% projection (vertically flipped)
-result = proj0deg(settings.cutoff+1:end-settings.cutoff,settings.cutoff+1:end-settings.cutoff)-fliplr(proj180deg(settings.cutoff+1:end-settings.cutoff,settings.cutoff+1:end-settings.cutoff));
-
-end
\ No newline at end of file
+function result = checkRotaxis(proj0deg,proj180deg,settings)
+
+% Compares two images which were recorded at an angular difference of 180°; one of the images is
+% flipped around the vertical axis; if the rotation axis is in the center
+% of the detector, both images should cancel each other out; in case of a
+% misaligned rotation axis, the images have to be shifted against each other
+% by the distance of the actual rotation axis to the center of the detector
+
+% Optional settings to be passed as 3rd argument
+%
+% shiftx - horizontal shift needed so that the images cancel each
+% other out
+% shifty - vertical shift needed so that the images cancel each
+% other out (this should be zero as otherwise the sample moved in vertical direction during the acquisition)
+% tiltAngle - tilt angle of the detector with respect to the rotation
+% axis
+% cutoff - cuts away the specified number of pixels from the edges of
+% the images
+
+if nargin<3
+ settings=struct;
+end
+
+defaults.shiftx = 0;
+defaults.shifty = 0;
+defaults.tiltAngle = 0;
+defaults.cutoff = 0;
+
+if (nargin == 0)
+ result = defaults;
+ return
+end
+
+settings=completeStruct(settings, defaults);
+
+% rotate the projections
+if settings.tiltAngle~=0
+ proj0deg = imrotate(proj0deg,settings.tiltAngle,'bilinear','crop');
+ proj180deg = imrotate(proj180deg,settings.tiltAngle,'bilinear','crop');
+end
+
+% shift the projections
+if mod(settings.shiftx,1)==0 || mod(settings.shifty,1)==0
+ proj0deg = circshiftSubPixel(proj0deg,[0 -settings.shiftx]);
+ proj180deg = circshiftSubPixel(proj180deg,[settings.shifty -settings.shiftx]);
+else
+ proj0deg = circshiftSubPixel(proj0deg,[0 -settings.shiftx]);
+ proj180deg = circshiftSubPixel(proj180deg,[settings.shifty -settings.shiftx]);
+end
+
+% show the difference between the first projection and the second
+% projection (vertically flipped)
+result = proj0deg(settings.cutoff+1:end-settings.cutoff,settings.cutoff+1:end-settings.cutoff)-fliplr(proj180deg(settings.cutoff+1:end-settings.cutoff,settings.cutoff+1:end-settings.cutoff));
+
+end
+
diff --git a/functions/tomography/create3Dphantom.m b/functions/tomography/create3Dphantom.m
index a398603427d3c673ef9ba3bc779bd1708f49f024..a98b206ea841e362575720d506f6311283d1b953 100755
--- a/functions/tomography/create3Dphantom.m
+++ b/functions/tomography/create3Dphantom.m
@@ -1,102 +1,102 @@
-function phantom = create3Dphantom(outputSize,settings)
-
-% this function creates a 3D phantom in which a specified number of objects
-% with a specified size and geometry are positioned; the positions of these
-% objects can be either predefined or random and also the additive strength
-% of each object can be chosen as a fixed or random value
-
-%{
-% Example: create a phantom with spheres of random size and strength
-settings=create3Dphantom;
-settings.geometry='sphere';
-settings.numObjects=30;
-settings.objectStrength=1;
-settings.distrStrength='random';
-settings.sizeObjects=100;
-settings.distrSize='random';
-phantom=create3Dphantom(512,settings);
-
-figure; imagesc(phantom(:,:,randi(size(phantom,3)))); colorbar;
-
-%}
-
-
-if (nargin < 2)
- settings = struct;
-end
-
-% defaults
-defaults.sizeObjects = 0; % size of the objects; 0 leads to a size of 1/10 of the outputSize
-defaults.distrSize = 'fixed'; % switch between 'fixed' for fixed size or 'random' with a random size distribution with maximum settings.sizeObjects and minimum (settings.sizeObjects)/2
-defaults.numObjects = 10; % number of objects within the 3D volume
-defaults.geometry = 'sphere'; % geometry of the object; choices: sphere, cube and cylinder (with a random orientation along the three main axes)
-defaults.positions = 0; % positions of the object in space; 0 leads to a random distribution
-defaults.objectStrength = 1; % additive strength of the objects;
-defaults.distrStrength = 'fixed'; % switch between 'fixed' for fixed object strength or 'random' with a random strength between settings.objectStrength and (settings.objectStrength)/2 for each object
-
-% return default settings
-if (nargin == 0)
- phantom=defaults;
- return
-end
-
-settings = completeStruct(settings, defaults);
-
-% set size of the objects to 1/10 of the output size if no other settings
-% were specified
-if settings.sizeObjects == 0
- settings.sizeObjects = floor(outputSize/10);
-end
-
-% determine the radius for each of the objects (either random or fixed value)
-if strcmp(settings.distrSize,'random')
- radius=(settings.sizeObjects-rand(settings.numObjects,1)*0.5*settings.sizeObjects)/2;
-else
- radius=(repmat(settings.sizeObjects,[settings.numObjects 1]))/2;
-end
-
-
-% determine the additive strength for each of the objects (either random or fixed value)
-if strcmp(settings.distrStrength,'random')
- delta=settings.objectStrength-rand(settings.numObjects,1)*0.5*settings.objectStrength;
-else
- delta=repmat(settings.objectStrength,[settings.numObjects 1]);
-end
-
-% determine random positions of objects within the 3D array if no positions
-% were given
-if settings.positions == 0
- settings.positions = floor((rand(settings.numObjects,3))*(outputSize-2*max(radius)))+max(radius);
-end
-
-% create 3D object
-phantom = zeros(outputSize,outputSize,outputSize,'single');
-[X,Y,Z]=meshgrid(1:outputSize);
-
-% place the individual objects on their positions in the 3D volume
-for objectIdx=1:settings.numObjects
- % position of the center of the object
- posX = settings.positions(objectIdx,1);
- posY = settings.positions(objectIdx,2);
- posZ = settings.positions(objectIdx,3);
-
- switch settings.geometry
- case 'sphere'
- phantom(((X-posX).^2+(Y-posY).^2+(Z-posZ).^2)<=radius(objectIdx)^2)=phantom(((X-posX).^2+(Y-posY).^2+(Z-posZ).^2)<=radius(objectIdx)^2)+delta(objectIdx);
- case 'cube'
- phantom((abs(X-posX)<radius(objectIdx))&(abs(Y-posY)<radius(objectIdx))&(abs(Z-posZ)<radius(objectIdx)))=phantom((abs(X-posX)<radius(objectIdx))&(abs(Y-posY)<radius(objectIdx))&(abs(Z-posZ)<radius(objectIdx)))+delta(objectIdx);
- case 'cylinder'
- orientation = rand*3; % determine orientation of the cylinder
- if orientation < 1
- phantom((((Y-posY).^2+(Z-posZ).^2)<=radius(objectIdx)^2)&(abs(X-posX)<radius(objectIdx)))=phantom((((Y-posY).^2+(Z-posZ).^2)<=radius(objectIdx)^2)&(abs(X-posX)<radius(objectIdx)))+delta(objectIdx);
- elseif orientation < 2
- phantom((((X-posX).^2+(Z-posZ).^2)<=radius(objectIdx)^2)&(abs(Y-posY)<radius(objectIdx)))=phantom((((X-posX).^2+(Z-posZ).^2)<=radius(objectIdx)^2)&(abs(Y-posY)<radius(objectIdx)))+delta(objectIdx);
- else
- phantom((((Y-posY).^2+(X-posX).^2)<=radius(objectIdx)^2)&(abs(Z-posZ)<radius(objectIdx)))=phantom((((Y-posY).^2+(X-posX).^2)<=radius(objectIdx)^2)&(abs(Z-posZ)<radius(objectIdx)))+delta(objectIdx);
- end
- otherwise
- error('Geometry not known: possible values are sphere, cube or cylinder')
- end
-end
-
-end
+function phantom = create3Dphantom(outputSize,settings)
+
+% this function creates a 3D phantom in which a specified number of objects
+% with a specified size and geometry are positioned; the positions of these
+% objects can be either predefined or random and also the additive strength
+% of each object can be chosen as a fixed or random value
+
+%{
+% Example: create a phantom with spheres of random size and strength
+settings=create3Dphantom;
+settings.geometry='sphere';
+settings.numObjects=30;
+settings.objectStrength=1;
+settings.distrStrength='random';
+settings.sizeObjects=100;
+settings.distrSize='random';
+phantom=create3Dphantom(512,settings);
+
+figure; imagesc(phantom(:,:,randi(size(phantom,3)))); colorbar;
+
+%}
+
+
+if (nargin < 2)
+ settings = struct;
+end
+
+% defaults
+defaults.sizeObjects = 0; % size of the objects; 0 leads to a size of 1/10 of the outputSize
+defaults.distrSize = 'fixed'; % switch between 'fixed' for fixed size or 'random' with a random size distribution with maximum settings.sizeObjects and minimum (settings.sizeObjects)/2
+defaults.numObjects = 10; % number of objects within the 3D volume
+defaults.geometry = 'sphere'; % geometry of the object; choices: sphere, cube and cylinder (with a random orientation along the three main axes)
+defaults.positions = 0; % positions of the object in space; 0 leads to a random distribution
+defaults.objectStrength = 1; % additive strength of the objects;
+defaults.distrStrength = 'fixed'; % switch between 'fixed' for fixed object strength or 'random' with a random strength between settings.objectStrength and (settings.objectStrength)/2 for each object
+
+% return default settings
+if (nargin == 0)
+ phantom=defaults;
+ return
+end
+
+settings = completeStruct(settings, defaults);
+
+% set size of the objects to 1/10 of the output size if no other settings
+% were specified
+if settings.sizeObjects == 0
+ settings.sizeObjects = floor(outputSize/10);
+end
+
+% determine the radius for each of the objects (either random or fixed value)
+if strcmp(settings.distrSize,'random')
+ radius=(settings.sizeObjects-rand(settings.numObjects,1)*0.5*settings.sizeObjects)/2;
+else
+ radius=(repmat(settings.sizeObjects,[settings.numObjects 1]))/2;
+end
+
+
+% determine the additive strength for each of the objects (either random or fixed value)
+if strcmp(settings.distrStrength,'random')
+ delta=settings.objectStrength-rand(settings.numObjects,1)*0.5*settings.objectStrength;
+else
+ delta=repmat(settings.objectStrength,[settings.numObjects 1]);
+end
+
+% determine random positions of objects within the 3D array if no positions
+% were given
+if settings.positions == 0
+ settings.positions = floor((rand(settings.numObjects,3))*(outputSize-2*max(radius)))+max(radius);
+end
+
+% create 3D object
+phantom = zeros(outputSize,outputSize,outputSize,'single');
+[X,Y,Z]=meshgrid(1:outputSize);
+
+% place the individual objects on their positions in the 3D volume
+for objectIdx=1:settings.numObjects
+ % position of the center of the object
+ posX = settings.positions(objectIdx,1);
+ posY = settings.positions(objectIdx,2);
+ posZ = settings.positions(objectIdx,3);
+
+ switch settings.geometry
+ case 'sphere'
+ phantom(((X-posX).^2+(Y-posY).^2+(Z-posZ).^2)<=radius(objectIdx)^2)=phantom(((X-posX).^2+(Y-posY).^2+(Z-posZ).^2)<=radius(objectIdx)^2)+delta(objectIdx);
+ case 'cube'
+ phantom((abs(X-posX)<radius(objectIdx))&(abs(Y-posY)<radius(objectIdx))&(abs(Z-posZ)<radius(objectIdx)))=phantom((abs(X-posX)<radius(objectIdx))&(abs(Y-posY)<radius(objectIdx))&(abs(Z-posZ)<radius(objectIdx)))+delta(objectIdx);
+ case 'cylinder'
+ orientation = rand*3; % determine orientation of the cylinder
+ if orientation < 1
+ phantom((((Y-posY).^2+(Z-posZ).^2)<=radius(objectIdx)^2)&(abs(X-posX)<radius(objectIdx)))=phantom((((Y-posY).^2+(Z-posZ).^2)<=radius(objectIdx)^2)&(abs(X-posX)<radius(objectIdx)))+delta(objectIdx);
+ elseif orientation < 2
+ phantom((((X-posX).^2+(Z-posZ).^2)<=radius(objectIdx)^2)&(abs(Y-posY)<radius(objectIdx)))=phantom((((X-posX).^2+(Z-posZ).^2)<=radius(objectIdx)^2)&(abs(Y-posY)<radius(objectIdx)))+delta(objectIdx);
+ else
+ phantom((((Y-posY).^2+(X-posX).^2)<=radius(objectIdx)^2)&(abs(Z-posZ)<radius(objectIdx)))=phantom((((Y-posY).^2+(X-posX).^2)<=radius(objectIdx)^2)&(abs(Z-posZ)<radius(objectIdx)))+delta(objectIdx);
+ end
+ otherwise
+ error('Geometry not known: possible values are sphere, cube or cylinder')
+ end
+end
+
+end
diff --git a/functions/tomography/radon3Darray.m b/functions/tomography/radon3Darray.m
index 87bf0ad522f5d72b1edd6cf24c2df8187b15928a..e93eccc29f221bda56349640dca4737fe85111f6 100755
--- a/functions/tomography/radon3Darray.m
+++ b/functions/tomography/radon3Darray.m
@@ -1,49 +1,49 @@
-%{
-% Example: Radon-transform of a cuboid-phatom
-vol = ones([132,112,84]);
-vol = padarray(vol, ([200,200,200]-size(vol))/2);
-angles = (0:179);
-settings.rotationAxis = 1;
-projs = radon3Darray(vol, angles, settings);
-
-figure; imagesc(projs(:,:,randi(numel(angles)))); colorbar;
-
-%}
-
-function projs = radon3Darray(vol, tomoAngles, settings)
-% performs the radon transform for an entire 3d object along a specified axis via a for-loop
-% over the respective dimension of the object; this results in 2d projections as
-%, e.g., acquired by area detectors
-
-if (nargin < 3)
- settings = struct;
-end
-
-% defaults
-defaults.outputSize = 0; % size of the projections, the default of 0 results in the default size of Matlab's radon transform
-defaults.rotationAxis = 3; % along which axis should the object be rotated; 1: along y-axis, 2: along x-axis, 3: along z-axis
-
-% return default settings
-if (nargin == 0)
- projs=defaults;
- return
-end
-
-settings = completeStruct(settings, defaults);
-
-% Input and output dimensions
-sizeVol = size(vol);
-numSlices = sizeVol(settings.rotationAxis);
-if settings.outputSize == 0
- sizeSlice = sizeVol(setdiff([1,2,3], settings.rotationAxis));
- settings.outputSize = 2*ceil(norm(sizeSlice-floor((sizeSlice-1)/2)-1))+3;
-end
-
-% Compute Radon transform for each slice
-projs = zeros([numSlices, settings.outputSize, numel(tomoAngles)], class(vol));
-for sliceIdx = 1:numSlices
- projs(sliceIdx,:,:) = radon(squeeze(getSlice(vol, sliceIdx, settings.rotationAxis)), tomoAngles);%, settings.outputSize);
-end
-
-end
-
+%{
+% Example: Radon-transform of a cuboid-phatom
+vol = ones([132,112,84]);
+vol = padarray(vol, ([200,200,200]-size(vol))/2);
+angles = (0:179);
+settings.rotationAxis = 1;
+projs = radon3Darray(vol, angles, settings);
+
+figure; imagesc(projs(:,:,randi(numel(angles)))); colorbar;
+
+%}
+
+function projs = radon3Darray(vol, tomoAngles, settings)
+% performs the radon transform for an entire 3d object along a specified axis via a for-loop
+% over the respective dimension of the object; this results in 2d projections as
+%, e.g., acquired by area detectors
+
+if (nargin < 3)
+ settings = struct;
+end
+
+% defaults
+defaults.outputSize = 0; % size of the projections, the default of 0 results in the default size of Matlab's radon transform
+defaults.rotationAxis = 3; % along which axis should the object be rotated; 1: along y-axis, 2: along x-axis, 3: along z-axis
+
+% return default settings
+if (nargin == 0)
+ projs=defaults;
+ return
+end
+
+settings = completeStruct(settings, defaults);
+
+% Input and output dimensions
+sizeVol = size(vol);
+numSlices = sizeVol(settings.rotationAxis);
+if settings.outputSize == 0
+ sizeSlice = sizeVol(setdiff([1,2,3], settings.rotationAxis));
+ settings.outputSize = 2*ceil(norm(sizeSlice-floor((sizeSlice-1)/2)-1))+3;
+end
+
+% Compute Radon transform for each slice
+projs = zeros([numSlices, settings.outputSize, numel(tomoAngles)], class(vol));
+for sliceIdx = 1:numSlices
+ projs(sliceIdx,:,:) = radon(squeeze(getSlice(vol, sliceIdx, settings.rotationAxis)), tomoAngles);%, settings.outputSize);
+end
+
+end
+
diff --git a/functions/tomography/ringremove.m b/functions/tomography/ringremove.m
index 0e07ea6a176624f051e859c01633808107d5eca3..ccf6758bb6b771ec5d28abc9407ec91d9c957012 100755
--- a/functions/tomography/ringremove.m
+++ b/functions/tomography/ringremove.m
@@ -1,106 +1,106 @@
-function sinoCorrected = ringremove(sino,settings)
-% this function can be used to remove ring artifacts in tomographic
-% reconstructions by reducing the amount of horizontal lines in the given
-% sinogram
-
-% ring removal is either based on the method proposed by Ketcham et al.
-% (Ketcham et al., Proc. SPIE 6318:63180O (2006)), denoted as 'additive' in
-% the script or on the method proposed by Muench et al. (Muench et al.,
-% Opt. Expr. 17(10):8567-8591 (2009)), denoted as 'wavelet'
-
-% large parts of the code for the wavelet-based approach are taken from the
-% original paper
-
-if nargin<2
- settings=struct;
-end
-
-defaults.method = 'additive'; % which ringremove method should be used? 'additive' or 'wavelet'
-
-defaults.additive_SmoothMethod='mean'; % smoothing method for averaged sinogram: moving averaging window ('mean'), moving median window ('median') or Gaussian filtering ('gauss')
-defaults.additive_WindowSize=10; % Window size for moving averaging or median window
-defaults.additive_GaussSigma=3; % sigma for the Gaussian filtering
-defaults.additive_Strength=1; % strength of the subtraction of the identified stripe artifacts from the original sinogram
-
-defaults.wavelet_Wname='sym8'; % wavelet base to be used for the decomposition of the sinogram
-defaults.wavelet_DecNum=3; % highest decomposition level
-defaults.wavelet_Sigma=1; % sigma for the Gaussian filtering of the horizontal part of the sinogram in Fourier space
-
-if (nargin == 0)
- sinoCorrected=defaults;
- return
-end
-
-settings=completeStruct(settings, defaults);
-
-
-switch settings.method
- case 'additive'
- numAngles = size(sino,2);
- % averaging of sinogram along the angle
- sinoMean = mean(sino,2);
-
- % smoothing of the sinogram
- switch settings.additive_SmoothMethod
- case 'mean'
- sinoMeanSmoothed = smoothdata(sinoMean,'movmean',settings.additive_WindowSize);
- case 'median'
- sinoMeanSmoothed = smoothdata(sinoMean,'movmedian',settings.additive_WindowSize);
- case 'gauss'
- sinoMeanSmoothed = imgaussfilt(sinoMean,settings.additive_GaussSigma);
- otherwise
- error('Choose between "mean", "median" and "gauss" for setting "additive_SmoothMethod"');
- end
- % identification of horizontal stripes
- correctionProfile = sinoMean-sinoMeanSmoothed;
- % creation of 2d correction mask
- correction2D = repmat(correctionProfile,[1 numAngles]);
- % sinogram correction
- sinoCorrected = sino-settings.additive_Strength*correction2D;
-
- case 'wavelet'
- numAngles = size(sino,2);
- sizeProjs = size(sino,1);
-
- % wavelet decomposition
- for decompIdx=1:settings.wavelet_DecNum
- [sino,Ch{decompIdx},Cv{decompIdx},Cd{decompIdx}]=dwt2(sino,settings.wavelet_Wname);
- end
-
- % FFT transform of horizontal frequency bands
- for decompIdx=1:settings.wavelet_DecNum
- % FFT
- fCh=fftshift(fft(Ch{decompIdx}'));
- [my,mx]=size(fCh);
-
- % damping of horizontal stripe information
- damp=1-exp(-[-floor(my/2):-floor(my/2)+my-1].^2/(2*settings.wavelet_Sigma^2));
- fCh=fCh.*repmat(damp',1,mx);
-
- % inverse FFT
- Ch{decompIdx}=ifft(ifftshift(fCh))';
- end
-
- % wavelet reconstruction
- sinoCorrected=sino;
- for decompIdx=settings.wavelet_DecNum:-1:1
- sinoCorrected=sinoCorrected(1:size(Cv{decompIdx},1),1:size(Cv{decompIdx},2));
- sinoCorrected=idwt2(sinoCorrected,Ch{decompIdx},Cv{decompIdx},Cd{decompIdx},settings.wavelet_Wname);
- end
-
- % make sure that corrected sinogram has the same size as the original sinogram
- if size(sinoCorrected,2) == numAngles+1
- sinoCorrected = sinoCorrected(:,1:end-1);
- end
- if size(sinoCorrected,1) == sizeProjs+1
- sinoCorrected = sinoCorrected(1:end-1,:);
- end
-
- otherwise
- error('Choose between ringremove based on "additive" or "wavelet" correction for setting "method"');
-end
-
-end
-
-
-
+function sinoCorrected = ringremove(sino,settings)
+% this function can be used to remove ring artifacts in tomographic
+% reconstructions by reducing the amount of horizontal lines in the given
+% sinogram
+
+% ring removal is either based on the method proposed by Ketcham et al.
+% (Ketcham et al., Proc. SPIE 6318:63180O (2006)), denoted as 'additive' in
+% the script or on the method proposed by Muench et al. (Muench et al.,
+% Opt. Expr. 17(10):8567-8591 (2009)), denoted as 'wavelet'
+
+% large parts of the code for the wavelet-based approach are taken from the
+% original paper
+
+if nargin<2
+ settings=struct;
+end
+
+defaults.method = 'additive'; % which ringremove method should be used? 'additive' or 'wavelet'
+
+defaults.additive_SmoothMethod='mean'; % smoothing method for averaged sinogram: moving averaging window ('mean'), moving median window ('median') or Gaussian filtering ('gauss')
+defaults.additive_WindowSize=10; % Window size for moving averaging or median window
+defaults.additive_GaussSigma=3; % sigma for the Gaussian filtering
+defaults.additive_Strength=1; % strength of the subtraction of the identified stripe artifacts from the original sinogram
+
+defaults.wavelet_Wname='sym8'; % wavelet base to be used for the decomposition of the sinogram
+defaults.wavelet_DecNum=3; % highest decomposition level
+defaults.wavelet_Sigma=1; % sigma for the Gaussian filtering of the horizontal part of the sinogram in Fourier space
+
+if (nargin == 0)
+ sinoCorrected=defaults;
+ return
+end
+
+settings=completeStruct(settings, defaults);
+
+
+switch settings.method
+ case 'additive'
+ numAngles = size(sino,2);
+ % averaging of sinogram along the angle
+ sinoMean = mean(sino,2);
+
+ % smoothing of the sinogram
+ switch settings.additive_SmoothMethod
+ case 'mean'
+ sinoMeanSmoothed = smoothdata(sinoMean,'movmean',settings.additive_WindowSize);
+ case 'median'
+ sinoMeanSmoothed = smoothdata(sinoMean,'movmedian',settings.additive_WindowSize);
+ case 'gauss'
+ sinoMeanSmoothed = imgaussfilt(sinoMean,settings.additive_GaussSigma);
+ otherwise
+ error('Choose between "mean", "median" and "gauss" for setting "additive_SmoothMethod"');
+ end
+ % identification of horizontal stripes
+ correctionProfile = sinoMean-sinoMeanSmoothed;
+ % creation of 2d correction mask
+ correction2D = repmat(correctionProfile,[1 numAngles]);
+ % sinogram correction
+ sinoCorrected = sino-settings.additive_Strength*correction2D;
+
+ case 'wavelet'
+ numAngles = size(sino,2);
+ sizeProjs = size(sino,1);
+
+ % wavelet decomposition
+ for decompIdx=1:settings.wavelet_DecNum
+ [sino,Ch{decompIdx},Cv{decompIdx},Cd{decompIdx}]=dwt2(sino,settings.wavelet_Wname);
+ end
+
+ % FFT transform of horizontal frequency bands
+ for decompIdx=1:settings.wavelet_DecNum
+ % FFT
+ fCh=fftshift(fft(Ch{decompIdx}'));
+ [my,mx]=size(fCh);
+
+ % damping of horizontal stripe information
+ damp=1-exp(-[-floor(my/2):-floor(my/2)+my-1].^2/(2*settings.wavelet_Sigma^2));
+ fCh=fCh.*repmat(damp',1,mx);
+
+ % inverse FFT
+ Ch{decompIdx}=ifft(ifftshift(fCh))';
+ end
+
+ % wavelet reconstruction
+ sinoCorrected=sino;
+ for decompIdx=settings.wavelet_DecNum:-1:1
+ sinoCorrected=sinoCorrected(1:size(Cv{decompIdx},1),1:size(Cv{decompIdx},2));
+ sinoCorrected=idwt2(sinoCorrected,Ch{decompIdx},Cv{decompIdx},Cd{decompIdx},settings.wavelet_Wname);
+ end
+
+ % make sure that corrected sinogram has the same size as the original sinogram
+ if size(sinoCorrected,2) == numAngles+1
+ sinoCorrected = sinoCorrected(:,1:end-1);
+ end
+ if size(sinoCorrected,1) == sizeProjs+1
+ sinoCorrected = sinoCorrected(1:end-1,:);
+ end
+
+ otherwise
+ error('Choose between ringremove based on "additive" or "wavelet" correction for setting "method"');
+end
+
+end
+
+
+