<div dir="ltr"><div><div><div><div><div><div>Hi Matt, Hi Kent,<br><br>Thank you for the suggestions!<br><br></div>I have tried the NeighborhoodIterator and then I did a little test comparing that with directly calling GetPixel() for each of the 27 neighbors.<br>
<br></div>Basically, I randomly pick 10,000,000 points inside the "inner-region" of the an 3D volume. For each point, take average of all its 27 neighbors, and set its value to the average value.<br><br></div>Since the random points are all in [1 ,size[i]-2], range, accessing their +/-1 neighbors are safe. So I did not check boundary when accessing.<br>
<br></div>I just used linux command "time" to time the code, as a rough test.<br><br>Both takes about 5.5seconds (Release build, again Release version of ITK 4.3.0). There is not dramatic difference... <br><br></div>
<div>When using NeighborhoodIterator, i'm not sure if the boundary checking is automatically turned on. I guess if that is by default on, then there is not much speed gain with it.<br></div><div><br>I paste the code below. Could I know if there's something I didn't do right? Thank you!<br>
</div><div><br><br>The code are:<br><br></div>1. using NeighborhoodIterator<br><br>#include "itkImage.h"<br>#include "itkImageFileReader.h"<br>#include "itkImageFileWriter.h"<br><br>#include "itkNeighborhoodIterator.h"<br>
<br>#include <vnl/vnl_random.h><br><br>int main( int argc, char ** argv )<br>{<br> if ( argc < 3 )<br> {<br> std::cerr << "Usage: " << std::endl;<br> std::cerr << argv[0]<br>
<< " inputImage numberOfPoints\n";<br> return -1;<br> }<br><br> typedef float PixelType;<br> typedef itk::Image< PixelType, 3 > ImageType;<br>
<br> typedef itk::ImageFileReader< ImageType > ImageReaderType;<br> typename ImageReaderType::Pointer reader = ImageReaderType::New();<br> reader->SetFileName(argv[1]);<br> reader->Update();<br> typename ImageType::Pointer image = reader->GetOutput();<br>
<br> ImageType::SizeType size = image->GetLargestPossibleRegion().GetSize();<br> <br> typedef itk::NeighborhoodIterator< ImageType > NeighborhoodIteratorType;<br><br> NeighborhoodIteratorType::RadiusType radius;<br>
radius.Fill(1);<br> NeighborhoodIteratorType it(radius, image, image->GetLargestPossibleRegion());<br><br> ImageType::IndexType index;<br> long n = ::atol(argv[2]);<br> std::cout<<n<<std::endl;<br><br>
vnl_random rg;<br><br> PixelType s = static_cast<PixelType>(it.Size());<br> std::cout<<s<<std::endl;<br><br> for (long ip = 0; ip < n; ++ip)<br> {<br> index[0] = rg.lrand32(1, size[0] - 1);<br>
index[1] = rg.lrand32(1, size[1] - 1);<br> index[2] = rg.lrand32(1, size[2] - 1);<br><br> it.SetLocation(index);<br> //std::cout<<index<<std::endl;<br><br> PixelType sum = it.GetCenterPixel();<br>
for (unsigned i = 0; i < it.Size(); i++)<br> {<br> sum += it.GetPixel(i);<br> }<br><br> it.SetCenterPixel( sum/s );<br><br> //std::cout<<sum/s<<std::endl;<br> }<br>
<br> typedef itk::ImageFileWriter< ImageType > WriterType;<br><br> typename WriterType::Pointer writer = WriterType::New();<br> writer->SetFileName( "o.nrrd" );<br> writer->SetInput(image);<br> writer->UseCompressionOff();<br>
writer->Update();<br><br> return 0;<br>}<br><br><br><br><br></div>2. using direct access<br>#include "itkImage.h"<br>#include "itkImageFileReader.h"<br>#include "itkImageFileWriter.h"<br>
<br>#include <vnl/vnl_random.h><br><br>int main( int argc, char ** argv )<br>{<br> if ( argc < 3 )<br> {<br> std::cerr << "Usage: " << std::endl;<br> std::cerr << argv[0]<br>
<< " inputImage numberOfPoints\n";<br> return -1;<br> }<br><br> typedef float PixelType;<br> typedef itk::Image< PixelType, 3 > ImageType;<br>
<br> typedef itk::ImageFileReader< ImageType > ImageReaderType;<br> typename ImageReaderType::Pointer reader = ImageReaderType::New();<br> reader->SetFileName(argv[1]);<br> reader->Update();<br> typename ImageType::Pointer image = reader->GetOutput();<br>
<br> ImageType::SizeType size = image->GetLargestPossibleRegion().GetSize();<br> <br> ImageType::IndexType index;<br> ImageType::IndexType nbhdindex;<br><br> long n = ::atol(argv[2]);<br> std::cout<<n<<std::endl;<br>
<br> vnl_random rg;<br><br> PixelType s = 27.0;<br><br> for (long ip = 0; ip < n; ++ip)<br> {<br> PixelType sum = 0.0;<br><br> int cx = rg.lrand32(1, size[0] - 2);<br> int cy = rg.lrand32(1, size[1] - 2);<br>
int cz = rg.lrand32(1, size[2] - 2);<br><br> index[0] = cx;<br> index[1] = cy;<br> index[2] = cz;<br><br> for (int iz = cz - 1; iz <= cz + 1; ++iz)<br> {<br> nbhdindex[2] = iz;<br>
for (int iy = cy - 1; iy <= cy + 1; ++iy)<br> {<br> nbhdindex[1] = iy;<br> for (int ix = cx - 1; ix <= cx + 1; ++ix)<br> {<br> nbhdindex[0] = ix;<br>
<br> sum += image->GetPixel(nbhdindex);<br> }<br> }<br> }<br><br> image->SetPixel(index, sum/s);<br> }<br><br> typedef itk::ImageFileWriter< ImageType > WriterType;<br>
<br> typename WriterType::Pointer writer = WriterType::New();<br> writer->SetFileName( "p.nrrd" );<br> writer->SetInput(image);<br> writer->UseCompressionOff();<br> writer->Update();<br><br> return 0;<br>
}<br><br><br><br><br><br><br><br><div><div><br><br></div></div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Mon, May 13, 2013 at 10:02 AM, Williams, Norman K <span dir="ltr"><<a href="mailto:norman-k-williams@uiowa.edu" target="_blank">norman-k-williams@uiowa.edu</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">I would think an itk::NeighborhoodIterator with a radius of [1,1,1] would<br>
do. Are there any requirements other than direct adjacency to the<br>
neighborhood you want to traverse?<br>
<br>
I haven't looked closely at the code but I think it's about as well<br>
optimized as it could be.<br>
--<br>
Kent Williams <a href="mailto:norman-k-williams@uiowa.edu">norman-k-williams@uiowa.edu</a><br>
<div><div class="h5"><br>
<br>
<br>
<br>
<br>
<br>
On 5/12/13 8:43 PM, "Gao, Yi" <<a href="mailto:gaoyi.cn@gmail.com">gaoyi.cn@gmail.com</a>> wrote:<br>
<br>
>Dear all,<br>
><br>
><br>
>Could I have some suggestions on ways to fast access the nerighbors of a<br>
>set of random voxels? The scenario is as follows:<br>
><br>
><br>
>I have to access some voxels in the volume, these voxels are quite<br>
>randomly distributed. For each of them I know their (x, y, z) indexes.<br>
><br>
>At each of these voxel, I need to access all its 26 neighbors. I want<br>
>this step to be faster.<br>
><br>
><br>
>So i'm hoping to have something like:<br>
><br>
><br>
>1. for each (x, y, z), I get the pointer to that voxel: ptr.<br>
><br>
>2. For each of its neighbors, I use things like ptr[1], ptr[nx], ptr[1+<br>
>nx*ny] to access (x+1, y, z), (x, y+1, z) and (x+1, y, z+1)....<br>
><br>
><br>
>Could I know if there is some good way to do that?<br>
><br>
>Thank you in advance for any hint!<br>
><br>
>Best,<br>
>yi<br>
><br>
<br>
<br>
<br>
</div></div>________________________________<br>
Notice: This UI Health Care e-mail (including attachments) is covered by the Electronic Communications Privacy Act, 18 U.S.C. 2510-2521, is confidential and may be legally privileged. If you are not the intended recipient, you are hereby notified that any retention, dissemination, distribution, or copying of this communication is strictly prohibited. Please reply to the sender that you have received the message in error, then delete it. Thank you.<br>
________________________________<br>
</blockquote></div><br></div>