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