[Insight-users] PeriodicBoundayCondition problem...

devieill at irit.fr devieill at irit.fr
Tue Jan 19 09:06:49 EST 2010


Dear all,

I just do not quite understand why the PeriodicBoundaryCondition
does not initialize properly on this code example (see below).
Basically after compilation the program ends with a segmentation fault:

$ ./bin/test3
 Values for the image
  0   1   2   3
  4   5   6   7
  8   9  10  11
 12  13  14  15

 16  17  18  19
 20  21  22  23
 24  25  26  27
 28  29  30  31

 32  33  34  35
 36  37  38  39
 40  41  42  43
 44  45  46  47

 48  49  50  51
 52  53  54  55
 56  57  58  59
 60  61  62  63

 Values for the psf
  1   0   0
  0   0   0
  0   0   0

  0   0   0
  0   0   0
  0   0   0

  0   0   0
  0   0   0
  0   0   0

 index is (0,0,0) in a size of (3,3,3)
0
 index is (0,0,0) in a size of (3,3,3)
Segmentation fault (core dumped)



I use gcc :

$ gcc -v
Using built-in specs.
Target: x86_64-redhat-linux
Configured with: ../configure --prefix=/usr --mandir=/usr/share/man
--infodir=/usr/share/info
--with-bugurl=http://bugzilla.redhat.com/bugzilla --enable-bootstrap
--enable-shared --enable-threads=posix --enable-checking=release
--with-system-zlib --enable-__cxa_atexit --disable-libunwind-exceptions
--enable-gnu-unique-object
--enable-languages=c,c++,objc,obj-c++,java,fortran,ada
--enable-java-awt=gtk --disable-dssi --enable-plugin
--with-java-home=/usr/lib/jvm/java-1.5.0-gcj-1.5.0.0/jre
--enable-libgcj-multifile --enable-java-maintainer-mode
--with-ecj-jar=/usr/share/java/eclipse-ecj.jar --disable-libjava-multilib
--with-ppl --with-cloog --with-tune=generic --with-arch_32=i686
--build=x86_64-redhat-linux
Thread model: posix
gcc version 4.4.2 20091222 (Red Hat 4.4.2-20) (GCC)

and the itk 3.16 version.

Any ideas on this matter ?

Regards,
François

CODE
__________________________________________________________________

typedef double                                             PixelType;
typedef itk::Image< PixelType, 3 >                         WorkingImageType;
typedef itk::PeriodicBoundaryCondition<WorkingImageType>   BConditionType;
typedef itk::ConstNeighborhoodIterator< WorkingImageType,
					BConditionType >   ConstPBCNeighborhoodIteratorType;
typedef itk::NeighborhoodIterator< WorkingImageType,
				   BConditionType >        PBCNeighborhoodIteratorType;
typedef itk::ConstNeighborhoodIterator< WorkingImageType> 
ConstNeighborhoodIteratorType;
typedef itk::NeighborhoodIterator< WorkingImageType>      
NeighborhoodIteratorType;


typedef itk::ImageRegionConstIterator< WorkingImageType > 
ImageConstIteratorType;
typedef itk::ImageRegionIterator< WorkingImageType >       ImageIteratorType;

WorkingImageType::Pointer createAndInitialiseImage3D(
						     unsigned int dimx,
						     unsigned int dimy,
						     unsigned int dimz) {
  WorkingImageType::Pointer res = WorkingImageType::New();
  WorkingImageType::SizeType size;
  size[0] = dimx;   size[1] = dimy;   size[2] = dimz;
  WorkingImageType::IndexType start;
  start[0] = 0;  start[1] = 0;  start[2] = 0;
  WorkingImageType::RegionType region;
  region.SetSize(size);
  region.SetIndex(start);
  res->SetRegions(region);
  res->Allocate();
  return res;
}

int main( int argc, char *argv[]) {

  const unsigned int img_lg = 4;
  const unsigned int psf_lg = 3;

  WorkingImageType::Pointer src = createAndInitialiseImage3D( img_lg,
img_lg, img_lg);
  WorkingImageType::Pointer res = createAndInitialiseImage3D( img_lg,
img_lg, img_lg);
  WorkingImageType::Pointer psf = createAndInitialiseImage3D( psf_lg,
psf_lg, psf_lg);

  src->FillBuffer(5);
  psf->FillBuffer(0);
  res->FillBuffer(0);
  WorkingImageType::IndexType pixelIndex;
  int val = 0;
  for (unsigned int i = 0; i < img_lg; ++i) {
    pixelIndex[0] = i;
    for (unsigned int j = 0; j < img_lg; ++j) {
      pixelIndex[1] = j;
      for (unsigned int k = 0; k < img_lg; ++k) {
	pixelIndex[2] = k;
	src->SetPixel(pixelIndex, val++);
      }
    }
  }
  std::cout << " Values for the image " << std::endl;
  WorkingImageType::PixelType pixelValue;
  for (unsigned int i = 0; i < img_lg; ++i) {
    pixelIndex[0] = i;
    for (unsigned int j = 0; j < img_lg; ++j) {
      pixelIndex[1] = j;
      for (unsigned int k = 0; k < img_lg; ++k) {
	pixelIndex[2] = k;
	pixelValue = src->GetPixel(pixelIndex);
	std::cout << std::setw(3) << pixelValue << " ";
      }
      std::cout << std::endl;
    }
    std::cout << std::endl;
  }
  std::cout << " Values for the psf " << std::endl;
  for (unsigned int i = 0; i < psf_lg; ++i) {
    pixelIndex[0] = i;
    for (unsigned int j = 0; j < psf_lg; ++j) {
      pixelIndex[1] = j;
      for (unsigned int k = 0; k < psf_lg; ++k) {
	pixelIndex[2] = k;
	if ( i+j+k == 0 ) {
	  psf->SetPixel(pixelIndex, 1);
	}
	std::cout << std::setw(3) << psf->GetPixel(pixelIndex) << " ";
      }
      std::cout << std::endl;
    }
    std::cout << std::endl;
  }


  WorkingImageType::Pointer input_image  = src;
  WorkingImageType::Pointer input_filter = psf;
  WorkingImageType::Pointer output_image = res;

  NeighborhoodIteratorType::RadiusType radius;
  radius.Fill( 0 );
  for (unsigned int i = 0; i < WorkingImageType::ImageDimension; ++i)
{radius[i] = 1;}
  ConstPBCNeighborhoodIteratorType it_PBC;
  ConstNeighborhoodIteratorType it;
  NeighborhoodIteratorType it_out;

  ImageConstIteratorType it_psf(psf, psf->GetRequestedRegion() );

  it_PBC = ConstPBCNeighborhoodIteratorType(radius,
					input_image,
					output_image->GetLargestPossibleRegion());
  it = ConstNeighborhoodIteratorType(radius, 			     input_image,				    
output_image->GetLargestPossibleRegion());
 output_image->GetLargestPossibleRegion());

  WorkingImageType::IndexType regionIndex;
  WorkingImageType::SizeType  regionSize;

  it_PBC.NeedToUseBoundaryConditionOn();

  for ( it.GoToBegin(), it_PBC.GoToBegin(); !it.IsAtEnd(); ++it, ++it_PBC) {
    float accum = 0.0;
    it_psf.GoToBegin();
    regionIndex = it.GetIndex();
    regionSize  = it.GetSize();
    std::cout << " index is "
	      << "(" << regionIndex[0]
	      << "," << regionIndex[1]
	      << "," << regionIndex[2]
	      << ")"
	      << " in a size of "
	      << "(" << regionSize[0]
	      << "," << regionSize[1]
	      << "," << regionSize[2]
	      << ")"
	      << std::endl;

    std::cout << it.GetPixel(0) << " " << std::endl;

    regionIndex = it_PBC.GetIndex();
    regionSize  = it_PBC.GetSize();
    std::cout << " index is "
	      << "(" << regionIndex[0]
	      << "," << regionIndex[1]
	      << "," << regionIndex[2]
	      << ")"
	      << " in a size of "
	      << "(" << regionSize[0]
	      << "," << regionSize[1]
	      << "," << regionSize[2]
	      << ")"
	      << std::endl;


    std::cout << it_PBC.GetPixel(0) << " " << std::endl;

  }

  return 0;
}


More information about the Insight-users mailing list