int main()
{
  FloatVectorImageType::Pointer image = FloatVectorImageType::New();
  Testing::GetBlankImage(image.GetPointer(), 4);

  Mask::Pointer mask = Mask::New();
  Testing::GetFullyValidMask(mask.GetPointer());

  UnsignedCharScalarImageType::Pointer manualPriorityImage = UnsignedCharScalarImageType::New();
  Testing::GetBlankImage(manualPriorityImage.GetPointer());

  unsigned int patchRadius = 5;
  typedef PriorityConfidence ConfidencePriorityType;
  ConfidencePriorityType priorityConfidence(mask, patchRadius);
  
  PriorityManual<itk::Index<2>, UnsignedCharScalarImageType, ConfidencePriorityType>
      priority(manualPriorityImage, &priorityConfidence);

  itk::Index<2> filledPixel = {{0,0}};
  priority.Update(filledPixel);
  priority.SetManualPriorityImage(manualPriorityImage);

  itk::Index<2> queryPixel = {{0,0}};
  priority.ComputePriority(queryPixel);

  return EXIT_SUCCESS;
}

void WriteImagePixelsToRGBSpace(const FloatVectorImageType::Pointer image, const std::string& outputFileName)
{
  std::cout << "WriteImagePixelsToRGBSpace()" << std::endl;
  vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
  vtkSmartPointer<vtkUnsignedCharArray> colors = vtkSmartPointer<vtkUnsignedCharArray>::New();
  colors->SetName("Colors");
  colors->SetNumberOfComponents(3);

  itk::ImageRegionConstIterator<FloatVectorImageType> imageIterator(image, image->GetLargestPossibleRegion());
  while(!imageIterator.IsAtEnd())
    {
    FloatVectorImageType::PixelType pixel = imageIterator.Get();
    points->InsertNextPoint(pixel[0], pixel[1], pixel[2]);
    // TODO: Narrowing conversion warning
    //unsigned char color[3] = {pixel[0], pixel[1], pixel[2]};
    unsigned char color[3]; // TODO: Narrowing conversion warning
    colors->InsertNextTupleValue(color);
    ++imageIterator;
    }

  vtkSmartPointer<vtkPolyData> polyData = vtkSmartPointer<vtkPolyData>::New();
  polyData->SetPoints(points);
  polyData->GetPointData()->SetScalars(colors);

  vtkSmartPointer<vtkVertexGlyphFilter> glyphFilter = vtkSmartPointer<vtkVertexGlyphFilter>::New();
  glyphFilter->SetInputData(polyData);
  glyphFilter->Update();

  vtkSmartPointer<vtkXMLPolyDataWriter> writer = vtkSmartPointer<vtkXMLPolyDataWriter>::New();
  writer->SetInputConnection(glyphFilter->GetOutputPort());
  writer->SetFileName(outputFileName.c_str());
  writer->Write();

}

int main()
{
  FloatVectorImageType::Pointer image = FloatVectorImageType::New();
  Testing::GetBlankImage(image.GetPointer(), 4);

  Mask::Pointer mask = Mask::New();
  Testing::GetFullyValidMask(mask.GetPointer());

  unsigned int patchRadius = 5;
  PriorityDepth<FloatVectorImageType> priority(image, mask, patchRadius);

  itk::Index<2> filledPixel = {{0,0}};
  priority.Update(filledPixel);

  itk::Index<2> queryPixel = {{0,0}};
  priority.ComputePriority(queryPixel);

  return EXIT_SUCCESS;
}

// Convert a vector ITK image to a VTK image for display
void ITKImagetoVTKImage(FloatVectorImageType::Pointer image, vtkImageData* outputImage)
{
  std::cout << "ITKImagetoVTKImage()" << std::endl;
  if(image->GetNumberOfComponentsPerPixel() >= 3)
    {
    ITKImagetoVTKRGBImage(image, outputImage);
    }
  else
    {
    ITKImagetoVTKMagnitudeImage(image, outputImage);
    }
}

int main(int argc, char *argv[])
{
  if(argc < 3)
  {
    throw std::runtime_error("Required arguments: inputFileName outputPrefix");
  }

  std::string inputFileName = argv[1];
  std::string outputPrefix = argv[2];

  std::cout << "Input: " << inputFileName << std::endl;
  std::cout << "Output prefix: " << outputPrefix << std::endl;

  typedef itk::ImageFileReader<FloatVectorImageType> ReaderType;
  ReaderType::Pointer reader = ReaderType::New();
  reader->SetFileName(inputFileName);
  reader->Update();

  ITKHelpers::WriteVectorImageAsRGB(reader->GetOutput(), outputPrefix + "/Image.mha");

  WriteImagePixelsToRGBSpace(reader->GetOutput(), outputPrefix + "/ImageColors.vtp");

  WriteClusteredPixelsInRGBSpace(reader->GetOutput(), 20, outputPrefix + "/ImageColorsClustered.vtp");

  FloatVectorImageType::Pointer blurred = FloatVectorImageType::New();
  //float blurVariance = 2.0f; // almost no visible blurring
  //float blurVariance = 10.0f; // slight blurring of concrete
  //float blurVariance = 30.0f;

  // TODO: Update this call to the new API
  //Helpers::AnisotropicBlurAllChannels<FloatVectorImageType>(reader->GetOutput(), blurred, blurVariance);

  ITKHelpers::WriteVectorImageAsRGB(blurred.GetPointer(), outputPrefix + "/BlurredImage.mha");

  WriteImagePixelsToRGBSpace(blurred, outputPrefix + "/BlurredImageColors.vtp");

  WriteClusteredPixelsInRGBSpace(blurred, 20, outputPrefix + "/BlurredImageColorsClustered.vtp");

  return EXIT_SUCCESS;
}

void WriteClusteredPixelsInRGBSpace(const FloatVectorImageType::Pointer image, const unsigned int numberOfClusters, const std::string& outputFileName)
{
  std::cout << "WriteClusteredPixelsInRGBSpace()" << std::endl;
  ClusterColorsAdaptive clusterColors;
  clusterColors.SetNumberOfColors(numberOfClusters);
  clusterColors.ConstructFromImage(image);

  vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
  vtkSmartPointer<vtkUnsignedCharArray> colorsVTK = vtkSmartPointer<vtkUnsignedCharArray>::New();
  colorsVTK->SetName("Colors");
  colorsVTK->SetNumberOfComponents(3);

  vtkSmartPointer<vtkUnsignedIntArray> ids = vtkSmartPointer<vtkUnsignedIntArray>::New();
  ids->SetName("Ids");
  ids->SetNumberOfComponents(1);

  ColorMeasurementVectorType queryPoint;
  std::vector<ColorMeasurementVectorType> colors = clusterColors.GetColors();
  itk::ImageRegionConstIterator<FloatVectorImageType> imageIterator(image, image->GetLargestPossibleRegion());
  while(!imageIterator.IsAtEnd())
    {
    FloatVectorImageType::PixelType pixel = imageIterator.Get();
    queryPoint[0] = pixel[0];
    queryPoint[1] = pixel[1];
    queryPoint[2] = pixel[2];

    ClusterColors::TreeType::InstanceIdentifierVectorType neighbors;
    clusterColors.GetKDTree()->Search( queryPoint, 1u, neighbors );
    points->InsertNextPoint(pixel[0], pixel[1], pixel[2]);
    // TODO: Narrowing conversion warning
    //unsigned char color[3] = {colors[neighbors[0]][0], colors[neighbors[0]][1], colors[neighbors[0]][2]};
    unsigned char color[3]; // TODO: placeholder so it will compile
    colorsVTK->InsertNextTupleValue(color);
    ids->InsertNextValue(neighbors[0]);
    ++imageIterator;
    }

  vtkSmartPointer<vtkPolyData> polyData = vtkSmartPointer<vtkPolyData>::New();
  polyData->SetPoints(points);
  polyData->GetPointData()->SetScalars(colorsVTK);
  polyData->GetPointData()->AddArray(ids);

  vtkSmartPointer<vtkVertexGlyphFilter> glyphFilter = vtkSmartPointer<vtkVertexGlyphFilter>::New();
  glyphFilter->SetInputData(polyData);
  glyphFilter->Update();

  vtkSmartPointer<vtkXMLPolyDataWriter> writer = vtkSmartPointer<vtkXMLPolyDataWriter>::New();
  writer->SetInputConnection(glyphFilter->GetOutputPort());
  writer->SetFileName(outputFileName.c_str());
  writer->Write();

}

// Convert a vector ITK image to a VTK image for display
void ITKImagetoVTKRGBImage(FloatVectorImageType::Pointer image, vtkImageData* outputImage)
{
  // This function assumes an ND (with N>3) image has the first 3 channels as RGB and extra information in the remaining channels.
  
  //std::cout << "ITKImagetoVTKRGBImage()" << std::endl;
  if(image->GetNumberOfComponentsPerPixel() < 3)
    {
    std::stringstream ss;
    ss << "The input image has " << image->GetNumberOfComponentsPerPixel() << " components, but at least 3 are required.";
    throw std::runtime_error(ss.str());
    }

  // Setup and allocate the image data
  outputImage->SetNumberOfScalarComponents(3);
  outputImage->SetScalarTypeToUnsignedChar();
  outputImage->SetDimensions(image->GetLargestPossibleRegion().GetSize()[0],
                             image->GetLargestPossibleRegion().GetSize()[1],
                             1);

  outputImage->AllocateScalars();

  // Copy all of the input image pixels to the output image
  itk::ImageRegionConstIteratorWithIndex<FloatVectorImageType> imageIterator(image,image->GetLargestPossibleRegion());
  imageIterator.GoToBegin();

  while(!imageIterator.IsAtEnd())
    {
    unsigned char* pixel = static_cast<unsigned char*>(outputImage->GetScalarPointer(imageIterator.GetIndex()[0],
                                                                                     imageIterator.GetIndex()[1],0));
    for(unsigned int component = 0; component < 3; component++)
      {
      pixel[component] = static_cast<unsigned char>(imageIterator.Get()[component]);
      }

    ++imageIterator;
    }
}

// Convert a vector ITK image to a VTK image for display
void ITKImagetoVTKMagnitudeImage(FloatVectorImageType::Pointer image, vtkImageData* outputImage)
{
  std::cout << "ITKImagetoVTKMagnitudeImage()" << std::endl;
  // Compute the magnitude of the ITK image
  typedef itk::VectorMagnitudeImageFilter<
                  FloatVectorImageType, FloatScalarImageType >  VectorMagnitudeFilterType;

  // Create and setup a magnitude filter
  VectorMagnitudeFilterType::Pointer magnitudeFilter = VectorMagnitudeFilterType::New();
  magnitudeFilter->SetInput( image );
  magnitudeFilter->Update();

  // Rescale and cast for display
  typedef itk::RescaleIntensityImageFilter<
                  FloatScalarImageType, UnsignedCharScalarImageType > RescaleFilterType;

  RescaleFilterType::Pointer rescaleFilter = RescaleFilterType::New();
  rescaleFilter->SetOutputMinimum(0);
  rescaleFilter->SetOutputMaximum(255);
  rescaleFilter->SetInput( magnitudeFilter->GetOutput() );
  rescaleFilter->Update();

  // Setup and allocate the VTK image
  outputImage->SetNumberOfScalarComponents(1);
  outputImage->SetScalarTypeToUnsignedChar();
  outputImage->SetDimensions(image->GetLargestPossibleRegion().GetSize()[0],
                             image->GetLargestPossibleRegion().GetSize()[1],
                             1);

  outputImage->AllocateScalars();

  // Copy all of the scaled magnitudes to the output image
  itk::ImageRegionConstIteratorWithIndex<UnsignedCharScalarImageType> imageIterator(rescaleFilter->GetOutput(), rescaleFilter->GetOutput()->GetLargestPossibleRegion());
  imageIterator.GoToBegin();

  while(!imageIterator.IsAtEnd())
    {
    unsigned char* pixel = static_cast<unsigned char*>(outputImage->GetScalarPointer(imageIterator.GetIndex()[0],
                                                                                     imageIterator.GetIndex()[1],0));
    pixel[0] = imageIterator.Get();

    ++imageIterator;
    }
}

int main(int argc, char *argv[])
{
  // Verify arguments
  if(argc != 5)
    {
    std::cerr << "Required arguments: image imageMask patchRadius outputPrefix" << std::endl;
    return EXIT_FAILURE;
    }

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

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

  std::string outputPrefix = argv[4];

  // Output arguments
  std::cout << "Reading image: " << imageFilename << std::endl;
  std::cout << "Reading mask: " << maskFilename << std::endl;
  std::cout << "Patch radius: " << patchRadius << std::endl;
  //std::cout << "Output: " << outputFilename << std::endl;

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

  FloatVectorImageType::Pointer scaledImage = FloatVectorImageType::New();
  // Initialize
  Helpers::DeepCopy<FloatVectorImageType>(imageReader->GetOutput(), scaledImage);

  std::vector<float> maxValues = Helpers::MaxValuesVectorImage<float>(imageReader->GetOutput());

  // Scale all channels the same
//   for(unsigned int channel = 0; channel < imageReader->GetOutput()->GetNumberOfComponentsPerPixel(); ++channel)
//     {
//     Helpers::ScaleChannel<float>(imageReader->GetOutput(), channel, 1.0f, scaledImage);
//     }

  // Scale color channels
  for(unsigned int channel = 0; channel < 3; ++channel)
    {
    Helpers::ScaleChannel<float>(scaledImage, channel, 0.33f, scaledImage);
    }

  // Scale depth channel
  Helpers::ScaleChannel<float>(scaledImage, 3, 1.0f, scaledImage);

  Helpers::WriteImage<FloatVectorImageType>(scaledImage, "scaled.mha");

  typedef  itk::ImageFileReader< Mask > MaskReaderType;
  MaskReaderType::Pointer maskReader = MaskReaderType::New();
  maskReader->SetFileName(maskFilename.c_str());
  maskReader->Update();

  Mask::Pointer finalMask = Mask::New();
  ModifyMask(maskReader->GetOutput(), patchRadius, finalMask);

  cout.setf(ios::showpoint);

  std::vector<float> lambdas;
  for(unsigned int i = 0; i <= 10; ++i)
    {
    lambdas.push_back(0.1f * static_cast<float>(i));
    std::cout << "Using lambda " << lambdas[i] << std::endl;
    }

  std::shared_ptr<SelfPatchCompare> patchCompare(new SelfPatchCompare);
  patchCompare->SetNumberOfComponentsPerPixel(imageReader->GetOutput()->GetNumberOfComponentsPerPixel());
  //patchCompare->FunctionsToCompute.push_back(boost::bind(&SelfPatchCompare::SetPatchAverageAbsoluteSourceDifference,patchCompare,_1));
  patchCompare->FunctionsToCompute.push_back(boost::bind(&SelfPatchCompare::SetPatchColorDifference,patchCompare,_1));
  patchCompare->FunctionsToCompute.push_back(boost::bind(&SelfPatchCompare::SetPatchDepthDifference,patchCompare,_1));

  std::ofstream fout("scores.txt");

  fout.setf(ios::showpoint);
  for(unsigned int lambdaId = 0; lambdaId < lambdas.size(); ++lambdaId)
    {
    // Inpaint
    std::cout << "Inpainting with lambda = " << lambdas[lambdaId] << std::endl;

    PatchPair::DepthColorLambda = lambdas[lambdaId];

    CriminisiInpainting inpainting;
    //inpainting.SetDebugFunctionEnterLeave(true);
    inpainting.SetPatchRadius(patchRadius);
    inpainting.SetImage(scaledImage);
    inpainting.SetMask(finalMask);
    inpainting.SetMaxForwardLookPatches(3);
    inpainting.SetPatchCompare(patchCompare);

    inpainting.PatchSortFunction = &SortByDepthAndColor;
    //inpainting.PatchSortFunction = &SortByAverageAbsoluteDifference;

    //DepthAndColorDifference = ColorDifference * Lambda + (1.0 - Lambda) * DepthDifference;
    // When lambda = 0, only the depth is used
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
  unsigned int t = time(NULL);
  srand(t);

  itk::Size<2> size;
  size.Fill(100);

  itk::Index<2> index;
  index.Fill(0);

  itk::ImageRegion<2> region(index, size);
/*
  // Generate a random image (this method doesn't work with VectorImage)
  itk::RandomImageSource<FloatVectorImageType>::Pointer imageSource =
    itk::RandomImageSource<FloatVectorImageType>::New();
  imageSource->SetNumberOfThreads(1); // to produce non-random results
  imageSource->SetSize(size);
  imageSource->SetMin(0);
  imageSource->SetMax(100);
  imageSource->Update();
  FloatVectorImageType::Pointer image = imageSource->GetOutput();
*/
  // Generate a random image
  FloatVectorImageType::Pointer image = FloatVectorImageType::New();
  image->SetRegions(region);
  image->SetNumberOfComponentsPerPixel(3);
  image->Allocate();

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

  while(!imageIterator.IsAtEnd())
    {
    FloatVectorImageType::PixelType pixel;
    pixel.SetSize(3);
    pixel[0] = drand48();
    pixel[1] = drand48();
    pixel[2] = drand48();
    imageIterator.Set(pixel);
    ++imageIterator;
    }
  }

  // Generate a random membership image
  IntImageType::Pointer membershipImage = IntImageType::New();
  membershipImage->SetRegions(region);
  membershipImage->Allocate();
  membershipImage->FillBuffer(0);

  {
  itk::ImageRegionIterator<IntImageType> membershipImageIterator(membershipImage, membershipImage->GetLargestPossibleRegion());

  while(!membershipImageIterator.IsAtEnd())
    {
    IntImageType::PixelType pixel;
    pixel = rand() / 1000;
    membershipImageIterator.Set(pixel);
    ++membershipImageIterator;
    }
  }

  // Write the image
  itk::ImageFileWriter<FloatVectorImageType>::Pointer imageWriter =
    itk::ImageFileWriter<FloatVectorImageType>::New();
  imageWriter->SetFileName("image.mha");
  imageWriter->SetInput(image);
  imageWriter->Update();

//   // Generate a random mask
//   itk::RandomImageSource<Mask>::Pointer maskSource = itk::RandomImageSource<Mask>::New();
//   maskSource->SetNumberOfThreads(1); // to produce non-random results
//   maskSource->SetSize(size);
//   maskSource->SetMin(0);
//   maskSource->SetMax(255);
//   maskSource->Update();
//
//   // Threshold the mask
//   //typedef itk::ThresholdImageFilter <UnsignedCharImageType> ThresholdImageFilterType;
//   typedef itk::BinaryThresholdImageFilter <Mask, Mask> ThresholdImageFilterType;
//   ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
//   thresholdFilter->SetInput(maskSource->GetOutput());
//   thresholdFilter->SetLowerThreshold(0);
//   thresholdFilter->SetUpperThreshold(122);
//   thresholdFilter->SetOutsideValue(1);
//   thresholdFilter->SetInsideValue(0);
//   thresholdFilter->Update();
//   Mask::Pointer mask = thresholdFilter->GetOutput();

  std::cout << "Creating mask..." << std::endl;
  Mask::Pointer mask = Mask::New();
  mask->SetRegions(region);
  mask->Allocate();

  {
  itk::ImageRegionIterator<Mask> maskIterator(mask, mask->GetLargestPossibleRegion());

  while(!maskIterator.IsAtEnd())
    {
    int randomNumber = rand()%10;
//.........这里部分代码省略.........

void PartialPatchVectorComparison()
{
  std::cout << "PartialPatchVectorComparison()" << std::endl;

  const unsigned int dimension = 3;
  FloatVectorImageType::Pointer vectorImage = FloatVectorImageType::New();
  Testing::GetBlankImage<FloatVectorImageType>(vectorImage, dimension);

  // Make the left half of the image (0,0,0) and the right half (5,6,7)
  itk::ImageRegionIterator<FloatVectorImageType> imageIterator(vectorImage, vectorImage->GetLargestPossibleRegion());

  itk::VariableLengthVector<float> leftHalfPixel;
  leftHalfPixel.SetSize(dimension);
  leftHalfPixel.Fill(0);

  itk::VariableLengthVector<float> rightHalfPixel;
  rightHalfPixel.SetSize(dimension);
  rightHalfPixel[0] = 5;
  rightHalfPixel[1] = 6;
  rightHalfPixel[2] = 7;

  while(!imageIterator.IsAtEnd())
    {
    if(imageIterator.GetIndex()[0] < static_cast<int>(vectorImage->GetLargestPossibleRegion().GetSize()[0]/2))
      {
      imageIterator.Set(leftHalfPixel);
      }
    else
      {
      imageIterator.Set(rightHalfPixel);
      }

    ++imageIterator;
    }

  itk::Size<2> patchSize;
  patchSize.Fill(10);

  // Full patches differ
  std::cout << "Full patch different test." << std::endl;
  {
  itk::Index<2> sourceCorner;
  sourceCorner.Fill(0);
  itk::ImageRegion<2> sourceRegion(sourceCorner, patchSize);
  ImagePatchPixelDescriptor<FloatVectorImageType> sourcePatch(vectorImage, sourceRegion, true);

  itk::Index<2> targetCorner;
  targetCorner.Fill(vectorImage->GetLargestPossibleRegion().GetSize()[0]/2 + 4); // No magic about 4, just want a patch on the right side of the image
  itk::ImageRegion<2> targetRegion(targetCorner, patchSize);
  ImagePatchPixelDescriptor<FloatVectorImageType> targetPatch(vectorImage, targetRegion, true);

  PatchPair<FloatVectorImageType> patchPair(&sourcePatch, targetPatch);
  PatchDifferencePixelWiseSum<FloatVectorImageType, PixelDifference> vector_patchDifferencePixelWiseSum;
  vector_patchDifferencePixelWiseSum.SetImage(vectorImage);
  float difference = vector_patchDifferencePixelWiseSum.Difference(patchPair);

  float correctDifference = targetRegion.GetNumberOfPixels() * 18; // 18 = 5+6+7, the sum of the elements of 'rightHalfPixel'
  if(difference != correctDifference)
    {
    std::stringstream ss;
    ss << "Difference " << difference << " does not match correctDifference " << correctDifference;
    throw std::runtime_error(ss.str());
    }
  }

  // Full patches identical
  std::cout << "Identical patch test." << std::endl;
  {
  itk::Index<2> sourceCorner;
  sourceCorner.Fill(5);
  itk::ImageRegion<2> sourceRegion(sourceCorner, patchSize);
  ImagePatchPixelDescriptor<FloatVectorImageType> sourcePatch(vectorImage, sourceRegion, true);

  itk::Index<2> targetCorner;
  targetCorner.Fill(5);
  itk::ImageRegion<2> targetRegion(targetCorner, patchSize);
  ImagePatchPixelDescriptor<FloatVectorImageType> targetPatch(vectorImage, targetRegion, true);

  PatchPair<FloatVectorImageType> patchPair(&sourcePatch, targetPatch);
  PatchDifferencePixelWiseSum<FloatVectorImageType, PixelDifference> vector_patchDifferencePixelWiseSum;
  vector_patchDifferencePixelWiseSum.SetImage(vectorImage);
  float difference = vector_patchDifferencePixelWiseSum.Difference(patchPair);

  float correctDifference = 0;
  if(difference != correctDifference)
    {
    std::stringstream ss;
    ss << "Difference " << difference << " does not match correctDifference " << correctDifference;
    throw std::runtime_error(ss.str());
    }
  }
}

bool TestHasBracketOperator()
{
  {
  typedef itk::CovariantVector<int, 3> IntVectorType;

  static_assert(Helpers::HasBracketOperator<IntVectorType>::value,
              "TestHasBracketOperator for CovariantVector<int, 3> failed!");
  }

  {
  typedef itk::CovariantVector<unsigned char, 3> UCharVectorType;

  static_assert(Helpers::HasBracketOperator<UCharVectorType>::value,
              "TestHasBracketOperator for CovariantVector<unsigned char, 3> failed!");
  }

  {
  typedef itk::CovariantVector<float, 3> FloatVectorType;

  static_assert(Helpers::HasBracketOperator<FloatVectorType>::value,
              "TestHasBracketOperator for CovariantVector<float, 3> failed!");
  }

  {
  typedef itk::Image<itk::CovariantVector<float, 3>, 2> FloatVectorImageType;

  static_assert(Helpers::HasBracketOperator<FloatVectorImageType::PixelType>::value,
              "TestHasBracketOperator for CovariantVector<float, 3> failed!");
  }

  {
  typedef itk::Image<itk::CovariantVector<float, 3>, 2> FloatVectorImageType;
  FloatVectorImageType::Pointer image = FloatVectorImageType::New();

  TestHasBracketOperator_Template(image.GetPointer());
  }

  {
  typedef itk::Image<itk::CovariantVector<float, 3>, 2> FloatVectorImageType;
  itk::SmartPointer<FloatVectorImageType> image = FloatVectorImageType::New();

  TestHasBracketOperator_Template(image.GetPointer());
  TestHasBracketOperator_ConstTemplate(image.GetPointer());
  }

  {
  typedef itk::Image<itk::CovariantVector<float, 3>, 2> FloatVectorImageType;

  TestHasBracketOperator_Template2<FloatVectorImageType>();
  }

  {
  typedef itk::VariableLengthVector<int> VectorType;

  static_assert(Helpers::HasBracketOperator<VectorType>::value,
              "TestHasBracketOperator for VariableLengthVector failed!");
  }

  {
  typedef std::vector<int> VectorType;

  static_assert(Helpers::HasBracketOperator<VectorType>::value,
              "TestHasBracketOperator for std::vector failed!");
  }

  // This (intentionally) fails
//  {
//  static_assert(Helpers::HasBracketOperator<float>::value,
//              "TestHasBracketOperator for float failed!");
//  }

  return true;
}

int main(int argc, char *argv[])
{
  unsigned int t = time(NULL);
  srand(t);
  
  itk::Size<2> size;
  size.Fill(10);

  itk::Index<2> index;
  index.Fill(0);

  itk::ImageRegion<2> region(index, size);
/*
  // Generate a random image (this method doesn't work with VectorImage)
  itk::RandomImageSource<FloatVectorImageType>::Pointer imageSource =
    itk::RandomImageSource<FloatVectorImageType>::New();
  imageSource->SetNumberOfThreads(1); // to produce non-random results
  imageSource->SetSize(size);
  imageSource->SetMin(0);
  imageSource->SetMax(100);
  imageSource->Update();
  FloatVectorImageType::Pointer image = imageSource->GetOutput();
*/
  // Generate a random image
  FloatVectorImageType::Pointer image = FloatVectorImageType::New();
  image->SetRegions(region);
  image->SetNumberOfComponentsPerPixel(3);
  image->Allocate();

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

  while(!imageIterator.IsAtEnd())
    {
    FloatVectorImageType::PixelType pixel;
    pixel.SetSize(3);
    pixel[0] = RandomFloat();
    pixel[1] = RandomFloat();
    pixel[2] = RandomFloat();
    imageIterator.Set(pixel);
    ++imageIterator;
    }
  }
  
  // Write the image
  itk::ImageFileWriter<FloatVectorImageType>::Pointer imageWriter =
    itk::ImageFileWriter<FloatVectorImageType>::New();
  imageWriter->SetFileName("image.mha");
  imageWriter->SetInput(image);
  imageWriter->Update();

  // Generate a random mask
  itk::RandomImageSource<Mask>::Pointer maskSource =
    itk::RandomImageSource<Mask>::New();
  maskSource->SetNumberOfThreads(1); // to produce non-random results
  maskSource->SetSize(size);
  maskSource->SetMin(0);
  maskSource->SetMax(255);
  maskSource->Update();

  // Threshold the mask
  //typedef itk::ThresholdImageFilter <UnsignedCharImageType> ThresholdImageFilterType;
  typedef itk::BinaryThresholdImageFilter <Mask, Mask> ThresholdImageFilterType;
  ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
  thresholdFilter->SetInput(maskSource->GetOutput());
  thresholdFilter->SetLowerThreshold(0);
  thresholdFilter->SetUpperThreshold(122);
  thresholdFilter->SetOutsideValue(1);
  thresholdFilter->SetInsideValue(0);
  thresholdFilter->Update();
  Mask::Pointer mask = thresholdFilter->GetOutput();

  // Write the mask
  itk::ImageFileWriter<Mask>::Pointer maskWriter =
    itk::ImageFileWriter<Mask>::New();
  maskWriter->SetFileName("mask.png");
  maskWriter->SetInput(mask);
  maskWriter->Update();

  unsigned int patchRadius = 1;
  // Create source patches
  itk::ImageRegionConstIterator<FloatVectorImageType> imageIterator(image, image->GetLargestPossibleRegion());
  std::vector<itk::ImageRegion<2> > sourcePatches;
  while(!imageIterator.IsAtEnd())
    {
    itk::Index<2> currentPixel = imageIterator.GetIndex();
    itk::ImageRegion<2> region = GetRegionInRadiusAroundPixel(currentPixel, patchRadius);
    if(image->GetLargestPossibleRegion().IsInside(region))
      {
      sourcePatches.push_back(region);
      }
    ++imageIterator;
    }

  itk::Size<2> targetSize;
  targetSize.Fill(patchRadius * 2 + 1);

  itk::Index<2> targetIndex;
  targetIndex.Fill(3);
  
//.........这里部分代码省略.........