1 /*=========================================================================
2  *
3  *  Copyright Insight Software Consortium
4  *
5  *  Licensed under the Apache License, Version 2.0 (the "License");
6  *  you may not use this file except in compliance with the License.
7  *  You may obtain a copy of the License at
8  *
9  *         http://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  *  Unless required by applicable law or agreed to in writing, software
12  *  distributed under the License is distributed on an "AS IS" BASIS,
13  *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  *  See the License for the specific language governing permissions and
15  *  limitations under the License.
16  *
17  *=========================================================================*/
18 
19 #include <iostream>
20 #include "itkTimeProbe.h"
21 #include "itkVectorImageToImageAdaptor.h"
22 #include "itkNthElementImageAdaptor.h"
23 #include "itkImageLinearIteratorWithIndex.h"
24 #include "itkShapedNeighborhoodIterator.h"
25 #include "itkImageRegionIteratorWithIndex.h"
26 #include "itkImageFileWriter.h"
27 #include "itkImageFileReader.h"
28 #include "itkTestingMacros.h"
29 #include "itkMath.h"
30 
31 
32 // This test tests:
33 // VectorImage
34 // -----------
35 //  - Iterators on the VectorImage,
36 //  - Timing tests comparing itk::VectorImage to similar itk::Image using
37 //    FixedArray and VariableLengthVector
38 //  - IO support for VectorImage.
39 //
40 
41 template< typename TPixel, unsigned int VDimension, typename TAdaptor, unsigned int VVectorLength >
testVectorImageAdaptor(typename TAdaptor::Pointer & vectorImageAdaptor,const typename itk::VectorImage<TPixel,VDimension>::Pointer & vectorImage,const typename itk::VectorImage<TPixel,VDimension>::RegionType & region,const unsigned int componentToExtract)42 bool testVectorImageAdaptor( typename TAdaptor::Pointer & vectorImageAdaptor,
43   const typename itk::VectorImage< TPixel, VDimension >::Pointer & vectorImage,
44   const typename itk::VectorImage< TPixel, VDimension >::RegionType & region,
45   const unsigned int componentToExtract )
46 {
47   std::cout << "---------------------------------------------------------------" << std::endl;
48   std::cout << "Testing " << typeid(TAdaptor).name();
49   std::cout << " to extract a component from the vector image" << std::endl;
50 
51   using PixelType = TPixel;
52   const unsigned int Dimension = VDimension;
53   const unsigned int VectorLength = VVectorLength;
54 
55   bool failed = false;
56 
57   using AdaptedImageType = itk::Image< PixelType , Dimension >;
58   typename AdaptedImageType::IndexType index;
59   index.Fill(10);
60   std::cout << "Before adaptor initialization, vectorImage->GetPixel("
61             << index << ")[" << componentToExtract
62             << "] = "
63             << vectorImage->GetPixel( index )[componentToExtract]
64             << std::endl;
65 
66 
67   using AdaptorType = TAdaptor;
68   vectorImageAdaptor->SetImage( vectorImage );
69   vectorImageAdaptor->Update();
70 
71   if(   (itk::Math::NotExactlyEquals(vectorImageAdaptor->GetPixel(index), vectorImage->GetPixel( index )[componentToExtract]))
72      || (itk::Math::NotExactlyEquals(vectorImage->GetPixel( index )[componentToExtract], componentToExtract) ))
73     {
74     std::cerr << "[FAILED]" << std::endl;
75     std::cerr << "vImageToImageAdaptor->GetPixel("
76               << index << ") = "
77               << vectorImageAdaptor->GetPixel(index)
78               << ", vectorImage->GetPixel("
79               << index << ")[" << componentToExtract
80               << "] = "
81               << vectorImage->GetPixel( index )[componentToExtract]
82               << std::endl;
83     failed = true;
84     }
85   else
86     {
87     std::cout << "[PASSED]" << std::endl;
88     }
89 
90     // Test adaptor with iterators:
91   itk::ImageRegionConstIteratorWithIndex< AdaptorType > adaptIt(vectorImageAdaptor, region);
92   adaptIt.GoToBegin();
93   bool itFailed = false;
94   using InternalPixelType = itk::FixedArray< PixelType, VectorLength >;
95   InternalPixelType f;
96   for( unsigned int i=0; i<VectorLength; i++ ) { f[i] = i; }
97   while (!adaptIt.IsAtEnd())
98     {
99     PixelType pixel = adaptIt.Get();
100     if (itk::Math::NotExactlyEquals(pixel, f[componentToExtract]))
101       {
102       itFailed = true;
103       std::cout << "adaptIt(" << adaptIt.GetIndex() << ") = " << adaptIt.Get()
104                 << ", but f[" << componentToExtract << "] = " << f[componentToExtract] << std::endl;
105       }
106     ++adaptIt;
107     }
108   if (itFailed)
109     {
110     std::cerr << "ImageRegionConstIteratorWithIndex on VectorImageAdaptor [FAILED]" << std::endl;
111     failed = true;
112     }
113   else
114     {
115     std::cout << "ImageRegionConstIteratorWithIndex on VectorImageAdaptor [PASSED]" << std::endl;
116     }
117 
118   return failed;
119 }
120 
121 template< typename TPixel, unsigned int VDimension >
testVectorImageBasicMethods()122 bool testVectorImageBasicMethods()
123 {
124   const unsigned int VectorLength = 3 * VDimension;
125 
126   using VectorImageType = itk::VectorImage< TPixel, VDimension >;
127 
128   std::cout << "Testing Get/SetPixel methods." << std::endl;
129 
130   typename VectorImageType::Pointer image = VectorImageType::New();
131   typename VectorImageType::SizeType  size = {{11, 9, 7}};
132   typename VectorImageType::RegionType region;
133   region.SetSize( size );
134   image->SetRegions( region );
135   image->SetNumberOfComponentsPerPixel( VectorLength );
136   image->Allocate();
137 
138   typename VectorImageType::PixelType f( VectorLength );
139   f.Fill(3.14f);
140 
141   image->FillBuffer( f );
142 
143   typename VectorImageType::IndexType idx = {{1,1,1}};
144 
145   typename VectorImageType::ConstPointer cimage(image);
146 
147   // test get methods
148 
149   TEST_EXPECT_EQUAL(f, image->GetPixel(idx));
150   TEST_EXPECT_EQUAL(f, cimage->GetPixel(idx));
151   TEST_EXPECT_EQUAL(f, (*image)[idx]);
152   TEST_EXPECT_EQUAL(f, (*cimage)[idx]);
153 
154 
155   // test get by reference methods
156 
157   typename VectorImageType::PixelType v( VectorLength );
158   v.Fill(2.22f);
159 
160   // note: this VectorImage method requires the compiler to perform return value
161   // optimization to work as expected
162   image->GetPixel(idx).Fill( 2.22f );
163   TEST_EXPECT_EQUAL(v, image->GetPixel(idx));
164   TEST_EXPECT_EQUAL(v, cimage->GetPixel(idx));
165   TEST_EXPECT_EQUAL(v, (*image)[idx]);
166   TEST_EXPECT_EQUAL(v, (*cimage)[idx]);
167 
168   v.Fill(2.23f);
169 
170   /*
171   // Cannot be used as an l-value
172   image->GetPixel(idx) = v;
173   TEST_EXPECT_EQUAL(v, image->GetPixel(idx));
174   TEST_EXPECT_EQUAL(v, cimage->GetPixel(idx));
175   TEST_EXPECT_EQUAL(v, (*image)[idx]);
176   TEST_EXPECT_EQUAL(v, (*cimage)[idx]);
177   */
178 
179   typename VectorImageType::PixelType temp = image->GetPixel(idx);
180   temp.Fill(2.24f);
181   v.Fill(2.24f);
182   TEST_EXPECT_EQUAL(v, image->GetPixel(idx));
183   TEST_EXPECT_EQUAL(v, cimage->GetPixel(idx));
184   TEST_EXPECT_EQUAL(v, (*image)[idx]);
185   TEST_EXPECT_EQUAL(v, (*cimage)[idx]);
186 
187   v.Fill(3.33f);
188 
189   (*image)[idx].Fill( 3.33f );
190   TEST_EXPECT_EQUAL(v, image->GetPixel(idx));
191   TEST_EXPECT_EQUAL(v, cimage->GetPixel(idx));
192   TEST_EXPECT_EQUAL(v, (*image)[idx]);
193   TEST_EXPECT_EQUAL(v, (*cimage)[idx]);
194 
195 
196   // test immutable access methods
197   v.Fill(3.33f);
198 
199   // comp error
200   // cimage->GetPixel(idx) = f;
201 
202   typename VectorImageType::PixelType temp2 = cimage->GetPixel(idx);
203   // the image actually get modified!!! :(
204   // The following line modifies the image and is considered a bug in
205   // the interface design that it works.
206   //temp2.Fill(3.44f);
207   TEST_EXPECT_EQUAL(v, image->GetPixel(idx));
208   TEST_EXPECT_EQUAL(v, cimage->GetPixel(idx));
209   TEST_EXPECT_EQUAL(v, (*image)[idx]);
210   TEST_EXPECT_EQUAL(v, (*cimage)[idx]);
211 
212   // test set method
213 
214   v.Fill(4.44f);
215 
216   image->SetPixel(idx, v);
217   TEST_EXPECT_EQUAL(v, image->GetPixel(idx));
218   TEST_EXPECT_EQUAL(v, cimage->GetPixel(idx));
219   TEST_EXPECT_EQUAL(v, (*image)[idx]);
220   TEST_EXPECT_EQUAL(v, (*cimage)[idx]);
221 
222   std::cout << "Testing Get/SetPixel methods [PASSED]" << std::endl;
223 
224   // test Graft method
225   typename VectorImageType::Pointer imageGraft = VectorImageType::New();
226   imageGraft->Graft(image);
227   TEST_EXPECT_EQUAL(image->GetPixelContainer(), imageGraft->GetPixelContainer());
228 
229   std::cout << "Testing Graft method [PASSED]" << std::endl;
230   return EXIT_SUCCESS;
231 }
232 
itkVectorImageTest(int,char * argv[])233 int itkVectorImageTest( int, char* argv[] )
234 {
235   bool failed = false;
236 
237   constexpr unsigned int Dimension = 3;
238   const unsigned int VectorLength = 2 * Dimension;
239   using PixelType = float;
240 
241   if ( testVectorImageBasicMethods<double, 3>() == EXIT_FAILURE )
242     {
243     std::cout << "Testing basic methods [FAILED]" << std::endl;
244     failed = true;
245     }
246 
247   {
248   // Test 1.
249   //
250   // Create an Image of VariableLengthVector, FixedArray, VectorImage of length 6 and compare
251   // times.
252   //
253   // Three images.. for crude timing analysis.
254 
255   using VariableLengthVectorImageType =
256       itk::Image< itk::VariableLengthVector< PixelType >, Dimension >;
257   using FixedArrayImageType = itk::Image< itk::FixedArray< PixelType, VectorLength >,
258                                       Dimension >;
259   using VectorImageType = itk::VectorImage< PixelType, Dimension >;
260 
261   // Using image of VariableLengthVector< PixelType >
262   {
263   using InternalPixelType = itk::VariableLengthVector< PixelType >;
264 
265   itk::TimeProbe clock;
266   clock.Start();
267 
268   VariableLengthVectorImageType::Pointer image = VariableLengthVectorImageType::New();
269   VariableLengthVectorImageType::IndexType start;
270   InternalPixelType f( VectorLength );
271   VariableLengthVectorImageType::SizeType  size;
272   for( unsigned int i=0; i<VectorLength; i++ ) { f[i] = i; }
273   start.Fill(0);
274   size.Fill(50);
275   VariableLengthVectorImageType::RegionType region;
276   region.SetSize( size );
277   region.SetIndex( start );
278   image->SetRegions( region );
279   image->Allocate();
280   image->FillBuffer( f );
281 
282   clock.Stop();
283   double timeTaken = clock.GetMean();
284   std::cout << "Allocating an image of itk::VariableLengthVector of length " <<  VectorLength
285           << " with image size " << size << " took " << timeTaken << " s." << std::endl;
286 
287     // Const iterator over the image...
288     {
289     clock.Start();
290     using IteratorType = itk::ImageRegionConstIterator< VariableLengthVectorImageType >;
291     IteratorType it( image, image->GetBufferedRegion() );
292     it.GoToBegin();
293     while( !it.IsAtEnd() )
294       {
295       it.Get();
296       ++it;
297       }
298     clock.Stop();
299     std::cout << "ConstIterator Get() over the entire image took : " <<
300       clock.GetMean() << " s." << std::endl;
301     }
302   }
303 
304 
305   // Using image of FixedArray< PixelType, VectorLength >
306   {
307   itk::TimeProbe clock;
308   clock.Start();
309 
310   using InternalPixelType = itk::FixedArray< PixelType, VectorLength >;
311 
312   FixedArrayImageType::Pointer image = FixedArrayImageType::New();
313   FixedArrayImageType::IndexType start;
314   InternalPixelType f;
315   FixedArrayImageType::SizeType  size;
316   for( unsigned int i=0; i<VectorLength; i++ ) { f[i] = i; }
317   start.Fill(0);
318   size.Fill(50);
319   FixedArrayImageType::RegionType region;
320   region.SetSize( size );
321   region.SetIndex( start );
322   image->SetRegions( region );
323   image->Allocate();
324   image->FillBuffer( f );
325 
326   clock.Stop();
327   double timeTaken = clock.GetMean();
328   std::cout << "Allocating an image of itk::FixedArray of length " <<  VectorLength
329           << " with image size " << size << " took " << timeTaken << " s." << std::endl;
330 
331     {
332     // Test and compare times with iterators
333     //
334     // First set some pixel
335     for( unsigned int i=0; i<VectorLength; i++ ) { f[i] = i*0.1; }
336     FixedArrayImageType::IndexType idx;
337     for (unsigned int i = 0; i < Dimension; i++)
338       {
339       idx[i] = 4;
340       }
341     idx[Dimension-1] = 12;
342     image->SetPixel( idx, f );
343     }
344 
345     {
346     clock.Start();
347     using IteratorType = itk::ImageRegionConstIterator< FixedArrayImageType >;
348     IteratorType it( image, image->GetBufferedRegion() );
349     it.GoToBegin();
350     while( !it.IsAtEnd() )
351       {
352       it.Get();
353       ++it;
354       }
355     clock.Stop();
356     std::cout << "ConstIterator Get() over the entire image took : " <<
357       clock.GetMean() << " s." << std::endl;
358     }
359   }
360 
361   // Using VectorImage< PixelType, Dimension >
362   {
363   itk::TimeProbe clock;
364   clock.Start();
365 
366   VectorImageType::Pointer vectorImage = VectorImageType::New();
367   VectorImageType::IndexType start;
368   itk::VariableLengthVector< PixelType > f( VectorLength );
369   VectorImageType::SizeType  size;
370   for( unsigned int i=0; i<VectorLength; i++ ) { f[i] = i; }
371   start.Fill(0);
372   size.Fill(50);
373   VectorImageType::RegionType region;
374   region.SetSize( size );
375   region.SetIndex( start );
376   vectorImage->SetVectorLength( VectorLength );
377   vectorImage->SetRegions( region );
378   vectorImage->Allocate();
379   vectorImage->FillBuffer( f );
380 
381   clock.Stop();
382   double timeTaken = clock.GetMean();
383   std::cout << "Allocating an image of itk::VectorImage with pixels length " <<  VectorLength
384      << " with image size " << size << " took " << timeTaken << " s." << std::endl;
385 
386 
387   // Iterator tests on the vector image.
388   //
389   // Const iterator over the vector image...
390   {
391     clock.Start();
392     using IteratorType = itk::ImageRegionConstIterator< VectorImageType >;
393     IteratorType it( vectorImage, vectorImage->GetBufferedRegion() );
394     it.GoToBegin();
395     while( !it.IsAtEnd() )
396       {
397       it.Get();
398       ++it;
399       }
400     clock.Stop();
401     std::cout << "ConstIterator Get() over the entire vectorImage took : " <<
402       clock.GetMean() << " s." << std::endl;
403   }
404 
405   const unsigned int componentToExtract = 2 * (Dimension -1);
406   using VectorImageToImageAdaptorType = itk::VectorImageToImageAdaptor< PixelType, Dimension >;
407   VectorImageToImageAdaptorType::Pointer vectorImageToImageAdaptor = VectorImageToImageAdaptorType::New();
408   vectorImageToImageAdaptor->SetExtractComponentIndex( componentToExtract );
409   if( vectorImageToImageAdaptor->GetExtractComponentIndex() != componentToExtract )
410     {
411     std::cerr << "[FAILED]" << std::endl;
412     failed = true;
413     }
414   bool adaptorTestResult;
415   adaptorTestResult = testVectorImageAdaptor< PixelType, Dimension, VectorImageToImageAdaptorType, VectorLength >( vectorImageToImageAdaptor,
416     vectorImage,
417     region,
418     componentToExtract );
419   if( adaptorTestResult )
420     {
421     failed = true;
422     }
423   using NthElementImageAdaptorType = itk::NthElementImageAdaptor< VectorImageType, PixelType >;
424   NthElementImageAdaptorType::Pointer nthElementImageAdaptor = NthElementImageAdaptorType::New();
425   nthElementImageAdaptor->SelectNthElement( componentToExtract );
426   adaptorTestResult = testVectorImageAdaptor< PixelType, Dimension, NthElementImageAdaptorType, VectorLength >( nthElementImageAdaptor,
427     vectorImage,
428     region,
429     componentToExtract );
430   if( adaptorTestResult )
431     {
432     failed = true;
433     }
434   }
435 
436   // Test with Region and Linear iterators...
437   {
438     // Create a  small image
439     VectorImageType::Pointer vectorImage = VectorImageType::New();
440     VectorImageType::IndexType start;
441     itk::VariableLengthVector< PixelType > f( VectorLength );
442     itk::VariableLengthVector< PixelType > ZeroPixel( VectorLength );
443     ZeroPixel.Fill( itk::NumericTraits< PixelType >::ZeroValue() );
444     for( unsigned int i=0; i<VectorLength; i++ ) { f[i] = i; }
445     start.Fill(0);
446     VectorImageType::SizeType  size;
447     size.Fill(11);
448     size[Dimension-1] = 5;
449     unsigned long midCtr = 1;
450     for (unsigned int i = 0; i < Dimension; i++) { midCtr *= size[i]; }
451     VectorImageType::RegionType region( start, size );
452     vectorImage->SetVectorLength( VectorLength );
453     vectorImage->SetRegions( region );
454     vectorImage->Allocate();
455     vectorImage->FillBuffer( ZeroPixel );
456 
457     start.Fill(3);
458     start[Dimension-1] = 2;
459     size.Fill(4);
460     size[Dimension-1] = 2;
461     VectorImageType::RegionType subRegion( start, size );
462     using ImageRegionIteratorType = itk::ImageRegionIterator< VectorImageType >;
463     ImageRegionIteratorType rit( vectorImage, subRegion );
464     rit.GoToBegin();
465 
466     while( !rit.IsAtEnd() )
467       {
468       rit.Set( f );
469       ++rit;
470       }
471 
472     using ConstIteratorType = itk::ImageRegionConstIterator< VectorImageType >;
473     ConstIteratorType cit( vectorImage, vectorImage->GetBufferedRegion() );
474     unsigned long ctr = 0;
475     cit.GoToBegin();
476     midCtr /= 2;
477     while( !cit.IsAtEnd() )
478       {
479       itk::VariableLengthVector< PixelType > value = cit.Get();
480       ++cit;
481       if( ctr == midCtr )
482         {
483         if( value != f )
484           {
485           std::cerr <<
486             "ImageRegionConstIteratorTest on VectorImage [FAILED]" << std::endl;
487           failed = true;
488           }
489         }
490       ++ctr;
491       }
492     std::cout << "ImageRegionConstIteratorTest on VectorImage [PASSED]" << std::endl;
493 
494 
495     {
496     // Test itkImageLinearIteratorWithIndex
497     using LinearConstIteratorType = itk::ImageLinearConstIteratorWithIndex< VectorImageType >;
498     using LinearIteratorType = itk::ImageLinearIteratorWithIndex< VectorImageType >;
499 
500     LinearConstIteratorType lcit( vectorImage, vectorImage->GetBufferedRegion() );
501     lcit.SetDirection( Dimension-1 );
502     lcit.GoToBegin();
503     itk::VariableLengthVector< PixelType > value;
504     while( !lcit.IsAtEnd() )
505       {
506       while( !lcit.IsAtEndOfLine() )
507         {
508         value = lcit.Get();
509         if( subRegion.IsInside( lcit.GetIndex() ) )
510           {
511           if( value!=f )
512             {
513             std::cerr <<
514               "ImageLinearConstIteratorWithIndex on VectorImage [FAILED]" << std::endl;
515             failed = true;
516             }
517           }
518         else
519           {
520           if( value!=ZeroPixel )
521             {
522             std::cerr <<
523               "ImageLinearConstIteratorWithIndex on VectorImage [FAILED]" << std::endl;
524             failed = true;
525             }
526           }
527         ++lcit;
528         }
529       lcit.NextLine();
530       }
531 
532     VectorImageType::IndexType idx;
533     idx.Fill(1);
534     LinearIteratorType lit( vectorImage, vectorImage->GetBufferedRegion() );
535     lit.SetIndex( idx );
536     lit.Set( f );
537 
538     lcit.SetIndex( idx );
539     value = lcit.Get();
540     if( value != f )
541       {
542       std::cerr <<
543         "ImageLinearConstIteratorWithIndex on VectorImage [FAILED]" << std::endl;
544       failed = true;
545       }
546 
547     std::cout << "ImageLinearConstIteratorWithIndex on VectorImage [PASSED]" << std::endl;
548     std::cout << "ImageLinearIteratorWithIndex on VectorImage [PASSED]" << std::endl;
549     }
550 
551   }
552 
553   }
554 
555 
556   // Test IO support.
557   {
558   // Create an image using itk::Vector
559   using VectorPixelType = itk::Vector< PixelType, VectorLength >;
560   using VectorImageType = itk::Image< itk::Vector< PixelType, VectorLength >,
561                                       Dimension >;
562   VectorImageType::Pointer image = VectorImageType::New();
563   VectorImageType::IndexType start;
564   start.Fill(0);
565   VectorImageType::SizeType size;
566   size.Fill(5);
567   VectorImageType::RegionType region( start, size );
568   image->SetRegions( region );
569   image->Allocate();
570 
571   using IteratorType = itk::ImageRegionIteratorWithIndex< VectorImageType >;
572   IteratorType it( image, region );
573   it.GoToBegin();
574 
575   while( !it.IsAtEnd() )
576     {
577     VectorPixelType f;
578     for (unsigned int i = 0; i < Dimension; i++)
579       {
580       f[i] = it.GetIndex()[i];
581       f[Dimension+i] = it.GetIndex()[i];
582       }
583     it.Set( f );
584     ++it;
585     }
586 
587   using WriterType = itk::ImageFileWriter< VectorImageType >;
588   WriterType::Pointer writer = WriterType::New();
589   writer->SetInput( image );
590   writer->SetFileName( argv[2] );
591   writer->Update();
592 
593   writer->SetFileName( argv[1] );
594   writer->Update();
595   }
596 
597   {
598   // Now read it as a itk::VectorImage.
599   using VectorImageType = itk::VectorImage< PixelType, Dimension >;
600   using ReaderType = itk::ImageFileReader< VectorImageType >;
601   ReaderType::Pointer reader = ReaderType::New();
602   reader->SetFileName( argv[1] );
603   reader->Update();
604 
605   VectorImageType::Pointer vectorImage = reader->GetOutput();
606 
607   using IteratorType = itk::ImageRegionConstIteratorWithIndex< VectorImageType >;
608   IteratorType cit( vectorImage, vectorImage->GetBufferedRegion() );
609   cit.GoToBegin();
610 
611   bool failed1 = false;
612   while( !cit.IsAtEnd() )
613     {
614     for (unsigned int i = 0; i < Dimension; i++)
615       {
616       if (itk::Math::NotExactlyEquals(cit.Get()[i], cit.GetIndex()[i]) ||
617           itk::Math::NotExactlyEquals(cit.Get()[i+Dimension], cit.GetIndex()[i]))
618         {
619         failed1 = true;
620         }
621       }
622     ++cit;
623     }
624 
625   if( failed1 )
626     {
627     std::cerr << "Read VectorImage [FAILED]" << std::endl;
628     failed = true;
629     }
630   else
631     {
632     std::cout << "Read VectorImage [PASSED]" << std::endl;
633     }
634 
635 
636   // Now write this out this VectorImage and read it again
637   using WriterType = itk::ImageFileWriter< VectorImageType >;
638   WriterType::Pointer writer = WriterType::New();
639   writer->SetInput( vectorImage );
640   writer->SetFileName(argv[1]);
641   writer->Update();
642   }
643 
644 
645   {
646   // Now read it as a itk::VectorImage.
647   using VectorImageType = itk::VectorImage< PixelType, Dimension >;
648   using ReaderType = itk::ImageFileReader< VectorImageType >;
649   ReaderType::Pointer reader = ReaderType::New();
650   reader->SetFileName( argv[1] );
651   reader->Update();
652 
653   VectorImageType::Pointer vectorImage = reader->GetOutput();
654 
655   using IteratorType = itk::ImageRegionConstIteratorWithIndex< VectorImageType >;
656   IteratorType cit( vectorImage, vectorImage->GetBufferedRegion() );
657   cit.GoToBegin();
658 
659   bool failed1 = false;
660   while( !cit.IsAtEnd() )
661     {
662     for (unsigned int i = 0; i < Dimension; i++)
663       {
664       if (itk::Math::NotExactlyEquals(cit.Get()[i], cit.GetIndex()[i]) ||
665           itk::Math::NotExactlyEquals(cit.Get()[i+Dimension], cit.GetIndex()[i]))
666         {
667         failed1 = true;
668         }
669       }
670     ++cit;
671     }
672 
673   if( failed1 )
674     {
675     std::cerr << "Write VectorImage [FAILED]" << std::endl;
676     failed = true;
677     }
678   else
679     {
680     std::cout << "Write VectorImage [PASSED]" << std::endl;
681     }
682 
683 
684     {
685     // Check support for Neighborhood Iterators
686     //
687     // 1. Test ConstNeighborhoodIterator
688     //
689     std::cout << "Testing ConstNeighborhoodIterator...." << std::endl;
690 
691     using ConstNeighborhoodIteratorType = itk::ConstNeighborhoodIterator<VectorImageType>;
692     ConstNeighborhoodIteratorType::RadiusType radius;
693     radius.Fill(1);
694 
695     ConstNeighborhoodIteratorType::RegionType region
696                             = vectorImage->GetBufferedRegion();
697     ConstNeighborhoodIteratorType::SizeType size;
698     size.Fill(4);
699     ConstNeighborhoodIteratorType::IndexType index;
700     index.Fill(1);
701     region.SetIndex(index);
702     region.SetSize( size );
703 
704     ConstNeighborhoodIteratorType cNit(radius, vectorImage, region);
705 
706     // Move Iterator to a point and see if it reads out the right value
707     //
708     unsigned int centerIndex = 1;
709     for (unsigned int i = 0; i < Dimension; i++)
710       { centerIndex *= (radius[i]*2+1); }
711     centerIndex /= 2;
712 
713     ConstNeighborhoodIteratorType::IndexType location;
714     for (unsigned int i = 0; i < Dimension; i++)
715       {
716       location[i] = i+1;
717       }
718     cNit.SetLocation( location );
719 
720     if( cNit.GetPixel(centerIndex) != vectorImage->GetPixel( location ) )
721       {
722       std::cerr << "  SetLocation [FAILED]" << std::endl;
723       failed=true;
724       }
725 
726     // Test GoToBegin()
727     cNit.GoToBegin();
728     if( cNit.GetPixel(centerIndex) != vectorImage->GetPixel( index ) )
729       {
730       std::cerr << "  GoToBegin [FAILED]" << std::endl;
731       failed=true;
732       }
733 
734     // Test GoToEnd()
735     cNit.GoToEnd();
736     --cNit;
737     ConstNeighborhoodIteratorType::IndexType endIndex;
738     endIndex.Fill(4);
739     if( cNit.GetPixel(centerIndex) != vectorImage->GetPixel( endIndex ) )
740       {
741       std::cerr << "  GoToEnd [FAILED]" << std::endl;
742       failed=true;
743       }
744 
745     // Test IsAtEnd()
746     if( !((++cNit).IsAtEnd()) )
747       {
748       std::cerr << "  IsAtEnd() [FAILED]" << std::endl;
749       failed = true;
750       }
751     cNit.GoToBegin();
752     unsigned int numPixelsTraversed = 1;
753     for (unsigned int i = 0; i < Dimension; i++)
754       {
755       numPixelsTraversed *= size[i];
756       }
757     while (! cNit.IsAtEnd())
758       {
759       ++cNit;
760       --numPixelsTraversed;
761       }
762     if( numPixelsTraversed ) { std::cerr << "  IsAtEnd() [FAILED]" << std::endl; }
763 
764     // Test operator-
765     --cNit;
766     ConstNeighborhoodIteratorType::OffsetType offset;
767     offset.Fill(1);
768     cNit -= offset;
769     itk::VariableLengthVector< PixelType > pixel = cNit.GetCenterPixel();
770     itk::VariableLengthVector< PixelType > correctAnswer( VectorLength );
771     correctAnswer.Fill( 3 );
772     if( pixel != correctAnswer )
773       {
774       std::cerr << "  operator- [FAILED]" << std::endl;
775       failed = true;
776       }
777 
778     // Test GetNeighborhood()
779     cNit.SetLocation( location );
780     ConstNeighborhoodIteratorType::NeighborhoodType
781                     neighborhood = cNit.GetNeighborhood();
782     //const unsigned int neighborhoodSize = neighborhood.Size();
783     //for( unsigned int i=0; i< neighborhoodSize; i++)
784     //  { std::cout << neighborhood[i] << std::endl; }
785     if( (itk::Math::NotExactlyEquals(neighborhood[0][0], 0)) || (itk::Math::NotExactlyEquals(neighborhood[0][2*Dimension-1], (Dimension-1))))
786       {
787       std::cerr << "  GetNeighborhood() on ConstNeighborhoodIterator [FAILED]" << std::endl;
788       failed = true;
789       }
790 
791     //
792     // 2. Test NeighborhoodIterator on VectorImage
793     //
794 
795     std::cout << "Testing NeighborhoodIterator..." << std::endl;
796 
797     using NeighborhoodIteratorType = itk::NeighborhoodIterator< VectorImageType >;
798     NeighborhoodIteratorType nit(radius, vectorImage, region);
799     nit.SetLocation( location );
800     itk::VariableLengthVector< PixelType > p( VectorLength );
801     p.Fill( 100.0 );
802     nit.SetNext( 1, 1, p );
803 
804     // Test SetNext()
805     NeighborhoodIteratorType::IndexType index1;
806     index1 = location; index1[1] = location[1] + 1;
807     nit.SetLocation( index1 );
808 
809     if( nit.GetCenterPixel() != p )
810       {
811       std::cerr << "  SetNext() [FAILED]" << std::endl;
812       failed = true;
813       }
814     for (unsigned i = 0; i < Dimension; i++)
815       {
816       p[i] = p[Dimension + i] = (float)index1[i];
817       }
818     nit.SetCenterPixel( p );
819     if( nit.GetCenterPixel() != p )
820       {
821       std::cerr << "  SetCenterPixel() [FAILED]" << std::endl;
822       failed = true;
823       }
824 
825     // Test SetNeighborhood() and GetPrevious()
826     nit.SetLocation( index1 );
827     nit.SetNeighborhood( neighborhood );
828     for (unsigned i = 0; i < Dimension; i++)
829       {
830       p[i] = p[Dimension + i] = i + 1;
831       }
832     p[Dimension-1] = p[2*Dimension - 1] = Dimension - 1;
833     if( nit.GetPrevious( Dimension-1, 1 ) != p )
834       {
835       std::cerr << "  SetNeighborhood() or GetPrevious() [FAILED]" << std::endl;
836       failed = true;
837       }
838 
839     if (Dimension == 3)
840       {
841 
842       //
843       // 3. Testing ConstShapedNeighborhoodIterator on VectorImage
844       //
845 
846       // Go back to original image, where pixel values tell us the indices
847       std::cout << "Testing ConstShapedNeighborhoodIterator on VectorImage..."
848                                                                     << std::endl;
849       reader->SetFileName( "dummy.nrrd");
850       reader->SetFileName( argv[1] );
851       reader->Update();
852       vectorImage = reader->GetOutput();
853 
854       using ConstShapedNeighborhoodIteratorType = itk::ConstShapedNeighborhoodIterator<VectorImageType>;
855       ConstShapedNeighborhoodIteratorType cSnit( radius, vectorImage, region );
856       cSnit.SetLocation( location );
857       ConstShapedNeighborhoodIteratorType::OffsetType offset1;
858       offset1[0] = 0; offset1[1] = 0; offset1[2] = 0;
859       cSnit.ActivateOffset( offset1 ); //activate the center
860       // activate the top plane
861       offset1[0] = 0; offset1[1] = 0; offset1[2] = -1;
862       cSnit.ActivateOffset( offset1 );
863       offset1[0] = 0; offset1[1] = 1; offset1[2] = -1;
864       cSnit.ActivateOffset( offset1 );
865       offset1[0] = 0; offset1[1] = -1; offset1[2] = -1;
866       cSnit.ActivateOffset( offset1 );
867       offset1[0] = 1; offset1[1] = 0; offset1[2] = -1;
868       cSnit.ActivateOffset( offset1 );
869       offset1[0] = 1; offset1[1] = 1; offset1[2] = -1;
870       cSnit.ActivateOffset( offset1 );
871       offset1[0] = 1; offset1[1] = -1; offset1[2] = -1;
872       cSnit.ActivateOffset( offset1 );
873       offset1[0] = -1; offset1[1] = 0; offset1[2] = -1;
874       cSnit.ActivateOffset( offset1 );
875       offset1[0] = -1; offset1[1] = 1; offset1[2] = -1;
876       cSnit.ActivateOffset( offset1 );
877       offset1[0] = -1; offset1[1] = -1; offset1[2] = -1;
878       cSnit.ActivateOffset( offset1 );
879 
880       ConstShapedNeighborhoodIteratorType::IndexListType l
881                                     = cSnit.GetActiveIndexList();
882       ConstShapedNeighborhoodIteratorType::IndexListType::const_iterator
883                                                           ali = l.begin();
884       while (ali != l.end())
885         {
886         std::cout << *ali << " ";
887         ++ali;
888         }
889       std::cout << std::endl;
890 
891       ConstShapedNeighborhoodIteratorType::ConstIterator ci = cSnit.Begin();
892       while (! ci.IsAtEnd())
893         {
894         ConstShapedNeighborhoodIteratorType::OffsetType offset2 = ci.GetNeighborhoodOffset();
895         if( (offset2[0] == -1) && (offset2[1]== -1) && (offset2[2]== -1) )
896           {
897           if( ci.GetNeighborhoodIndex() != 0 )
898             {
899             failed = true;
900             std::cerr << "GetNeighborhoodOffset() on ConstShapedNeighborhoodIterato [FAILED]"
901                                                                                 << std::endl;
902             }
903           if( (itk::Math::NotExactlyEquals(ci.Get()[0], 0)) || (itk::Math::NotExactlyEquals(ci.Get()[1], 1)) || (itk::Math::NotExactlyEquals(ci.Get()[2] , 2)) )
904             {
905             failed=true;
906             std::cerr
907               << "ConstShapedNeighborhoodIterator returned incorrect index [FAILED]"
908                                                                           << std::endl;
909             }
910           }
911         ci++;
912         }
913 
914       //
915       // 4. Test ShapedNeighborhoodIterator
916       //
917       using ShapedNeighborhoodIteratorType = itk::ShapedNeighborhoodIterator<VectorImageType>;
918       ShapedNeighborhoodIteratorType sNit( radius, vectorImage, region );
919 
920       offset1[0] = 0; offset1[1] = 0; offset1[2] = 0;
921       sNit.ActivateOffset( offset1 ); //activate the center
922       // activate the top plane
923       offset1[0] = 0; offset1[1] = 0; offset1[2] = -1;
924       sNit.ActivateOffset( offset1 );
925       offset1[0] = 0; offset1[1] = 1; offset1[2] = -1;
926       sNit.ActivateOffset( offset1 );
927       offset1[0] = 0; offset1[1] = -1; offset1[2] = -1;
928       sNit.ActivateOffset( offset1 );
929       offset1[0] = 1; offset1[1] = 0; offset1[2] = -1;
930       sNit.ActivateOffset( offset1 );
931       offset1[0] = 1; offset1[1] = 1; offset1[2] = -1;
932       sNit.ActivateOffset( offset1 );
933       offset1[0] = 1; offset1[1] = -1; offset1[2] = -1;
934       sNit.ActivateOffset( offset1 );
935       offset1[0] = -1; offset1[1] = 0; offset1[2] = -1;
936       sNit.ActivateOffset( offset1 );
937       offset1[0] = -1; offset1[1] = 1; offset1[2] = -1;
938       sNit.ActivateOffset( offset1 );
939       offset1[0] = -1; offset1[1] = -1; offset1[2] = -1;
940       sNit.ActivateOffset( offset1 );
941 
942       sNit.SetLocation( location );
943       ShapedNeighborhoodIteratorType::Iterator shit = sNit.Begin();
944       p[0] = p[3] = 10; p[1] = p[4] = 20; p[2] = p[5] = 30;
945       shit.Set( p );
946       index[0]=location[0]-1; index[1]=location[1]-1; index[2]=location[2]-1;
947       cNit.SetLocation( index );
948       if( cNit.GetCenterPixel() != p )
949         {
950         std::cerr << "ShapedNeighborhoodIterator Set() [FAILED]" << std::endl;
951         failed=true;
952         }
953       }
954 
955 
956     }  // End Testing Neighborhood Iterators on VectorImage
957 
958   }
959 
960   if( failed )
961     {
962     return EXIT_FAILURE;
963     }
964 
965   return EXIT_SUCCESS;
966 }
967