C++ imagetype::Pointer类代码示例

FuzzyCMeans::FuzzyCMeans(QVector<ClusterPoint2D> &points, QVector<ClusterCentroid2D> &clusters, float fuzzy, imageType::Pointer myImage, int numCluster)
{
    this->Eps = std::pow(10, -5);
        this->isConverged=false;
        this->Points = points;
        this->Clusters = clusters;
        this->myImageHeight = myImage->GetBufferedRegion().GetSize()[0];
        this->myImageWidth = myImage->GetBufferedRegion().GetSize()[1];
        this->myImage = myImage;
        U= boost::numeric::ublas::matrix<double>(Points.size(),Clusters.size());
        this->Fuzzyness = fuzzy;

        double diff;


        imageType::SizeType size;
        size[0]=myImageWidth; // x axis
        size[1]=myImageHeight; // y
        imageType::RegionType region;
        region.SetSize(size);
        imageType::Pointer image = imageType::New();
        image->SetRegions(region);
        image->Allocate(); // immagine create

        // Iterate through all points to create initial U matrix
        for (int i = 0; i < Points.size()-1; i++)
        {
            ClusterPoint2D p = Points.at(i);
            double sum = 0.0;

            for (int j = 0; j < Clusters.size()-1; j++)
            {
                ClusterCentroid2D c = Clusters.at(j);
                diff = std::sqrt(std::pow(CalculateEuclideanDistance(p, c), 2.0));
                //U(i, j) = (diff == 0) ? Eps : diff;
                if (diff==0){
                    U(i,j)=Eps;
                }
                else{
                    U(i,j)=diff;
                }
                sum += U(i, j);
            }
         }

        this->RecalculateClusterMembershipValues();
}

// ------------------------------------------------------------------------
InitialTransformExtractor::TranslationType InitialTransformExtractor::ComputeTranslation(
    const ImageType::Pointer &image,
    const OptionsData &options)
{
    // set the offsets
    TranslationType translation;
    translation.Fill(0);
    translation[0] -= options.roiOffset[0];
    translation[1] -= options.roiOffset[1];
    translation[2] -= options.roiOffset[2];


    ImageType::PointType origin = image->GetOrigin();
    translation[0] -= origin[0];
    translation[1] -= origin[1];
    translation[2] -= origin[2];

    return translation;
}

void TestDeepCopyUnsignedCharScalar()
{
  itk::Index<2> corner = {{0,0}};
  itk::Size<2> size = {{10,10}};
  itk::ImageRegion<2> region(corner, size);

  typedef itk::Image<unsigned char, 2> ImageType;
  ImageType::Pointer imageIn = ImageType::New();
  imageIn->SetRegions(region);
  imageIn->Allocate();

  ImageType::Pointer imageOut = ImageType::New();
  ITKHelpers::DeepCopy(imageIn.GetPointer(), imageOut.GetPointer());
}

void TestBlurAllChannels()
{

  {
  typedef itk::VectorImage<float, 2> ImageType;
  ImageType::Pointer image = ImageType::New();
  ImageType::Pointer blurred = ImageType::New();

  float sigma = 2.0f;
  ITKHelpers::BlurAllChannels(image.GetPointer(), blurred.GetPointer(), sigma);
  }

  {
  typedef itk::VectorImage<float, 2> ImageType;
  ImageType::Pointer image = ImageType::New();
  ImageType::Pointer blurred = ImageType::New();

  float sigma = 2.0f;
  ITKHelpers::BlurAllChannels(image.GetPointer(), blurred.GetPointer(), sigma);
  }

}

// ------------------------------------------------------------------------
void OpenCVValve::ConvertImage(const ImageType::Pointer &input, MatPtr &mat)
{
	// cast the image to uchar 
	typedef itk::Image<unsigned char, 2> OutputImageType;
	typedef itk::RescaleIntensityImageFilter<ImageType, OutputImageType> CasterType;
	CasterType::Pointer caster = CasterType::New();
	caster->SetOutputMaximum(255);
	caster->SetOutputMinimum(0);
	caster->SetInput(input);
	caster->Update();


	OutputImageType::Pointer output = caster->GetOutput();


	typedef itk::ImageFileWriter<OutputImageType> WriterType;
	WriterType::Pointer writer = WriterType::New();
	writer->SetImageIO(itk::PNGImageIO::New());
	writer->SetInput(output);
	writer->SetFileName("test.png");
	writer->Update();



	ImageType::SizeType size = input->GetLargestPossibleRegion().GetSize();
	unsigned int rows = size[1];
	unsigned int cols = size[0];

	mat = new MatType(rows,cols, CV_8UC1);


	itk::ImageRegionConstIterator<OutputImageType> it(output, output->GetLargestPossibleRegion());
	it.GoToBegin();
	while(!it.IsAtEnd())
	{
		OutputImageType::IndexType index = it.GetIndex();
		unsigned char val = it.Get();
		mat->at<unsigned char>(cv::Point(index[0], index[1])) = val;
		++it;
	}
}

void TestUpsample()
{
  typedef itk::Image<unsigned char, 2> ImageType;
  ImageType::Pointer image = ImageType::New();

  itk::Index<2> corner = {{0,0}};
  itk::Size<2> size = {{100,100}};
  itk::ImageRegion<2> region(corner, size);

  image->SetRegions(region);
  image->Allocate();

  ImageType::Pointer upsampled = ImageType::New();

  ITKHelpers::Upsample(image.GetPointer(), 2, upsampled.GetPointer());

  std::cout << "TestUpsample() output size: "
            << upsampled->GetLargestPossibleRegion().GetSize() << std::endl;
}

bool VolumeProcess:: RunFillingZeroOnBouandary(int Bx,int By,int Bz)
{
	ImageType::Pointer tempImg = ImageType::New();
	tempImg->SetRegions( m_outputImage->GetLargestPossibleRegion() );
	tempImg->Allocate();
	tempImg->FillBuffer(0);

	ImageType::RegionType fullRegion = m_outputImage->GetBufferedRegion();

	int numColumns = fullRegion.GetSize(0);
	int numRows = fullRegion.GetSize(1);
	int numStacks = fullRegion.GetSize(2);
	int i, j,k; 

	itk::ImageRegionIteratorWithIndex< ImageType > itr( m_outputImage, fullRegion );

	for(itr.GoToBegin(); !itr.IsAtEnd(); ++itr)
	{
		
		ImageType::IndexType index = itr.GetIndex();
	    ImageType::PixelType pix = m_outputImage->GetPixel(index);
        i = index[0];
		j = index[1];
		k = index[2];
		if (i < Bx || i > numColumns -Bx || j < By || j > numRows-By ||k <Bz ||k > numStacks-Bz )
		    tempImg->SetPixel(itr.GetIndex(), 0);
		else 
			tempImg->SetPixel(itr.GetIndex(), pix);
     }
		//Copy temp img back to image 
	itk::ImageRegionIterator< ImageType > itr1( tempImg, tempImg->GetLargestPossibleRegion());
	itk::ImageRegionIterator< ImageType > itr2( m_outputImage, m_outputImage->GetLargestPossibleRegion());
	for(itr1.GoToBegin(), itr2.GoToBegin() ; !itr1.IsAtEnd(); ++itr1, ++itr2)
	{
	  itr2.Set( itr1.Get() );
	}

	if(debug)
		std::cerr << "RunFillingZero Done" << std::endl;

	return true;
  
}

  void setUp(void)
  {
    // Load Image Data
    m_Image = dynamic_cast<mitk::Image*>(mitk::IOUtil::Load(GetTestDataFilePath("Pic3D.nrrd"))[0].GetPointer());
    mitk::CastToItkImage(m_Image,m_ItkImage);

    // Create a single mask with only one pixel within the regions
    mitk::Image::Pointer mask1 = dynamic_cast<mitk::Image*>(mitk::IOUtil::Load(GetTestDataFilePath("Pic3D.nrrd"))[0].GetPointer());
    mitk::CastToItkImage(mask1,m_ItkMask);
    m_ItkMask->FillBuffer(0);
    MaskType::IndexType index;
    index[0]=88;index[1]=81;index[2]=13;
    m_ItkMask->SetPixel(index, 1);
    MITK_INFO << "Pixel Value: "<<m_ItkImage->GetPixel(index);
    mitk::CastToMitkImage(m_ItkMask, m_Mask);

    // Create a mask with a covered region
    mitk::Image::Pointer lmask1 = dynamic_cast<mitk::Image*>(mitk::IOUtil::Load(GetTestDataFilePath("Pic3D.nrrd"))[0].GetPointer());
    mitk::CastToItkImage(lmask1,m_ItkMask1);
    m_ItkMask1->FillBuffer(0);
    int range=2;
    for (int x = 88-range;x < 88+range+1;++x)
    {
      for (int y=81-range;y<81+range+1;++y)
      {
        for (int z=13-range;z<13+range+1;++z)
        {
          index[0] = x;
          index[1] = y;
          index[2] = z;
          //MITK_INFO << "Pixel: " <<m_ItkImage->GetPixel(index);
          m_ItkMask1->SetPixel(index, 1);
        }
      }
    }
    mitk::CastToMitkImage(m_ItkMask1, m_Mask1);

    m_GradientImage=GenerateGradientWithDimXImage<unsigned char>(5,5,5);
    m_GradientMask = GenerateMaskImage<unsigned char>(5,5,5);
  }

bool TestDeepCopyFloatScalar()
{
  itk::Index<2> corner = {{0,0}};
  itk::Size<2> size = {{10,10}};
  itk::ImageRegion<2> region(corner, size);

  typedef itk::Image<float, 2> ImageType;
  ImageType::Pointer imageIn = ImageType::New();
  imageIn->SetRegions(region);
  imageIn->Allocate();

  ImageType::Pointer imageOut = ImageType::New();
  ITKHelpers::DeepCopy(imageIn.GetPointer(), imageOut.GetPointer());

  return true;
}

bool TestRandomImage()
{
  typedef itk::Image<unsigned char, 2> ImageType;
  ImageType::Pointer image = ImageType::New();

  itk::Index<2> corner = {{0,0}};
  itk::Size<2> size = {{100,100}};
  itk::ImageRegion<2> region(corner, size);

  image->SetRegions(region);
  image->Allocate();
  image->FillBuffer(0);

  ITKHelpers::RandomImage(image.GetPointer());
  ITKHelpers::WriteImage(image.GetPointer(), "random.png");

  return true;
}

void TestDrawRectangle()
{
  typedef itk::Image<unsigned char, 2> ImageType;
  ImageType::Pointer image = ImageType::New();

  itk::Index<2> corner = {{0,0}};
  itk::Size<2> size = {{100,100}};
  itk::ImageRegion<2> region(corner, size);

  image->SetRegions(region);
  image->Allocate();
  image->FillBuffer(0);

  itk::Index<2> corner0 = {{10,10}};
  itk::Index<2> corner1 = {{30,30}};
  ITKHelpers::DrawRectangle(image.GetPointer(), 255,
                            corner0, corner1);

  ITKHelpers::WriteImage(image.GetPointer(), "rectangle.png");
}

bool TestComputeGradientsInRegion()
{
  itk::Index<2> imageCorner = {{0,0}};
  itk::Size<2> imageSize = {{100,100}};
  itk::ImageRegion<2> imageRegion(imageCorner, imageSize);

  typedef itk::Image<unsigned char, 2> ImageType;
  ImageType::Pointer image = ImageType::New();
  image->SetRegions(imageRegion);
  image->Allocate();

  itk::ImageRegionIterator<ImageType> imageIterator(image, image->GetLargestPossibleRegion());

  while(!imageIterator.IsAtEnd())
  {
    imageIterator.Set(rand() % 255);

    ++imageIterator;
  }

  ITKHelpers::WriteImage(image.GetPointer(), "Image.mha");

  itk::Index<2> regionCorner = {{50,50}};
  itk::Size<2> regionSize = {{10,10}};
  itk::ImageRegion<2> region(regionCorner, regionSize);

  typedef itk::Image<itk::CovariantVector<float, 2>, 2> GradientImageType;
  GradientImageType::Pointer gradientImage = GradientImageType::New();

  std::cout << "Computing gradients..." << std::endl;
  ITKHelpers::ComputeGradientsInRegion(image.GetPointer(), region, gradientImage.GetPointer());

  std::cout << "Writing..." << std::endl;
  ITKHelpers::WriteImage(gradientImage.GetPointer(), "GradientImage.mha");

  return true;
}

// Run with: Data/trashcan.mha Data/trashcan_mask.mha 15 Data/trashcan.vtp Intensity filled.mha
int main(int argc, char *argv[])
{
  // Verify arguments
  if(argc != 6)
    {
    std::cerr << "Required arguments: image.mha imageMask.mha patch_half_width normals.vts output.mha" << std::endl;
    std::cerr << "Input arguments: ";
    for(int i = 1; i < argc; ++i)
      {
      std::cerr << argv[i] << " ";
      }
    return EXIT_FAILURE;
    }

  // Parse arguments
  std::string imageFilename = argv[1];
  std::string maskFilename = argv[2];

  std::stringstream ssPatchRadius;
  ssPatchRadius << argv[3];
  unsigned int patch_half_width = 0;
  ssPatchRadius >> patch_half_width;

  std::string normalsFileName = argv[4];

  std::string outputFilename = argv[5];

  // Output arguments
  std::cout << "Reading image: " << imageFilename << std::endl;
  std::cout << "Reading mask: " << maskFilename << std::endl;
  std::cout << "Patch half width: " << patch_half_width << std::endl;
  std::cout << "Reading normals: " << normalsFileName << std::endl;
  std::cout << "Output: " << outputFilename << std::endl;

  vtkSmartPointer<vtkXMLStructuredGridReader> structuredGridReader = vtkSmartPointer<vtkXMLStructuredGridReader>::New();
  structuredGridReader->SetFileName(normalsFileName.c_str());
  structuredGridReader->Update();
  
  typedef FloatVectorImageType ImageType;

  typedef  itk::ImageFileReader<ImageType> ImageReaderType;
  ImageReaderType::Pointer imageReader = ImageReaderType::New();
  imageReader->SetFileName(imageFilename);
  imageReader->Update();

  ImageType::Pointer image = ImageType::New();
  ITKHelpers::DeepCopy(imageReader->GetOutput(), image.GetPointer());

  Mask::Pointer mask = Mask::New();
  mask->Read(maskFilename);

  std::cout << "hole pixels: " << mask->CountHolePixels() << std::endl;
  std::cout << "valid pixels: " << mask->CountValidPixels() << std::endl;

  typedef ImagePatchPixelDescriptor<ImageType> ImagePatchPixelDescriptorType;
  typedef FeatureVectorPixelDescriptor FeatureVectorPixelDescriptorType;

  // Create the graph
  typedef boost::grid_graph<2> VertexListGraphType;
  boost::array<std::size_t, 2> graphSideLengths = { { imageReader->GetOutput()->GetLargestPossibleRegion().GetSize()[0],
                                                      imageReader->GetOutput()->GetLargestPossibleRegion().GetSize()[1] } };
  VertexListGraphType graph(graphSideLengths);
  typedef boost::graph_traits<VertexListGraphType>::vertex_descriptor VertexDescriptorType;

  // Get the index map
  typedef boost::property_map<VertexListGraphType, boost::vertex_index_t>::const_type IndexMapType;
  IndexMapType indexMap(get(boost::vertex_index, graph));

  // Create the priority map
  typedef boost::vector_property_map<float, IndexMapType> PriorityMapType;
  PriorityMapType priorityMap(num_vertices(graph), indexMap);

  // Create the node fill status map. Each pixel is either filled (true) or not filled (false).
  typedef boost::vector_property_map<bool, IndexMapType> FillStatusMapType;
  FillStatusMapType fillStatusMap(num_vertices(graph), indexMap);

  // Create the boundary status map. A node is on the current boundary if this property is true. 
  // This property helps the boundaryNodeQueue because we can mark here if a node has become no longer
  // part of the boundary, so when the queue is popped we can check this property to see if it should
  // actually be processed.
  typedef boost::vector_property_map<bool, IndexMapType> BoundaryStatusMapType;
  BoundaryStatusMapType boundaryStatusMap(num_vertices(graph), indexMap);

  // Create the descriptor map. This is where the data for each pixel is stored.
  typedef boost::vector_property_map<ImagePatchPixelDescriptorType, IndexMapType> ImagePatchDescriptorMapType;
  ImagePatchDescriptorMapType imagePatchDescriptorMap(num_vertices(graph), indexMap);

  // Create the descriptor map. This is where the data for each pixel is stored.
  typedef boost::vector_property_map<FeatureVectorPixelDescriptorType, IndexMapType> FeatureVectorDescriptorMapType;
  FeatureVectorDescriptorMapType featureVectorDescriptorMap(num_vertices(graph), indexMap);

  // Create the patch inpainter. The inpainter needs to know the status of each pixel to determine if they should be inpainted.
  typedef MaskedGridPatchInpainter<FillStatusMapType> InpainterType;
  InpainterType patchInpainter(patch_half_width, fillStatusMap);

  // Create the priority function
  typedef PriorityRandom PriorityType;
  PriorityType priorityFunction;

//.........这里部分代码省略.........

int
main(int argc, char* argv[])
{
	boost::filesystem::path inputFile;
	boost::filesystem::path outputFile;

	Algorithm algorithm;
	std::string algorithmName;

	po::options_description desc("Allowed options");
	desc.add_options()
		("help", "produce help message")
		("input,i", po::value<boost::filesystem::path>(&inputFile), "input file")
		("algorithm,a", po::value<std::string>(&algorithmName)->default_value("itk-implementation"), "itk-implementation")
		("output,o", po::value<boost::filesystem::path>(&outputFile), "output mask file")
		;

	po::variables_map vm;
	po::store(po::parse_command_line(argc, argv, desc), vm);
	po::notify(vm);

	algorithm = getAlgorithmFromString(algorithmName);

	if (vm.count("help")) {
		std::cout << desc << "\n";
		return 1;
	}

	if (vm.count("input") == 0) {
		std::cout << "Missing input filename\n" << desc << "\n";
		return 1;
	}

	if (vm.count("output") == 0) {
		std::cout << "Missing output filename\n" << desc << "\n";
		return 1;
	}

	typedef itk::ImageFileReader<ImageType>  ImageReaderType;

	BOOST_LOG_TRIVIAL(info) << "Loading inputs ...";
	ImageReaderType::Pointer imageReader = ImageReaderType::New();
	imageReader->SetFileName(inputFile.string());
	imageReader->Update();

	ImageType::Pointer image = imageReader->GetOutput();

	LabeledImageType::Pointer outputImage;

	switch (algorithm) {
	case Algorithm::ItkImplementation: {
		BOOST_LOG_TRIVIAL(info) << "Running ITK version of watershed transformation ...";
			WatershedFilterType::Pointer filter = WatershedFilterType::New();
			filter->SetInput(image);
			filter->MarkWatershedLineOff();

			boost::timer::cpu_timer computationTimer;
			computationTimer.start();
			filter->Update();
			BOOST_LOG_TRIVIAL(info) << "Computation time: " << computationTimer.format(9, "%w");

			outputImage = filter->GetOutput();
		}
		break;
	default:
		BOOST_LOG_TRIVIAL(info) << "Allocating output ...";
		outputImage = LabeledImageType::New();
		outputImage->SetRegions(image->GetLargestPossibleRegion());
		outputImage->Allocate();
		outputImage->SetSpacing(image->GetSpacing());
		BOOST_LOG_TRIVIAL(info) << "Running CUDA version of watershed transformation (topographical distance)...";
		boost::timer::cpu_timer computationTimer;
		computationTimer.start();
		watershedTransformation2(
			cugip::const_view(*(image.GetPointer())),
			cugip::view(*(outputImage.GetPointer())),
			Options()
			);
		BOOST_LOG_TRIVIAL(info) << "Computation time: " << computationTimer.format(9, "%w");
		//BOOST_LOG_TRIVIAL(error) << "Unknown algorithm";
	}

	BOOST_LOG_TRIVIAL(info) << "Saving output `" << outputFile << "` ...";
	WriterType::Pointer writer = WriterType::New();
	writer->SetFileName(outputFile.string());
	writer->SetInput(outputImage);
	try {
		writer->Update();
	} catch (itk::ExceptionObject & error) {
		std::cerr << "Error: " << error << std::endl;
		return EXIT_FAILURE;
	}

	return 0;
}

void LabelsVolumeGenerator::postProcessTargetITK()
{
    typedef unsigned char       PixelType;
    typedef Image<PixelType, 3> ImageType;

    // create a new image from the target data

    ImageType::Pointer input = ImageType::New();

    ImageType::IndexType start;
    start.Fill(0);

    ImageType::SizeType size;
    for (int i = 0; i < 3; ++i) size[i] = m_targetVolume->dimensions[i];

    ImageType::RegionType region;
    region.SetSize(size);
    region.SetIndex(start);

    input->SetRegions(region);
    input->SetSpacing(m_targetVolume->spacing.data());
    input->Allocate();

    memcpy(input->GetBufferPointer(), m_targetVolume->data, m_bufferSize);

    // create a grayscale dilation filter and a structuring element

    typedef BinaryBallStructuringElement<PixelType, 3> StructureType;
    typedef GrayscaleDilateImageFilter<ImageType, ImageType, StructureType> DilateFilterType;

    DilateFilterType::Pointer filter = DilateFilterType::New();

    StructureType structure;
    structure.SetRadius(1);

    filter->SetKernel(structure);

    // set up progress reporting
    if (m_progressReporter)
    {
        CStyleCommand::Pointer command = CStyleCommand::New();
        command->SetClientData(m_progressReporter);
        command->SetCallback(ProgressCallback);
        filter->AddObserver(ProgressEvent(), command);
        m_progressReporter->start("Post-Processing Label volume",
                                  "Dilating label field...");
    }

    // hook up the filters and run

    filter->SetInput(input);
    filter->Update();

    if (m_progressReporter)
        m_progressReporter->finish();

    // copy back into the target volume data
    memcpy(m_targetVolume->data, filter->GetOutput()->GetBufferPointer(), m_bufferSize);

    // threshold and gaussian blur to put a smooth version into the alpha channel

    typedef BinaryThresholdImageFilter<ImageType, ImageType> ThresholderType;
//    typedef DiscreteGaussianImageFilter<ImageType, ImageType> SmoothFilterType;
    typedef SmoothingRecursiveGaussianImageFilter<ImageType, ImageType> SmoothFilterType;

    ThresholderType::Pointer thresholder = ThresholderType::New();
    thresholder->SetLowerThreshold(1);
    thresholder->SetUpperThreshold(255);
    thresholder->SetInsideValue(255);
    thresholder->SetOutsideValue(0);

    SmoothFilterType::Pointer smoother = SmoothFilterType::New();
//    smoother->SetVariance(0.05); // in physical units
//    smoother->SetMaximumKernelWidth(5);
    smoother->SetSigma(0.2);

    // set up progress reporting (again)
    if (m_progressReporter)
    {
        CStyleCommand::Pointer command = CStyleCommand::New();
        command->SetClientData(m_progressReporter);
        command->SetCallback(ProgressCallback);
        smoother->AddObserver(ProgressEvent(), command);
        m_progressReporter->start("Post-Processing Label volume",
                                  "Smoothing alpha mask...");
    }

    // hook up the filters and run

    thresholder->SetInput(input);
    smoother->SetInput(thresholder->GetOutput());
    smoother->Update();

    // copy back into the target volume data
    memcpy(m_targetMask->data, smoother->GetOutput()->GetBufferPointer(), m_bufferSize);

    if (m_progressReporter)
        m_progressReporter->finish();
}