text: /*=========================================================================
   *
   *  Copyright Insight Software Consortium
   *
   *  Licensed under the Apache License, Version 2.0 (the "License");
   *  you may not use this file except in compliance with the License.
   *  You may obtain a copy of the License at
   *
   *         http://www.apache.org/licenses/LICENSE-2.0.txt
   *
   *  Unless required by applicable law or agreed to in writing, software
   *  distributed under the License is distributed on an "AS IS" BASIS,
   *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   *  See the License for the specific language governing permissions and
   *  limitations under the License.
   *
   *=========================================================================*/
  
  #include <iostream>
  
  #include "itkDeformableMesh3DFilter.h"
  #include "itkGradientRecursiveGaussianImageFilter.h"
  #include "itkGradientImageFilter.h"
  #include "itkGradientToMagnitudeImageFilter.h"
  #include "itkGradientVectorFlowImageFilter.h"
  #include "itkShrinkImageFilter.h"
  #include "itkBinaryMask3DMeshSource.h"
  
  int itkDeformableTest(int , char *[])
  {
    const itk::Size<3u>::SizeValueType WIDTH = 32;
    const itk::Size<3u>::SizeValueType HEIGHT = 32;
    const itk::Size<3u>::SizeValueType DEPTH = 32;
  
    // Define the dimension of the images
    const unsigned int myDimension = 3;
  
    // Declare the types of the images
    typedef itk::Image<double, myDimension>           myImageType;
  
    // Declare the types of the output images
    typedef itk::Image<unsigned short, myDimension>   binaryImageType;
  
    // Declare the type of the index to access images
    typedef itk::Index<myDimension>       myIndexType;
  
    // Declare the type of the size
    typedef itk::Size<myDimension>        mySizeType;
  
    // Declare the type of the Region
    typedef itk::ImageRegion<myDimension> myRegionType;
  
    // Declare the type of the Mesh
    typedef itk::Mesh<double>   DMesh;
  
    // Declare the type of the gradient image
    typedef itk::CovariantVector<double, myDimension> myGradientType;
    typedef itk::Image<myGradientType, myDimension>   myGradientImageType;
    typedef itk::CovariantVector<double, 2>           double2DVector;
    typedef itk::CovariantVector<double, 3>           double3DVector;
  
    typedef itk::BinaryMask3DMeshSource<binaryImageType,DMesh>  myMeshSource;
    typedef itk::LaplacianImageFilter<myImageType, myImageType> myLaplacianFilterType;
    typedef itk::GradientVectorFlowImageFilter<myGradientImageType, myGradientImageType>
                                                myGVFFilterType;
  
    typedef itk::GradientImageFilter<myImageType, double, double>
                                                myGFilterType;
  
    typedef itk::GradientToMagnitudeImageFilter<myGradientImageType, myImageType>
                                                myGToMFilterType;
  
    typedef itk::DeformableMesh3DFilter<DMesh, DMesh>  DFilter;
  
    binaryImageType::Pointer       biimg=binaryImageType::New();
    myGradientImageType::Pointer   gdimg=myGradientImageType::New();
  
    binaryImageType::SizeType      bisize={{WIDTH,HEIGHT,DEPTH}};
    binaryImageType::IndexType     biindex;
    binaryImageType::RegionType    biregion;
    biindex.Fill(0);
    biregion.SetSize(bisize);
    biregion.SetIndex(biindex);
  
    myGradientImageType::SizeType   gdsize={{WIDTH,HEIGHT,DEPTH}};
    myGradientImageType::IndexType  gdindex;
    myGradientImageType::RegionType gdregion;
    gdindex.Fill(0);
    gdregion.SetSize(gdsize);
    gdregion.SetIndex(gdindex);
  
    biimg->SetLargestPossibleRegion( biregion );
    biimg->SetBufferedRegion( biregion );
    biimg->SetRequestedRegion( biregion );
    biimg->Allocate();
  
    gdimg->SetLargestPossibleRegion( gdregion );
    gdimg->SetBufferedRegion( gdregion );
    gdimg->SetRequestedRegion( gdregion );
    gdimg->Allocate();
  
    // Create the image
    myImageType::Pointer inputImage  = myImageType::New();
  
    mySizeType size={{WIDTH,HEIGHT,DEPTH}};
    myIndexType start;
    start.Fill(0);
  
    myRegionType region;
    region.SetIndex( start );
    region.SetSize( size );
  
    // Initialize Image A
    inputImage->SetLargestPossibleRegion( region );
    inputImage->SetBufferedRegion( region );
    inputImage->SetRequestedRegion( region );
    inputImage->Allocate();
  
    itk::ImageRegionIteratorWithIndex <myImageType> it(inputImage, region);
    it.GoToBegin();
    itk::ImageRegionIteratorWithIndex <binaryImageType> bit(biimg, biregion);
    bit.GoToBegin();
  
    /////////////////////////////////////////////////////////////////////////
  
  
    while( !it.IsAtEnd() )
    {
      it.Set( 0.0 );
      bit.Set( 0 );
      ++it;
      ++bit;
    }
  
    size[0] = 16;
    size[1] = 16;
    size[2] = 16;
  
    start[0] = 8;
    start[1] = 8;
    start[2] = 8;
  
    // Create one iterator for an internal region
    region.SetSize( size );
    region.SetIndex( start );
    biregion.SetSize( size );
    biregion.SetIndex( start );
    itk::ImageRegionIteratorWithIndex <myImageType> itb( inputImage, region );
    itk::ImageRegionIteratorWithIndex <binaryImageType> bitb( biimg, biregion );
  
    // Initialize the content the internal region
    while( !itb.IsAtEnd() )
    {
      itb.Set( 100.0 );
      bitb.Set ( 255 );
      ++itb;
      ++bitb;
    }
  
  
    //////////////////////////////////////////////////////////////////////////
  
    itk::ShrinkImageFilter< myImageType, myImageType >::Pointer dshrink;
    dshrink = itk::ShrinkImageFilter< myImageType, myImageType >::New();
    dshrink->SetInput( inputImage );
    dshrink->SetNumberOfThreads(4);
  
    unsigned int dfactors[3] = { 1, 1, 1 };
    dshrink->SetShrinkFactors(dfactors);
    dshrink->UpdateLargestPossibleRegion();
  
    myImageType::RegionType drequestedRegion;
    drequestedRegion = dshrink->GetOutput()->GetRequestedRegion();
  
    typedef itk::GradientRecursiveGaussianImageFilter<
                                              myImageType,
                                              myGradientImageType
                                              >  myFilterType;
  
  
    // Create a  Filter
    myFilterType::Pointer grfilter = myFilterType::New();
    myGFilterType::Pointer gfilter = myGFilterType::New();
    myGToMFilterType::Pointer gtomfilter = myGToMFilterType::New();
  
    // Connect the input images
    grfilter->SetInput( dshrink->GetOutput() );
  
    // Set sigma
    grfilter->SetSigma( 1.0 );
  
    myLaplacianFilterType::Pointer m_LFilter = myLaplacianFilterType::New();
    myGVFFilterType::Pointer m_GVFFilter = myGVFFilterType::New();
  
    m_GVFFilter->SetInput(gfilter->GetOutput());
    m_GVFFilter->SetLaplacianFilter(m_LFilter);
    m_GVFFilter->SetNoiseLevel(500);
  
    gtomfilter->SetInput(grfilter->GetOutput());
  
    gfilter->SetInput(gtomfilter->GetOutput());
    gfilter->Update();
  
    std::cout << "The gradient map created!" << std::endl;
  
  //  the gvf is temproraily disabled since the problem related with gradientimagefilter
  //  m_GVFFilter->Update();
  
  //  std::cout << "GVF created! " << std::endl;
  
  ////////////////////////////////////////////////////////////////////////////////////////
  // construct the deformable mesh
  
    itk::ShrinkImageFilter< binaryImageType, binaryImageType >::Pointer shrink;
    shrink = itk::ShrinkImageFilter< binaryImageType, binaryImageType >::New();
    shrink->SetInput( biimg );
    shrink->SetNumberOfThreads(4);
  
    unsigned int factors[3] = { 1, 1, 1 };
    shrink->SetShrinkFactors(factors);
    shrink->UpdateLargestPossibleRegion();
  
    binaryImageType::RegionType requestedRegion;
    requestedRegion = shrink->GetOutput()->GetRequestedRegion();
  
    myMeshSource::Pointer m_bmmeshsource = myMeshSource::New();
  
    DFilter::Pointer m_dfilter = DFilter::New();
    m_dfilter->SetInput(m_bmmeshsource->GetOutput());
  //  m_dfilter->SetGradient(m_GVFFilter->GetOutput());
    m_dfilter->SetGradient(gfilter->GetOutput());
  
    m_bmmeshsource->SetInput( shrink->GetOutput() );
    m_bmmeshsource->SetObjectValue( 255 );
    m_bmmeshsource->Update();
  
    std::cout << "Deformable mesh created using Marching Cube!" << std::endl;
  
    double2DVector m_stiff;
    m_stiff[0] = 0.0001;
    m_stiff[1] = 0.1;
  
    double3DVector m_scale;
    m_scale[0] = 1;
    m_scale[1] = 1;
    m_scale[2] = 1;
    m_dfilter->SetStiffness(m_stiff);
    m_dfilter->SetGradientMagnitude(1.0);
    m_dfilter->SetTimeStep(0.01);
    m_dfilter->SetStepThreshold(10);
    m_dfilter->SetScale(m_scale);
    m_dfilter->SetObjectLabel(1);
    m_dfilter->SetPotentialOn(0);
    m_dfilter->SetPotentialMagnitude(1.0);
  
    std::cout << "Deformable mesh fitting...";
    m_dfilter->Update();
  
    /** raise coverage */
    myGradientImageType::Pointer grad_tmp = m_dfilter->GetGradient();
    std::cout << grad_tmp << std::endl;
    std::cout << m_dfilter->GetStepThreshold() << std::endl;
  
    double2DVector stiff_tmp = m_dfilter->GetStiffness();
    std::cout << stiff_tmp <<std::endl;
  
    std::cout << m_dfilter->GetTimeStep() << std::endl;
  
    DMesh::Pointer norm_tmp = m_dfilter->GetNormals();
    std::cout << norm_tmp << std::endl;
    std::cout << m_dfilter;
  
    std::cout << "Mesh Source: " << m_bmmeshsource;
  
  // All objects should be automatically destroyed at this point
  
    std::cout << "[TEST DONE]" << std::endl;
    return EXIT_SUCCESS;
  
  }
  