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();
}

int main( int argc, char* argv[] )
{
	std::string input_file;
	std::string output_file;
	float variance;
	int patchRadius;
	int searchRadius;
	po::options_description desc("Allowed options");
	desc.add_options()
		("help", "produce help message")
		("input,i", po::value<std::string>(&input_file), "input file")
		("output,o", po::value<std::string>(&output_file), "output file")
		("variance,v", po::value<float>(&variance)->default_value(10.0), "variance")
		("patch-radius,p", po::value<int>(&patchRadius)->default_value(2), "patch radius")
		("search-radius,s", po::value<int>(&searchRadius)->default_value(4), "search radius")
		;

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

	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;
	}


	const unsigned int Dimension = 3;

	//typedef float                              PixelType;
	//typedef itk::Image< PixelType, Dimension > ImageType;

	typedef itk::ImageFileReader<ImageType>  ReaderType;
	typedef itk::ImageFileWriter<ImageType> WriterType;

	ReaderType::Pointer reader = ReaderType::New();
	reader->SetFileName(input_file);
	reader->Update();

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

	ImageType::Pointer output_image = ImageType::New();
	output_image->SetRegions(image->GetLargestPossibleRegion());
	output_image->Allocate();
	output_image->SetSpacing(image->GetSpacing());

	//denoise(image, output_image);
	denoise(
		image->GetPixelContainer()->GetBufferPointer(),
		output_image->GetPixelContainer()->GetBufferPointer(),
		image->GetLargestPossibleRegion().GetSize()[0],
		image->GetLargestPossibleRegion().GetSize()[1],
		image->GetLargestPossibleRegion().GetSize()[2],
		variance,
		patchRadius,
		searchRadius);


	WriterType::Pointer writer = WriterType::New();
	writer->SetFileName(output_file);
	writer->SetInput(output_image);

	try
	{
		writer->Update();
	}
	catch( itk::ExceptionObject & error )
	{
		std::cerr << "Error: " << error << std::endl;
		return EXIT_FAILURE;
	}

	return EXIT_SUCCESS;
}

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

	typedef itk::Image<float,3> ImageType;

	typedef itk::ImageFileReader<ImageType> ReaderType;

	typedef itk::ImageFileWriter<ImageType> WriterType;

	ReaderType::Pointer reader = ReaderType::New();

	reader->SetFileName(argv[1]);

	typedef itk::HessianRecursiveGaussianImageFilter<ImageType> HessianFilterType;
	typedef HessianFilterType::OutputImageType HessianImageType;
	typedef HessianImageType::PixelType HessianType;

	HessianFilterType::Pointer hessianFilter = HessianFilterType::New();

	hessianFilter->SetSigma(0.05);
	hessianFilter->SetInput(reader->GetOutput());

	hessianFilter->Update();

	HessianImageType::Pointer hessianImage = hessianFilter->GetOutput();
	hessianImage->DisconnectPipeline();

	typedef itk::ImageRegionConstIteratorWithIndex<HessianImageType> IteratorType;

	IteratorType hessianIterator(hessianImage,hessianImage->GetLargestPossibleRegion());

	typedef HessianType::EigenValuesArrayType EigenvalueArrayType;
	EigenvalueArrayType expectation;
	expectation.Fill(0.0);
	unsigned count=0;

	ImageType::Pointer small = ImageType::New();
	small->CopyInformation(reader->GetOutput());

	small->SetRegions(reader->GetOutput()->GetLargestPossibleRegion());
	small->Allocate();
	small->FillBuffer(0.0);
	ImageType::Pointer medium = ImageType::New();
	medium->CopyInformation(reader->GetOutput());
	medium->SetRegions(reader->GetOutput()->GetLargestPossibleRegion());

	medium->Allocate();
	medium->FillBuffer(0.0);
	ImageType::Pointer large = ImageType::New();
	large->CopyInformation(reader->GetOutput());
	large->SetRegions(reader->GetOutput()->GetLargestPossibleRegion());
	large->Allocate();
	large->FillBuffer(0.0);

	while(!hessianIterator.IsAtEnd()){

		auto hessian = hessianIterator.Get();
		EigenvalueArrayType eigenValues;
		hessian.ComputeEigenValues(eigenValues);

		std::sort(eigenValues.Begin(),eigenValues.End(),[](double & a, double &b){
			return std::abs(a) > std::abs(b);
		});

		//std::cout << eigenValues[0] << "\t"<< eigenValues[1] << "\t" << eigenValues[2] <<std::endl;

		if(eigenValues[0]<0){
			expectation[0]+=eigenValues[0];
			expectation[1]+=eigenValues[1];
			expectation[2]+=eigenValues[2];

			count++;

			large->SetPixel(hessianIterator.GetIndex(),std::abs(eigenValues[0]));
			medium->SetPixel(hessianIterator.GetIndex(),std::abs(eigenValues[1]));
			small->SetPixel(hessianIterator.GetIndex(),std::abs(eigenValues[2]));
		}
		++hessianIterator;
	}
	expectation[0]=expectation[0]/count;
	expectation[1]=expectation[1]/count;
	expectation[2]=expectation[2]/count;

	WriterType::Pointer smallWriter = WriterType::New();
	smallWriter->SetInput(small);
	smallWriter->SetFileName("small.mha");
	smallWriter->Update();

	WriterType::Pointer mediumWriter = WriterType::New();
	mediumWriter->SetInput(medium);
	mediumWriter->SetFileName("medium.mha");
	mediumWriter->Update();

	WriterType::Pointer largeWriter = WriterType::New();
	largeWriter->SetInput(large);
	largeWriter->SetFileName("large.mha");
	largeWriter->Update();

	std::cout << "Expectation: \t" << expectation<< std::endl;

}

void QmitkBasicImageProcessing::StartButtonClicked()
{
  if(!m_SelectedImageNode->GetNode()) return;

  this->BusyCursorOn();

  mitk::Image::Pointer newImage;

  try
  {
    newImage = dynamic_cast<mitk::Image*>(m_SelectedImageNode->GetNode()->GetData());
  }
  catch ( std::exception &e )
  {
  QString exceptionString = "An error occured during image loading:\n";
  exceptionString.append( e.what() );
    QMessageBox::warning( NULL, "Basic Image Processing", exceptionString , QMessageBox::Ok, QMessageBox::NoButton );
    this->BusyCursorOff();
    return;
  }

  // check if input image is valid, casting does not throw exception when casting from 'NULL-Object'
  if ( (! newImage) || (newImage->IsInitialized() == false) )
  {
    this->BusyCursorOff();

    QMessageBox::warning( NULL, "Basic Image Processing", "Input image is broken or not initialized. Returning.", QMessageBox::Ok, QMessageBox::NoButton );
    return;
  }

  // check if operation is done on 4D a image time step
  if(newImage->GetDimension() > 3)
  {
    mitk::ImageTimeSelector::Pointer timeSelector = mitk::ImageTimeSelector::New();
    timeSelector->SetInput(newImage);
    timeSelector->SetTimeNr( ((QmitkSliderNavigatorWidget*)m_Controls->sliceNavigatorTime)->GetPos() );
    timeSelector->Update();
    newImage = timeSelector->GetOutput();
  }



  // check if image or vector image
  ImageType::Pointer itkImage = ImageType::New();
  VectorImageType::Pointer itkVecImage = VectorImageType::New();

  int isVectorImage = newImage->GetPixelType().GetNumberOfComponents();

  if(isVectorImage > 1)
  {
    CastToItkImage( newImage, itkVecImage );
  }
  else
  {
    CastToItkImage( newImage, itkImage );
  }

  std::stringstream nameAddition("");

  int param1 = m_Controls->sbParam1->value();
  int param2 = m_Controls->sbParam2->value();
  double dparam1 = m_Controls->dsbParam1->value();
  double dparam2 = m_Controls->dsbParam2->value();
  double dparam3 = m_Controls->dsbParam3->value();

  try{

  switch (m_SelectedAction)
  {

  case GAUSSIAN:
    {
      GaussianFilterType::Pointer gaussianFilter = GaussianFilterType::New();
      gaussianFilter->SetInput( itkImage );
      gaussianFilter->SetVariance( param1 );
      gaussianFilter->UpdateLargestPossibleRegion();
      newImage = mitk::ImportItkImage(gaussianFilter->GetOutput())->Clone();
      nameAddition << "_Gaussian_var_" << param1;
      std::cout << "Gaussian filtering successful." << std::endl;
      break;
    }

  case MEDIAN:
    {
      MedianFilterType::Pointer medianFilter = MedianFilterType::New();
      MedianFilterType::InputSizeType size;
      size.Fill(param1);
      medianFilter->SetRadius( size );
      medianFilter->SetInput(itkImage);
      medianFilter->UpdateLargestPossibleRegion();
      newImage = mitk::ImportItkImage(medianFilter->GetOutput())->Clone();
      nameAddition << "_Median_radius_" << param1;
      std::cout << "Median Filtering successful." << std::endl;
      break;
    }

  case TOTALVARIATION:
    {
      if(isVectorImage > 1)
      {
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
	/* % MINIMALPATH Recherche du chemin minimum de Haut vers le bas et de
	% bas vers le haut tel que dÈcrit par Luc Vincent 1998
	% [sR,sC,S] = MinimalPath(I,factx)
	%
	%   I     : Image d'entrÔøΩe dans laquelle on doit trouver le
	%           chemin minimal
	%   factx : Poids de linearite [1 10]
	%
	% Programme par : Ramnada Chav
	% ModifiÈ le 16 novembre 2007 */

	    typedef itk::ImageDuplicator< ImageType > DuplicatorType3D;
	    typedef itk::InvertIntensityImageFilter <ImageType> InvertIntensityImageFilterType;
	    typedef itk::StatisticsImageFilter<ImageType> StatisticsImageFilterType;

	    ImageType::Pointer inverted_image = image;

//	    if (invert)
//	    {
//	        StatisticsImageFilterType::Pointer statisticsImageFilterInput = StatisticsImageFilterType::New();
//	        statisticsImageFilterInput->SetInput(image);
//	        statisticsImageFilterInput->Update();
//	        double maxIm = statisticsImageFilterInput->GetMaximum();
//	        InvertIntensityImageFilterType::Pointer invertIntensityFilter = InvertIntensityImageFilterType::New();
//	        invertIntensityFilter->SetInput(image);
//	        invertIntensityFilter->SetMaximum(maxIm);
//	        invertIntensityFilter->Update();
//	        inverted_image = invertIntensityFilter->GetOutput();
//	    }


	    ImageType::SizeType sizeImage = image->GetLargestPossibleRegion().GetSize();
	    int m = sizeImage[0]; // x to change because we are in AIL
	    int n = sizeImage[2]; // y
	    int p = sizeImage[1]; // z

	    // create image with high values J1
	    DuplicatorType3D::Pointer duplicator = DuplicatorType3D::New();
	    duplicator->SetInputImage(inverted_image);
	    duplicator->Update();
	    ImageType::Pointer J1 = duplicator->GetOutput();
	    typedef itk::ImageRegionIterator<ImageType> ImageIterator3D;
	    ImageIterator3D vItJ1( J1, J1->GetBufferedRegion() );
	    vItJ1.GoToBegin();
	    while ( !vItJ1.IsAtEnd() )
	    {
	        vItJ1.Set(100000000);
	        ++vItJ1;
	    }

	    // create image with high values J2
	    DuplicatorType3D::Pointer duplicatorJ2 = DuplicatorType3D::New();
	    duplicatorJ2->SetInputImage(inverted_image);
	    duplicatorJ2->Update();
	    ImageType::Pointer J2 = duplicatorJ2->GetOutput();
	    ImageIterator3D vItJ2( J2, J2->GetBufferedRegion() );
	    vItJ2.GoToBegin();
	    while ( !vItJ2.IsAtEnd() )
	    {
	        vItJ2.Set(100000000);
	        ++vItJ2;
	    }

	    DuplicatorType3D::Pointer duplicatorCPixel = DuplicatorType3D::New();
	    duplicatorCPixel->SetInputImage(inverted_image);
	    duplicatorCPixel->Update();
	    ImageType::Pointer cPixel = duplicatorCPixel->GetOutput();

	    ImageType::IndexType index;

	    // iterate on slice from slice 1 (start=0) to slice p-2. Basically, we avoid first and last slices.
	    // IMPORTANT: first slice of J1 and last slice of J2 must be set to 0...
	    for (int x=0; x<m; x++)
	    {
	        for (int y=0; y<n; y++)
	        {
	            index[0] = x; index[1] = 0; index[2] = y;
	            J1->SetPixel(index, 0.0);
	        }
	    }
	    for (int slice=1; slice<p; slice++)
	    {
	        // 1. extract pJ = the (slice-1)th slice of the image J1
	        Matrice pJ = Matrice(m,n);
	        for (int x=0; x<m; x++)
	        {
	            for (int y=0; y<n; y++)
	            {
	                index[0] = x; index[1] = slice-1; index[2] = y;
	                pJ(x,y) = J1->GetPixel(index);
	            }
	        }

	        // 2. extract cP = the (slice)th slice of the image cPixel
	        Matrice cP = Matrice(m,n);
	        for (int x=0; x<m; x++)
	        {
	            for (int y=0; y<n; y++)
//.........这里部分代码省略.........

void NrrdTensorImageReader
::GenerateData()
{
    if ( m_FileName == "")
    {
        throw itk::ImageFileReaderException(__FILE__, __LINE__, "Sorry, the filename is empty!");
    }
    else
    {
        try
        {
            const std::string& locale = "C";
            const std::string& currLocale = setlocale( LC_ALL, NULL );