Hi Luis,<br>      thanks for the reply again! Unfortunately I cannot post the data or the results, it has got commercial value (**proprietary**)! <br><br>You are absolutely correct with regards to messing around with the the time step though! It shouldn&#39;t be much of a problem running with a smaller time step, the only caveat being a longer time for computations! That being said that is exactly what I did for testing too and I did not find anything suspect qualitatively (like you said) or quantitatively (in a rms perspective between the voxels). Maybe I&#39;ll just code up a implicit method for double checking. In terms of development though I am certainly going to stick with the smaller time step restriction!<br>
<br>Please forgive my thread etiquette about this next question but I noticed a strange memory issue when I run anisotropic diffusion on a dataset that is 1024^3. The dataset is 2GB but I noticed that 14GB was being used! Isn&#39;t this odd? I have not probed the code to find out where and why, so if there is an obvious answer my apologies, just seeing if there is a quick and dirty answer!<br>
<br>Cheers,<br><br>C.S.N<br>
<br><div class="gmail_quote">On Wed, Feb 17, 2010 at 5:24 PM, Luis Ibanez <span dir="ltr">&lt;<a href="mailto:luis.ibanez@kitware.com" target="_blank">luis.ibanez@kitware.com</a>&gt;</span> wrote:<br><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">

Hi Natarajan,<br>
<br>
Thanks for letting us know of your experiments.<br>
<br>
I would be curious to see if the images that you used<br>
as input had enough high contrast edges and thin spatial<br>
features that would make visible the artifact that we are<br>
concerned about.<br>
<br>
Being pragmatic, you should be able to obtain a similar<br>
degree of smoothing by running the filter with the smaller<br>
time-step that is currently defined in ITK, and double the<br>
number of iterations.<br>
<br>
The limit in the time-step is the one at which the conditions<br>
of stability break down. Playing in that borderline may not<br>
be a good way of using the code.<br>
<br>
<br>
 BTW:Thanks for posting the link to the report.<br>
<br>
<br>
            Regards,<br>
<br>
<br>
                   Luis<br>
<br>
<br>
----------------------------------------------------------------------------------<br>
<div><div></div><div>On Tue, Feb 16, 2010 at 11:14 PM, Natarajan CS &lt;<a href="mailto:csnataraj@gmail.com" target="_blank">csnataraj@gmail.com</a>&gt; wrote:<br>
&gt; Hello Luis,<br>
&gt;      just an update on the time step question. For my data I don&#39;t see any<br>
&gt; undesirable side effects like what you described. I also tested<br>
&gt; corresponding values for a 2D image, same thing! It seems that &lt;=<br>
&gt; 1/2^dimension should work! However, this comes with the caveat that image<br>
&gt; analysis is not my expertise!! It would be nice to get a affirmation from a<br>
&gt; more reliable source!<br>
&gt;<br>
&gt; Nonetheless, I really appreciate you guys making ITK open source and being<br>
&gt; great with answering questions!<br>
&gt;<br>
&gt; Cheers,<br>
&gt;<br>
&gt; C.S.N<br>
&gt;<br>
&gt;<br>
&gt; On Mon, Feb 15, 2010 at 5:01 PM, Natarajan CS &lt;<a href="mailto:csnataraj@gmail.com" target="_blank">csnataraj@gmail.com</a>&gt; wrote:<br>
&gt;&gt;<br>
&gt;&gt; Hi Luis,<br>
&gt;&gt;       sure, sounds like a good plan. I am heading out right now, but<br>
&gt;&gt; hopefully I should be able to get it done by tomorrow! I will try with time<br>
&gt;&gt; steps at 1/2^N, 1/2^(N+1), and 1.5/2^N and send out another e-mail.<br>
&gt;&gt;<br>
&gt;&gt; Cheers,<br>
&gt;&gt;<br>
&gt;&gt; C.S.N<br>
&gt;&gt;<br>
&gt;&gt; P.S :- Just in case anyone is interested, the paper is available as a tech<br>
&gt;&gt; report from Berkeley<br>
&gt;&gt;<br>
&gt;&gt; <a href="http://digitalassets.lib.berkeley.edu/techreports/ucb/text/CSD-88-483.pdf" target="_blank">http://digitalassets.lib.berkeley.edu/techreports/ucb/text/CSD-88-483.pdf</a><br>
&gt;&gt;<br>
&gt;&gt;<br>
&gt;&gt; On Mon, Feb 15, 2010 at 3:58 PM, Luis Ibanez &lt;<a href="mailto:luis.ibanez@kitware.com" target="_blank">luis.ibanez@kitware.com</a>&gt;<br>
&gt;&gt; wrote:<br>
&gt;&gt;&gt;<br>
&gt;&gt;&gt; Hi Natarajan,<br>
&gt;&gt;&gt;<br>
&gt;&gt;&gt; You are right, the code is actually testing for<br>
&gt;&gt;&gt;<br>
&gt;&gt;&gt;                        1/2^(N+1)<br>
&gt;&gt;&gt;<br>
&gt;&gt;&gt; not for<br>
&gt;&gt;&gt;<br>
&gt;&gt;&gt;                           1/2^(N)<br>
&gt;&gt;&gt;<br>
&gt;&gt;&gt; as I incorrectly wrote in my email.<br>
&gt;&gt;&gt;<br>
&gt;&gt;&gt;<br>
&gt;&gt;&gt; This filter is intended to be the implementation of<br>
&gt;&gt;&gt; the Perona-Malik method, as described in:<br>
&gt;&gt;&gt;<br>
&gt;&gt;&gt;  * Pietro Perona and Jalhandra Malik, ``Scale-space and edge detection<br>
&gt;&gt;&gt; using<br>
&gt;&gt;&gt;  * anisotropic diffusion,&#39;&#39; IEEE Transactions on Pattern Analysis Machine<br>
&gt;&gt;&gt;  * Intelligence, vol. 12, pp. 629-639, 1990.<br>
&gt;&gt;&gt;<br>
&gt;&gt;&gt;<br>
&gt;&gt;&gt; I don&#39;t have at hand the Perona Malik paper,<br>
&gt;&gt;&gt; so let&#39;s do this:<br>
&gt;&gt;&gt;<br>
&gt;&gt;&gt;<br>
&gt;&gt;&gt; Could you please modify the code locally in your ITK code,<br>
&gt;&gt;&gt; and run it with a couple of cases, using a timestep that is<br>
&gt;&gt;&gt; in the range:<br>
&gt;&gt;&gt;<br>
&gt;&gt;&gt;                    1/2^(N+1)  &lt;  t  &lt; 1/2^(N)<br>
&gt;&gt;&gt;<br>
&gt;&gt;&gt; and let us know if the output looks correct ?<br>
&gt;&gt;&gt;<br>
&gt;&gt;&gt; ---<br>
&gt;&gt;&gt;<br>
&gt;&gt;&gt; Typically when you use a timestep that is too large, you<br>
&gt;&gt;&gt; will see intensity artifacts that look like rhomboids in the<br>
&gt;&gt;&gt; borders of the objects in the image.<br>
&gt;&gt;&gt;<br>
&gt;&gt;&gt;<br>
&gt;&gt;&gt;     Please let us know what you find,<br>
&gt;&gt;&gt;<br>
&gt;&gt;&gt;<br>
&gt;&gt;&gt;            Thanks<br>
&gt;&gt;&gt;<br>
&gt;&gt;&gt;<br>
&gt;&gt;&gt;                 Luis<br>
&gt;&gt;&gt;<br>
&gt;&gt;&gt;<br>
&gt;&gt;&gt;<br>
&gt;&gt;&gt; ---------------------------------------------------------------------------<br>
&gt;&gt;&gt; On Mon, Feb 15, 2010 at 3:34 PM, Natarajan CS &lt;<a href="mailto:csnataraj@gmail.com" target="_blank">csnataraj@gmail.com</a>&gt;<br>
&gt;&gt;&gt; wrote:<br>
&gt;&gt;&gt; &gt; Hi Luis,<br>
&gt;&gt;&gt; &gt;      thanks for the quick reply. Nope, I do not believe I have enabled<br>
&gt;&gt;&gt; &gt; this<br>
&gt;&gt;&gt; &gt; flag. I am confused, I did take a look at the code snippet you<br>
&gt;&gt;&gt; &gt; suggested and<br>
&gt;&gt;&gt; &gt; figured that time step is &lt;= 1/2^(N+1) not 1/2^N (Judging from line 75<br>
&gt;&gt;&gt; &gt; and/through 81). In reality I believe this has to be 1/2^N as you<br>
&gt;&gt;&gt; &gt; suggested,<br>
&gt;&gt;&gt; &gt; at least if the implementation is same as Perona and Malik..<br>
&gt;&gt;&gt; &gt;<br>
&gt;&gt;&gt; &gt; just in case my metadata looks as follows :<br>
&gt;&gt;&gt; &gt;<br>
&gt;&gt;&gt; &gt; NDims = 3<br>
&gt;&gt;&gt; &gt; DimSize = 256 256 256<br>
&gt;&gt;&gt; &gt; ElementType = MET_SHORT<br>
&gt;&gt;&gt; &gt; ElementDataFile = raw256.raw<br>
&gt;&gt;&gt; &gt;<br>
&gt;&gt;&gt; &gt;<br>
&gt;&gt;&gt; &gt; Again, thanks for the help!<br>
&gt;&gt;&gt; &gt;<br>
&gt;&gt;&gt; &gt; Cheers,<br>
&gt;&gt;&gt; &gt;<br>
&gt;&gt;&gt; &gt; C.S.N<br>
&gt;&gt;&gt; &gt;<br>
&gt;&gt;&gt; &gt;<br>
&gt;&gt;&gt; &gt;<br>
&gt;&gt;&gt; &gt; On Mon, Feb 15, 2010 at 2:00 PM, Luis Ibanez &lt;<a href="mailto:luis.ibanez@kitware.com" target="_blank">luis.ibanez@kitware.com</a>&gt;<br>
&gt;&gt;&gt; &gt; wrote:<br>
&gt;&gt;&gt; &gt;&gt;<br>
&gt;&gt;&gt; &gt;&gt; Hi Natarajan,<br>
&gt;&gt;&gt; &gt;&gt;<br>
&gt;&gt;&gt; &gt;&gt;<br>
&gt;&gt;&gt; &gt;&gt;<br>
&gt;&gt;&gt; &gt;&gt;                      &quot;Glimpsing at the Source, leaves no doubt&quot;<br>
&gt;&gt;&gt; &gt;&gt;<br>
&gt;&gt;&gt; &gt;&gt;<br>
&gt;&gt;&gt; &gt;&gt; This is what the file:<br>
&gt;&gt;&gt; &gt;&gt;<br>
&gt;&gt;&gt; &gt;&gt;<br>
&gt;&gt;&gt; &gt;&gt; Insight/Code/BasicFilters/itkAnisotropicDiffusionImageFilter.txx<br>
&gt;&gt;&gt; &gt;&gt;<br>
&gt;&gt;&gt; &gt;&gt; has in lines 75-82:<br>
&gt;&gt;&gt; &gt;&gt;<br>
&gt;&gt;&gt; &gt;&gt;<br>
&gt;&gt;&gt; &gt;&gt;  if ( m_TimeStep &gt;  (minSpacing / vcl_pow(2.0,<br>
&gt;&gt;&gt; &gt;&gt; static_cast&lt;double&gt;(ImageDimension) + 1))  )<br>
&gt;&gt;&gt; &gt;&gt;    {<br>
&gt;&gt;&gt; &gt;&gt;    //    f-&gt;SetTimeStep(1.0 / vcl_pow(2.0,<br>
&gt;&gt;&gt; &gt;&gt; static_cast&lt;double&gt;(ImageDimension)));<br>
&gt;&gt;&gt; &gt;&gt;    itkWarningMacro( &lt;&lt; &quot;Anisotropic diffusion unstable time step: &quot;<br>
&gt;&gt;&gt; &gt;&gt;                     &lt;&lt; m_TimeStep &lt;&lt; std::endl<br>
&gt;&gt;&gt; &gt;&gt;                     &lt;&lt; &quot;Stable time step for this image must be<br>
&gt;&gt;&gt; &gt;&gt; smaller<br>
&gt;&gt;&gt; &gt;&gt; than &quot;<br>
&gt;&gt;&gt; &gt;&gt;                     &lt;&lt; minSpacing / vcl_pow(2.0,<br>
&gt;&gt;&gt; &gt;&gt; static_cast&lt;double&gt;(ImageDimension+1)));<br>
&gt;&gt;&gt; &gt;&gt;    }<br>
&gt;&gt;&gt; &gt;&gt;<br>
&gt;&gt;&gt; &gt;&gt;<br>
&gt;&gt;&gt; &gt;&gt; In summary:<br>
&gt;&gt;&gt; &gt;&gt;<br>
&gt;&gt;&gt; &gt;&gt;                    Timestep    &lt;     pixelSpacing  /  ( 2^N )<br>
&gt;&gt;&gt; &gt;&gt;<br>
&gt;&gt;&gt; &gt;&gt; where N is the image dimension (e.g. 3 for 3D images, and 2 for 2D<br>
&gt;&gt;&gt; &gt;&gt; images).<br>
&gt;&gt;&gt; &gt;&gt;<br>
&gt;&gt;&gt; &gt;&gt;<br>
&gt;&gt;&gt; &gt;&gt;<br>
&gt;&gt;&gt; &gt;&gt;<br>
&gt;&gt;&gt; &gt;&gt; Did you enable the flag:<br>
&gt;&gt;&gt; &gt;&gt;<br>
&gt;&gt;&gt; &gt;&gt;               UseImageSpacing  ?<br>
&gt;&gt;&gt; &gt;&gt;<br>
&gt;&gt;&gt; &gt;&gt; If so, what is the pixel spacing of your image ?<br>
&gt;&gt;&gt; &gt;&gt;<br>
&gt;&gt;&gt; &gt;&gt;<br>
&gt;&gt;&gt; &gt;&gt;      Regards,<br>
&gt;&gt;&gt; &gt;&gt;<br>
&gt;&gt;&gt; &gt;&gt;<br>
&gt;&gt;&gt; &gt;&gt;<br>
&gt;&gt;&gt; &gt;&gt;            Luis<br>
&gt;&gt;&gt; &gt;&gt;<br>
&gt;&gt;&gt; &gt;&gt;<br>
&gt;&gt;&gt; &gt;&gt; -------------------------------------------<br>
&gt;&gt;&gt; &gt;&gt; On Mon, Feb 15, 2010 at 12:27 PM, Natarajan CS &lt;<a href="mailto:csnataraj@gmail.com" target="_blank">csnataraj@gmail.com</a>&gt;<br>
&gt;&gt;&gt; &gt;&gt; wrote:<br>
&gt;&gt;&gt; &gt;&gt; &gt; Hello all,<br>
&gt;&gt;&gt; &gt;&gt; &gt;       Not sure if this is a bug or a implementation detail I am<br>
&gt;&gt;&gt; &gt;&gt; &gt; overlooking,<br>
&gt;&gt;&gt; &gt;&gt; &gt; I am tending towards the latter, so any info would be helpful !<br>
&gt;&gt;&gt; &gt;&gt; &gt;  I get  the following warning from the<br>
&gt;&gt;&gt; &gt;&gt; &gt; GradientAnisotropicDiffusionImageFilter (The filter is being applied<br>
&gt;&gt;&gt; &gt;&gt; &gt; to<br>
&gt;&gt;&gt; &gt;&gt; &gt; a<br>
&gt;&gt;&gt; &gt;&gt; &gt; 256^3 dataset of type short)<br>
&gt;&gt;&gt; &gt;&gt; &gt;<br>
&gt;&gt;&gt; &gt;&gt; &gt; WARNING: In<br>
&gt;&gt;&gt; &gt;&gt; &gt;<br>
&gt;&gt;&gt; &gt;&gt; &gt;<br>
&gt;&gt;&gt; &gt;&gt; &gt; /hpc/soft/natac0/itk/include/InsightToolkit/BasicFilters/itkAnisotropicDiffusionImageFilter.txx,<br>
&gt;&gt;&gt; &gt;&gt; &gt; line 78<br>
&gt;&gt;&gt; &gt;&gt; &gt; GradientAnisotropicDiffusionImageFilter (0xbc4120): Anisotropic<br>
&gt;&gt;&gt; &gt;&gt; &gt; diffusion<br>
&gt;&gt;&gt; &gt;&gt; &gt; unstable time step: 0.07<br>
&gt;&gt;&gt; &gt;&gt; &gt; Stable time step for this image must be smaller than 0.0625<br>
&gt;&gt;&gt; &gt;&gt; &gt;<br>
&gt;&gt;&gt; &gt;&gt; &gt; i was under the impression that the ITK implimentation was a<br>
&gt;&gt;&gt; &gt;&gt; &gt; explicit<br>
&gt;&gt;&gt; &gt;&gt; &gt; method<br>
&gt;&gt;&gt; &gt;&gt; &gt; with the time step requirement for 3D images is &lt;= 1/8 and &lt;= 1/4<br>
&gt;&gt;&gt; &gt;&gt; &gt; for<br>
&gt;&gt;&gt; &gt;&gt; &gt; 2D. Am<br>
&gt;&gt;&gt; &gt;&gt; &gt; i completely off-track here?<br>
&gt;&gt;&gt; &gt;&gt; &gt;<br>
&gt;&gt;&gt; &gt;&gt; &gt; Appreciate any help!<br>
&gt;&gt;&gt; &gt;&gt; &gt;<br>
&gt;&gt;&gt; &gt;&gt; &gt; Thanks,<br>
&gt;&gt;&gt; &gt;&gt; &gt;<br>
&gt;&gt;&gt; &gt;&gt; &gt; C.S.N<br>
&gt;&gt;&gt; &gt;&gt; &gt;<br>
&gt;&gt;&gt; &gt;&gt; &gt; _____________________________________<br>
&gt;&gt;&gt; &gt;&gt; &gt; Powered by <a href="http://www.kitware.com" target="_blank">www.kitware.com</a><br>
&gt;&gt;&gt; &gt;&gt; &gt;<br>
&gt;&gt;&gt; &gt;&gt; &gt; Visit other Kitware open-source projects at<br>
&gt;&gt;&gt; &gt;&gt; &gt; <a href="http://www.kitware.com/opensource/opensource.html" target="_blank">http://www.kitware.com/opensource/opensource.html</a><br>
&gt;&gt;&gt; &gt;&gt; &gt;<br>
&gt;&gt;&gt; &gt;&gt; &gt; Kitware offers ITK Training Courses, for more information visit:<br>
&gt;&gt;&gt; &gt;&gt; &gt; <a href="http://www.kitware.com/products/protraining.html" target="_blank">http://www.kitware.com/products/protraining.html</a><br>
&gt;&gt;&gt; &gt;&gt; &gt;<br>
&gt;&gt;&gt; &gt;&gt; &gt; Please keep messages on-topic and check the ITK FAQ at:<br>
&gt;&gt;&gt; &gt;&gt; &gt; <a href="http://www.itk.org/Wiki/ITK_FAQ" target="_blank">http://www.itk.org/Wiki/ITK_FAQ</a><br>
&gt;&gt;&gt; &gt;&gt; &gt;<br>
&gt;&gt;&gt; &gt;&gt; &gt; Follow this link to subscribe/unsubscribe:<br>
&gt;&gt;&gt; &gt;&gt; &gt; <a href="http://www.itk.org/mailman/listinfo/insight-users" target="_blank">http://www.itk.org/mailman/listinfo/insight-users</a><br>
&gt;&gt;&gt; &gt;&gt; &gt;<br>
&gt;&gt;&gt; &gt;&gt; &gt;<br>
&gt;&gt;&gt; &gt;<br>
&gt;&gt;&gt; &gt;<br>
&gt;&gt;<br>
&gt;<br>
&gt;<br>
</div></div></blockquote></div><br>