<div dir="ltr">I guess the obvious common traps need to be checked - confirm you are in Release mode. Things like neighbourhood iterators can be very expensive in debug mode. Also, the neighbourhood iterator will be doing some sort of edge of image safety check, which your other code may not be doing (I haven't checked). The neighbourhood iterator isn't, I guess, really designed to be resized every time it is used. Normally you set it up once and apply it lots of times. I've no clue as to how expensive the initialisation is, but it may pay to keep a structure containing iterators that you've already allocated so that you can pick the right one and move it.<div><br></div><div>It may turn out that the design is really bad for your problem, with performance being destroyed by the creation phase, that isn't an issue in the normal use case. Continuing to use it will give you edge safety and dimension independence, but I guess you'll need to address the issue of initialization by being clever about the reuse of iterators.</div><div><br></div><div>One way to do that would be to use a pass to record the location and value of all of the skeleton points, sort them by brightness, then visit in order so that neighbourhoods can be reused without resizing. This assumes, of course that the skeletons contain much fewer voxels than the image.</div><div><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Mon, Mar 30, 2015 at 6:51 PM, Johannes Weber <span dir="ltr"><<a href="mailto:JohannesWeber@gmx.at" target="_blank">JohannesWeber@gmx.at</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><div style="font-family:Verdana;font-size:12.0px"><div>
<div>thank you for your answers!</div>
<div> </div>
<div>I tried out following idea of yours: "If you use neighbourhoods, which is the appropriate ITK structure, then avoid visiting all voxels with the neighbourhood iterator. Use a standard single voxel iterator to find the skeleton voxels, then move the neighbourhood iterator to that location."</div>
<div>this is fast until I use the SetRadius method of the neighborhood iterator... it goes up form about 2 seconds to 1 minute.</div>
<div>
<div>for(iterator.GoToBegin(); !iterator.IsAtEnd(); iterator++)<br>
{<br>
if(iterator.Get() > 0)<br>
{<br>
neighborhoodIterator.SetLocation(iterator.GetIndex());<br>
radius.Fill( (int) iterator.Get());<br>
neighborhoodIterator.SetRadius(radius);<br>
}</div>
<div> }</div>
</div>
<div> </div>
<div>Than I implemented what I want to do by direct accessing the image array:</div>
<div>float * imgArray = distanceRidge->GetBufferPointer();</div>
<div> </div>
<div>and for every ridge point (value > 0) in the image I go through the desired region with 3 for loops:</div>
<div> </div>
<div>
<div>for (int i1 = iStart; i1 <= iStop; i1++)<br>
{<br>
r1SquaredK = (i1 - i) * (i1 - i);</div>
<div> for (int j1 = jStart; j1 <= jStop; j1++)<br>
{<br>
r1SquaredJK = r1SquaredK + (j1 - j) * (j1 - j);<br>
if(r1SquaredJK <= rSquared)<br>
{<br>
for (int k1 = kStart; k1 <= kStop; k1++)<br>
{ <br>
r1Squared = r1SquaredJK + (k1 - k) * (k1 - k);<br>
if (r1Squared <= rSquared)<br>
{ <br>
s1 = imgArray[i1 + imageSize[0] * (j1 + imageSize[1] * k1)];<br>
if(rSquared > s1)<br>
{<br>
imgArray[i1 + imageSize[0] * (j1 + imageSize[1] * k1)] = rSquared; <br>
} <br>
}<br>
}<br>
}<br>
} <br>
}</div>
</div>
<div> </div>
<div>and this works fine! and it is really fast! with this approach I am getting there where I want to be.</div>
<div>I am doing the same with the neighborhood iterator, or I also tried out doing it with a image region iterator and instead of going through the desired region with 3 for loops I always set the region of the iterator (see my first email).</div>
<div>So what I am asking myself is now, if it is so much faster using directly the flat array where the image data is saved... is it also a proper way to do it like this? and is there a different way of implementing this, which is as fast as with flat array but using iterators from ITK or other functions/options from ITK?</div>
<div>
<div name="quote" style="margin:10px 5px 5px 10px;padding:10px 0 10px 10px;border-left:2px solid #c3d9e5;word-wrap:break-word">
<div style="margin:0 0 10px 0"><b>Gesendet:</b> Freitag, 27. März 2015 um 00:47 Uhr<br>
<b>Von:</b> "Richard Beare" <<a href="mailto:richard.beare@gmail.com" target="_blank">richard.beare@gmail.com</a>><br>
<b>An:</b> <a href="mailto:Richard.Beare@ieee.org" target="_blank">Richard.Beare@ieee.org</a><br>
<b>Cc:</b> "Johannes Weber" <<a href="mailto:JohannesWeber@gmx.at" target="_blank">JohannesWeber@gmx.at</a>>, "Insight Users" <<a href="mailto:insight-users@itk.org" target="_blank">insight-users@itk.org</a>><br>
<b>Betreff:</b> Re: [ITK-users] setting big spherical neighborhoods to center voxel value</div>
<div name="quoted-content">
<div>
<div>
<div>
<div>
<div>
<div>
<div>
<div>
<div>
<div>
<div>
<div>
<div>
<div>Here's another option - I haven't thought it through completely so could be far from the mark, but I think it is close.<br>
</div>
<div>consider the following:<br>
</div>
</div>
A - input mask</div>
B - distance transform of mask (positive values inside the mask)</div>
C - ridge lines in B (sort of skeleton)<br>
</div>
You are aiming to create an image of spheres, centres on voxels in C, with size set by the "brightness" at centre location in C.<br>
</div>
Suppose you can create the spatially variant dilation by a sphere using the parabolic tools I mentioned in the previous message - that gives you another mask:<br>
</div>
D - spatially variant dilation of C<br>
</div>
Now invert C and take the distance transform to give E.</div>
take a distance transform of D and add to E, then mask by D. The idea here is that if you add two distance transforms together you end up with a flat surface, at least in the case where the second DT is computed from peaks in the first. <br>
</div>
The result should be a sort of tube, of varying thickness, where the voxel value corresponds to the tube thickness (possibly with a uniform offset). Note that I don't think you will see sphere corresponding to the brightest skeleton points using this approach, but the result may be useful for your application.<br>
<br>
</div>
Also, a couple of notes regarding the direct approach.<br>
</div>
<div>If you use neighbourhoods, which is the appropriate ITK structure, then avoid visiting all voxels with the neighbourhood iterator. Use a standard single voxel iterator to find the skeleton voxels, then move the neighbourhood iterator to that location.</div>
<div> </div>
<div> </div>
</div>
There are also some optimizations you can use if you do some checks of the change in sphere size between neighbours. If the spheres were all identically sized then it is feasible to construct lists of voxels that are not in common between neighbouring spheres, and just visit those, which is a big saving. Much trickier when the sizes aren't the same, but I'd guess that for a given skeleton there'd be a relatively small number of combinations of neighbouring sizes.
<div>
<div>
<div>
<div>
<div>
<div>
<div>
<div>
<div> </div>
<div> </div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
<div class="gmail_extra">
<div class="gmail_quote">On Fri, Mar 27, 2015 at 8:02 AM, Richard Beare <span><<a href="http://richard.beare@gmail.com" target="_blank">richard.beare@gmail.com</a>></span> wrote:
<blockquote class="gmail_quote" style="margin:0 0 0 0.8ex;border-left:1.0px rgb(204,204,204) solid;padding-left:1.0ex">
<div>This isn't a complete solution, but might give you some ideas.
<div> </div>
<div>The part about using the skeleton points, which are distance transform voxels, as the source for the radius, can be implemented using spatially variant morphology. I have a short publication about doing this efficiently with parabolic functions. However this doesn't propagate the voxel value, only producing a mask. I have thought about combining with label dilation, to propagate a label value, but haven't done it yet.</div>
<div> </div>
<div>A queue based approach to something similar is discussed in:</div>
<div>"Labelled reconstruction of binary objects: a vector propagation algorithm", buckley and lagerstrom </div>
</div>
<div>
<div>
<div class="gmail_extra">
<div class="gmail_quote">On Thu, Mar 26, 2015 at 11:36 PM, Bradley Lowekamp <span><<a href="http://blowekamp@mail.nih.gov" target="_blank">blowekamp@mail.nih.gov</a>></span> wrote:
<blockquote class="gmail_quote" style="margin:0 0 0 0.8ex;border-left:1.0px rgb(204,204,204) solid;padding-left:1.0ex">
<div>Hello,
<div> </div>
<div>If you have a 3D image and you are visiting a neighborhood of size 20^3, you are doing 8000 visits per pixel there is no way to make this efficient. You have a big O algorithm problem.</div>
<div> </div>
<div>The Neighborhood iterator would be a the more ITK way of doing what you are trying to do [1] [2]. But that's is not the order of efficiency improvement you need.</div>
<div> </div>
<div>You need to revise your algorithm so you only visit each pixel once. Perhaps with region growing and queues, or auxiliary images to keep track of the current distance or other data.</div>
<div> </div>
<div>Hope that helps,</div>
<div>Brad</div>
<div> </div>
<div>[1] <a href="http://www.itk.org/Doxygen/html/classitk_1_1NeighborhoodIterator.html" target="_blank">http://www.itk.org/Doxygen/html/classitk_1_1NeighborhoodIterator.html</a></div>
<div>[2] <a href="http://itk.org/Wiki/ITK/Examples/Iterators/NeighborhoodIterator" target="_blank">http://itk.org/Wiki/ITK/Examples/Iterators/NeighborhoodIterator</a></div>
<div>
<div>
<div>On Mar 26, 2015, at 5:24 AM, <a href="http://JohannesWeber@gmx.at" target="_blank">JohannesWeber@gmx.at</a> wrote:</div>
<blockquote>
<div>
<div style="font-family:Verdana;font-size:12.0px">
<div>Hello everyone,</div>
<div> </div>
<div>I have the following problem: after calculating the distance map (e.g. with DanielssonDistanceMapImageFilter) I am getting rid of most of the voxels (= setting 0) after calculating a so called distance ridge (kind of a skeleton). Now I want for every voxel of this distance ridge that it is a center voxel for a spherical neighborhood with the radius equal to the distance value of the voxel, and all voxels included in this sphere are set to the distance value of the voxel of the distance ridge (the center voxel of the sphere). So that means the neighborhoods can become big, e.g. radius of 10, or 20 voxels. The problem is here the performance... I implemented it somehow, but the performance nowhere near it should be.</div>
<div>e.g. going through the image with a neighborhood iterator and vor every voxel bigger than 0 creating a neighborhood with the radius with the distance value of this voxel that seems to take very long alone to create and indexing the neighborhood.</div>
<div>another approach I tried is to extract all the voxels of the distance ridge, iterate through them and calculate for every ridge voxel the region and iterate through the region doing proper calculations:</div>
<div> </div>
<div>
<div>for (int rp = 0; rp < nRidgePoints; rp++)<br>
{<br>
ImageType::IndexType s1Index;<br>
const int i = ridgePointsIndex[0][rp];<br>
const int j = ridgePointsIndex[1][rp];<br>
const int k = ridgePointsIndex[2][rp];<br>
const float r = ridgePointsValues[rp];</div>
<div><br>
rSquared = (int) ((r * r) + 0.5f);<br>
rInt = (int) r;</div>
<div> if(rInt < r) rInt++;</div>
<div> iStart = i - rInt;<br>
if(iStart < 0) iStart = 0;<br>
iStop = i + rInt;<br>
if(iStop >= imageSize[0]) iStop = imageSize[0] - 1;</div>
<div> jStart = j - rInt;<br>
if(jStart < 0) jStart = 0;<br>
jStop = j + rInt;<br>
if(jStop >= imageSize[1]) jStop = imageSize[1] - 1;</div>
<div> kStart = k - rInt;<br>
if(kStart < 0) kStart = 0;<br>
kStop = k + rInt;<br>
if(kStop >= imageSize[2]) kStop = imageSize[2] - 1;</div>
<div> ImageType::IndexType index;<br>
ImageType::SizeType size;</div>
<div> index[0] = iStart;<br>
index[1] = jStart;<br>
index[2] = kStart;</div>
<div> size[0] = iStop - iStart + 1;<br>
size[1] = jStop - jStart + 1;<br>
size[2] = kStop - kStart + 1;</div>
<div> ImageType::RegionType region;<br>
region.SetIndex(index);<br>
region.SetSize(size);</div>
<div> ImageRegionIteratorWithIndexType iteratorWithIndex (distanceRidge, region);<br>
<br>
for (iteratorWithIndex.GoToBegin(); !iteratorWithIndex.IsAtEnd(); iteratorWithIndex++)<br>
{<br>
s1Index = iteratorWithIndex.GetIndex();<br>
r1SquaredK = (s1Index[0] - i) * (s1Index[0] - i);<br>
r1SquaredJK = r1SquaredK + (s1Index[1] - j) * (s1Index[1] - j);</div>
<div> if(r1SquaredJK <= rSquared)<br>
{<br>
r1Squared = r1SquaredJK + (s1Index[2] - k) * (s1Index[2] - k);<br>
if (r1Squared <= rSquared)<br>
{ <br>
s1 = iteratorWithIndex.Get();<br>
if (rSquared > s1)<br>
{<br>
iteratorWithIndex.Set(rSquared); <br>
}</div>
<div> } <br>
}<br>
}<br>
<br>
}</div>
<div> </div>
<div>so every approach I tried until now is very slow comparing to other implementations of the algorithm I want to do... would maybe spatial objects help me somehow? But I do not really understand how they work...</div>
<div>thanks for your help!</div>
<div> </div>
<div>greetings,</div>
<div>Johannes</div>
</div>
</div>
</div>
_____________________________________<br>
Powered by <a href="http://www.kitware.com" target="_blank">www.kitware.com</a><br>
<br>
Visit other Kitware open-source projects at<br>
<a href="http://www.kitware.com/opensource/opensource.html" target="_blank">http://www.kitware.com/opensource/opensource.html</a><br>
<br>
Kitware offers ITK Training Courses, for more information visit:<br>
<a href="http://www.kitware.com/products/protraining.php" target="_blank">http://www.kitware.com/products/protraining.php</a><br>
<br>
Please keep messages on-topic and check the ITK FAQ at:<br>
<a href="http://www.itk.org/Wiki/ITK_FAQ" target="_blank">http://www.itk.org/Wiki/ITK_FAQ</a><br>
<br>
Follow this link to subscribe/unsubscribe:<br>
<a href="http://public.kitware.com/mailman/listinfo/insight-users" target="_blank">http://public.kitware.com/mailman/listinfo/insight-users</a></blockquote>
</div>
</div>
</div>
<br>
_____________________________________<br>
Powered by <a href="http://www.kitware.com" target="_blank">www.kitware.com</a><br>
<br>
Visit other Kitware open-source projects at<br>
<a href="http://www.kitware.com/opensource/opensource.html" target="_blank">http://www.kitware.com/opensource/opensource.html</a><br>
<br>
Kitware offers ITK Training Courses, for more information visit:<br>
<a href="http://www.kitware.com/products/protraining.php" target="_blank">http://www.kitware.com/products/protraining.php</a><br>
<br>
Please keep messages on-topic and check the ITK FAQ at:<br>
<a href="http://www.itk.org/Wiki/ITK_FAQ" target="_blank">http://www.itk.org/Wiki/ITK_FAQ</a><br>
<br>
Follow this link to subscribe/unsubscribe:<br>
<a href="http://public.kitware.com/mailman/listinfo/insight-users" target="_blank">http://public.kitware.com/mailman/listinfo/insight-users</a><br>
</blockquote>
</div>
</div>
</div>
</div>
</blockquote>
</div>
</div>
</div>
</div>
</div>
</div></div></div>
</blockquote></div><br></div>