C++ Pointer::SetIdentity方法代码示例

	void initialiseFilters() {
		// resamplers
		transform = TransformType::New();
	  transform->SetIdentity();
    volumeInterpolator = VolumeInterpolatorType::New();
		maskVolumeInterpolator = MaskVolumeInterpolatorType::New();
		resampler = ResamplerType::New();
    resampler->SetInput( originalImage );
		resampler->SetInterpolator( volumeInterpolator );
		resampler->SetOutputSpacing( resamplerSpacing );
		resampler->SetSize( resamplerSize );
		resampler->SetTransform( transform );
		resampler->SetDefaultPixelValue( 127 );
		maskResampler = MaskResamplerType::New();
		maskResampler->SetInput( originalMask );
		maskResampler->SetInterpolator( maskVolumeInterpolator );
		maskResampler->SetOutputSpacing( resamplerSpacing );
		maskResampler->SetSize( resamplerSize );
		maskResampler->SetTransform( transform );
		
		// extract image filters
    sliceExtractor = SliceExtractorType::New();
    sliceExtractor->SetInput( resampler->GetOutput() );
		maskSliceExtractor = MaskSliceExtractorType::New();
    maskSliceExtractor->SetInput( maskResampler->GetOutput() );
    
    // masks
    for(unsigned int i=0; i<resamplerSize[2]; i++) {
  		masks2D.push_back( MaskType2D::New() );
    }		
	}

//.........这里部分代码省略.........
  roiFilter->SetRegionOfInterest(region);
  roiFilter->ReleaseDataFlagOn();
  roiFilter->ReleaseDataBeforeUpdateFlagOn();

  typedef itk::ConstantPadImageFilter<CharImageType, CharImageType> PadFilterType;

  PadFilterType::Pointer padFilter = PadFilterType::New();
  const PadFilterType::SizeValueType pad[3] = { 10, 10, 10 };

  padFilter->SetInput(roiFilter->GetOutput());
  padFilter->SetConstant(0);
  padFilter->SetPadLowerBound(pad);
  padFilter->SetPadUpperBound(pad);
  padFilter->ReleaseDataFlagOn();
  padFilter->ReleaseDataBeforeUpdateFlagOn();
  padFilter->Update();

  CharImageType::Pointer roiImage = padFilter->GetOutput();

  roiImage->DisconnectPipeline();
  roiFilter = nullptr;
  padFilter = nullptr;

  // Correct origin of real geometry (changed by cropping and padding)

  typedef Geometry3D::TransformType TransformType;

  TransformType::Pointer transform = TransformType::New();
  TransformType::OutputVectorType translation;

  for (unsigned int dim = 0; dim < 3; ++dim)
    translation[dim] = (int)minIndex[dim] - (int)pad[dim];

  transform->SetIdentity();
  transform->Translate(translation);
  geometry->Compose(transform, true);

  ProgressBar::GetInstance()->Progress(1);

  // Median

  MITK_INFO << "Median...";

  typedef itk::BinaryMedianImageFilter<CharImageType, CharImageType> MedianFilterType;

  MedianFilterType::Pointer medianFilter = MedianFilterType::New();
  CharImageType::SizeType radius = { 0 };

  medianFilter->SetRadius(radius);
  medianFilter->SetBackgroundValue(0);
  medianFilter->SetForegroundValue(1);
  medianFilter->SetInput(roiImage);
  medianFilter->ReleaseDataFlagOn();
  medianFilter->ReleaseDataBeforeUpdateFlagOn();
  medianFilter->Update();

  ProgressBar::GetInstance()->Progress(1);

  // Intelligent closing

  MITK_INFO << "Intelligent closing...";

  unsigned int surfaceRatio = (unsigned int)((1.0f - closing) * 100.0f);

  typedef itk::IntelligentBinaryClosingFilter<CharImageType, ShortImageType> ClosingFilterType;

/**
 * @brief      3D resample data to new grid size
 *
 * @param  M   Incoming data
 * @param  f   Resampling factor in all 3 dimensions
 * @param  im  Interpolation method (LINEAR|BSPLINE)
 *
 * @return     Resampled data
 */
template<class T> static Matrix<T> 
resample (const Matrix<T>& M, const Matrix<double>& f, const InterpMethod& im) {
	

	Matrix <T> res = M;
	
#ifdef HAVE_INSIGHT
	
	typedef typename itk::OrientedImage< T, 3 > InputImageType;
	typedef typename itk::OrientedImage< T, 3 > OutputImageType;
	typedef typename itk::IdentityTransform< double, 3 > TransformType;
	typedef typename itk::LinearInterpolateImageFunction< InputImageType, double > InterpolatorType;
	typedef typename itk::ResampleImageFilter< InputImageType, InputImageType > ResampleFilterType;
	
	typename InterpolatorType::Pointer linterp = InterpolatorType::New();
	
	TransformType::Pointer trafo = TransformType::New();
	trafo->SetIdentity();
	
	typename InputImageType::SpacingType space;
	space[0] = 1.0/f[0];
	space[1] = 1.0/f[1];
	space[2] = 1.0/f[2];
	
	typedef typename InputImageType::SizeType::SizeValueType SizeValueType;
	typename InputImageType::SizeType size; 
	size[0] = static_cast<SizeValueType>(res.Dim(0));
	size[1] = static_cast<SizeValueType>(res.Dim(1));
	size[2] = static_cast<SizeValueType>(res.Dim(2));
	
	typename itk::OrientedImage< T, 3 >::Pointer input = itk::OrientedImage< T, 3 >::New();
	typename itk::OrientedImage< T, 3 >::Pointer output = itk::OrientedImage< T, 3 >::New();
	
	typename itk::Image< T, 3 >::IndexType ipos;
	ipos[0] = 0; ipos[1] = 0; ipos[2] = 0;
	typename itk::Image< T, 3 >::IndexType opos;
	opos[0] = 0; opos[1] = 0; opos[2] = 0;
	
	typename itk::Image< T, 3 >::RegionType ireg;
	ireg.SetSize(size);
	ireg.SetIndex(ipos);
	input->SetRegions(ireg);
	input->Allocate();
	
	typename itk::Image< T, 3 >::RegionType oreg;
	oreg.SetSize(size);
	ireg.SetIndex(opos);
	output->SetRegions(oreg);
	output->Allocate();
	
	for (size_t z = 0; z < res.Dim(2); z++)
		for (size_t y = 0; y < res.Dim(1); y++)
			for (size_t x = 0; x < res.Dim(0); x++) {
				ipos[0] = x; ipos[1] = y; ipos[2] = z;
				input->SetPixel (ipos, res.At(x,y,z));
			}
	
	typename ResampleFilterType::Pointer rs = ResampleFilterType::New();
	rs->SetInput( input );
	rs->SetTransform( trafo );
	rs->SetInterpolator( linterp );
	rs->SetOutputOrigin ( input->GetOrigin());
	rs->SetOutputSpacing ( space );
	rs->SetOutputDirection ( input->GetDirection());
	rs->SetSize ( size );
	rs->Update ();
	
	output = rs->GetOutput();
	
	res = Matrix<T> (res.Dim(0)*f[0], res.Dim(1)*f[1], res.Dim(2)*f[2]);
	res.Res(0) = res.Res(0)/f[0];
	res.Res(1) = res.Dim(1)/f[1];
	res.Res(2) = res.Dim(2)/f[2];
	
	for (size_t z = 0; z < res.Dim(2); z++)
		for (size_t y = 0; y < res.Dim(1); y++)
			for (size_t x = 0; x < res.Dim(0); x++) {
				opos[0] = x; opos[1] = y; opos[2] = z;
				res.At(x,y,z) = output->GetPixel (opos);
			}
	
#else 
	
	printf ("ITK ERROR - Resampling not performed without ITK!\n");
	
#endif
	
	return res;
	
}

//.........这里部分代码省略.........
	catch( itk::ExceptionObject & err ) 
	  { 
	  std::cerr << "ERROR: ExceptionObject caught !" << std::endl; 
	  std::cerr << err << std::endl; 
	  return EXIT_FAILURE;
	  } 
      }
    }
    
    if (!image) {
	std::cerr << argv[0] << ": input output" << std::endl;
	exit(1);
    }
  

  typedef itk::SigmoidImageFilter< InputImageType, FilterImageType > SigmoidCasterType;
  SigmoidCasterType::Pointer sigmoidcaster = SigmoidCasterType::New();
  
  sigmoidcaster->SetInput( image );
  sigmoidcaster->SetOutputMaximum( 4000 );
  sigmoidcaster->SetOutputMinimum( 1000 );

  
  typedef itk::AccumulateImageFilter< FilterImageType, FilterImageType > AccumulateFilter;
  AccumulateFilter::Pointer accumulator = AccumulateFilter::New();
  accumulator->SetAccumulateDimension(1);
  accumulator->SetInput( sigmoidcaster->GetOutput() );

  typedef itk::ExtractImageFilter< FilterImageType, OutFilterImageType > ExtractFilter;
  ExtractFilter::Pointer extractor = ExtractFilter::New();
  extractor->SetInput( accumulator->GetOutput() );
  FilterImageType::Pointer accuOut = accumulator->GetOutput();
  accuOut->UpdateOutputInformation();
  FilterImageType::RegionType extractRegion = accuOut->GetLargestPossibleRegion();
  
  extractRegion.SetSize(1,0);
  
  extractor->SetExtractionRegion( extractRegion );

  typedef itk::ResampleImageFilter<OutFilterImageType, OutFilterImageType > ResampleFilter;
  ResampleFilter::Pointer resampler = ResampleFilter::New();
  resampler->SetInput( extractor->GetOutput() );
  
  typedef itk::BSplineInterpolateImageFunction< OutFilterImageType > InterpolatorType;
  InterpolatorType::Pointer interpolator = InterpolatorType::New();
  interpolator->SetSplineOrder(3);
  
  resampler->SetInterpolator( interpolator );
  OutFilterImageType::Pointer exOut = extractor->GetOutput();
  exOut->UpdateOutputInformation();
  
  typedef itk::CenteredRigid2DTransform< double > TransformType;
  TransformType::Pointer transform = TransformType::New();
  transform->SetIdentity();
  OutFilterImageType::PointType exOutCenter = exOut->GetOrigin();
  exOutCenter[0] += (exOut->GetLargestPossibleRegion().GetSize()[0]-1) * exOut->GetSpacing()[0] *.5;
  exOutCenter[1] += (exOut->GetLargestPossibleRegion().GetSize()[1]-1) * exOut->GetSpacing()[1] *.5;
  transform->SetCenter( exOutCenter );
  transform->SetAngleInDegrees( 180 );
  resampler->SetTransform( transform );
  resampler->SetOutputParametersFromImage( exOut );

  OutFilterImageType::SpacingType resampleSpacing = exOut->GetSpacing();
  resampleSpacing.Fill( std::min( resampleSpacing[0], resampleSpacing[1] ) );
  OutFilterImageType::SizeType resampleSize;
  resampleSize[0] = exOut->GetLargestPossibleRegion().GetSize()[0] * exOut->GetSpacing()[0] / resampleSpacing[0];
  resampleSize[1] = exOut->GetLargestPossibleRegion().GetSize()[1] * exOut->GetSpacing()[1] / resampleSpacing[1];
  resampler->SetSize( resampleSize );
  resampler->SetOutputSpacing( resampleSpacing );
  
  OutFilterImageType::Pointer result = resampler->GetOutput();
  
  sigmoidcaster->SetBeta( -500 );
  sigmoidcaster->SetAlpha( 5 );
  result->Update();

  int outDicomIndex = 0;
  itk::EncapsulateMetaData( dict, "0008|0008", string("DERIVED\\SECONDARY\\AXIAL"));
  
  boost::filesystem::path outpath = output_dir;
  outpath = outpath / "IM%06d";
  
  std::vector< itk::MetaDataDictionary* > dictArray;
  dictArray.push_back(&dict);
  
  itkBasic::writeDicomSeries( itkBasic::ImageRescale(itkBasic::ImageSharp(result, 0.5), -1000, 4000), outpath.string(), &dictArray, outDicomIndex);
//  itkBasic::ImageSave( itkBasic::ImageSharp(result, 0.5), boost::str( boost::format("%s.%s.png") % output_name % "lung" ), 1, 0); // Auto Level

  sigmoidcaster->SetBeta( 1000 );
  sigmoidcaster->SetAlpha( 300 );
  result->Update();
  itkBasic::writeDicomSeries( itkBasic::ImageRescale(itkBasic::ImageSharp(result, 0.5), -1000, 4000), outpath.string(), &dictArray, outDicomIndex);
//  itkBasic::ImageSave( itkBasic::ImageSharp(result, 0.5), boost::str( boost::format("%s.%s.png") % output_name % "bone" ), 1, 0); // Auto Level
  
  sigmoidcaster->SetBeta( 0 );
  sigmoidcaster->SetAlpha( 2000 );
  result->Update();
  itkBasic::writeDicomSeries( itkBasic::ImageRescale(itkBasic::ImageSharp(result, 0.5), -1000, 4000), outpath.string(), &dictArray, outDicomIndex);
//  itkBasic::ImageSave( itkBasic::ImageSharp(result, 0.5), boost::str( boost::format("%s.%s.png") % output_name % "normal" ), 1, 0); // Auto Level
}

bool mitk::NavigationDataLandmarkTransformFilter::FindCorrespondentLandmarks(LandmarkPointContainer& sources, const LandmarkPointContainer& targets) const
{
  if (sources.size() < 6 || targets.size() < 6)
    return false;
  //throw std::invalid_argument("ICP correspondence finding needs at least 6 landmarks");

  /* lots of type definitions */
  typedef itk::PointSet<mitk::ScalarType, 3> PointSetType;
  //typedef itk::BoundingBox<PointSetType::PointIdentifier, PointSetType::PointDimension> BoundingBoxType;

  typedef itk::EuclideanDistancePointMetric< PointSetType, PointSetType> MetricType;
  //typedef MetricType::TransformType TransformBaseType;
  //typedef MetricType::TransformType::ParametersType ParametersType;
  //typedef TransformBaseType::JacobianType JacobianType;
  //typedef itk::Euler3DTransform< double > TransformType;
  typedef itk::VersorRigid3DTransform< double > TransformType;
  typedef TransformType ParametersType;
  typedef itk::PointSetToPointSetRegistrationMethod< PointSetType, PointSetType > RegistrationType;

  /* copy landmarks to itk pointsets for registration */
  PointSetType::Pointer sourcePointSet = PointSetType::New();
  unsigned int i = 0;
  for (LandmarkPointContainer::const_iterator it = sources.begin(); it != sources.end(); ++it)
  {
    PointSetType::PointType doublePoint;
    mitk::itk2vtk(*it, doublePoint); // copy mitk::ScalarType point into double point as workaround to ITK 3.10 bug
    sourcePointSet->SetPoint(i++, doublePoint /**it*/);
  }

  i = 0;
  PointSetType::Pointer targetPointSet = PointSetType::New();
  for (LandmarkPointContainer::const_iterator it = targets.begin(); it != targets.end(); ++it)
  {
    PointSetType::PointType doublePoint;
    mitk::itk2vtk(*it, doublePoint); // copy mitk::ScalarType point into double point as workaround to ITK 3.10 bug
    targetPointSet->SetPoint(i++, doublePoint /**it*/);
  }

  /* get centroid and extends of our pointsets */
  //BoundingBoxType::Pointer sourceBoundingBox = BoundingBoxType::New();
  //sourceBoundingBox->SetPoints(sourcePointSet->GetPoints());
  //sourceBoundingBox->ComputeBoundingBox();
  //BoundingBoxType::Pointer targetBoundingBox = BoundingBoxType::New();
  //targetBoundingBox->SetPoints(targetPointSet->GetPoints());
  //targetBoundingBox->ComputeBoundingBox();


  TransformType::Pointer transform = TransformType::New();
  transform->SetIdentity();
  //transform->SetTranslation(targetBoundingBox->GetCenter() - sourceBoundingBox->GetCenter());

  itk::LevenbergMarquardtOptimizer::Pointer optimizer = itk::LevenbergMarquardtOptimizer::New();
  optimizer->SetUseCostFunctionGradient(false);

  RegistrationType::Pointer registration = RegistrationType::New();

  // Scale the translation components of the Transform in the Optimizer
  itk::LevenbergMarquardtOptimizer::ScalesType scales(transform->GetNumberOfParameters());
  const double translationScale = 5000; //sqrtf(targetBoundingBox->GetDiagonalLength2())  * 1000; // dynamic range of translations
  const double rotationScale = 1.0; // dynamic range of rotations
  scales[0] = 1.0 / rotationScale;
  scales[1] = 1.0 / rotationScale;
  scales[2] = 1.0 / rotationScale;
  scales[3] = 1.0 / translationScale;
  scales[4] = 1.0 / translationScale;
  scales[5] = 1.0 / translationScale;
  //scales.Fill(0.01);
  unsigned long numberOfIterations = 80000;
  double gradientTolerance = 1e-10; // convergence criterion
  double valueTolerance = 1e-10; // convergence criterion
  double epsilonFunction = 1e-10; // convergence criterion
  optimizer->SetScales( scales );
  optimizer->SetNumberOfIterations( numberOfIterations );
  optimizer->SetValueTolerance( valueTolerance );
  optimizer->SetGradientTolerance( gradientTolerance );
  optimizer->SetEpsilonFunction( epsilonFunction );


  registration->SetInitialTransformParameters( transform->GetParameters() );
  //------------------------------------------------------
  // Connect all the components required for Registration
  //------------------------------------------------------
  MetricType::Pointer metric = MetricType::New();

  registration->SetMetric( metric );
  registration->SetOptimizer( optimizer );
  registration->SetTransform( transform );
  registration->SetFixedPointSet( targetPointSet );
  registration->SetMovingPointSet( sourcePointSet );

  try
  {
    //registration->StartRegistration();
    registration->Update();
  }
  catch( itk::ExceptionObject & e )
  {
    MITK_INFO << "Exception caught during ICP optimization: " << e;
    return false;
    //throw e;
//.........这里部分代码省略.........