[Insight-users] Operator '=' for Mesh

Luis Ibanez luis.ibanez at kitware.com
Mon Oct 4 16:21:28 EDT 2004


Hi Uday,

Thanks for posting the details of your code.

The problem is that you are inter-mixing two
different declarations of the MeshType.

A) In one part you use

    typedef itk::DefaultStaticMeshTraits<
              double, 3, 3,
              double,double
                  > TriangleMeshTraits;

    typedef itk::Mesh<double,3,
             TriangleMeshTraits> TriangleMeshType;


B) and later on, you define *another* mesh type as

  typedef double                            MeshPixelType;
  typedef itk::Mesh< MeshPixelType >        MeshType;




These two types (A) and (B) are incompatible. The second
is using the default template parameters defined for the
Mesh class.


You instantiate a SphereMeshSource using MeshType and
then try to assign its output to a smartPointer of
TriangleMeshType.

Please choose *only one* of the two mesh type and be
consistent all along your code.


Attached is a modified version that compiles on Linux.


   Regards,


     Luis



--------------------
Uday Kurkure wrote:
> Hello Luis,
> 
> Thanks for replying. I am using smart pointer. I am
> declaring it first and then using it. I am attaching
> the code. (The code is same as DeformableModel2.cxx, I
> am just adding mesh visualization to it)
> 
> Please see what I am missing.
> 
> Thanks,
> 
> -Uday
> 
> 
> 
> --- Luis Ibanez <luis.ibanez at kitware.com> wrote:
> 
> 
>>Hi Uday,
>>
>>You are probably trying to assing a Mesh raw pointer
>>to a variable of incompatible type.
>>
>>Please read the "Data Representation" Chapter of the
>>ITK Software Guide
>>
>>     http://www.itk.org/ItkSoftwareGuide.pdf
>>
>>Chapter 4, Sections 4.2, 4.3, pages pdf- 74 to 116.
>>
>>You should receive the Mesh pointer into a smart
>>pointer. Something like
>>
>>  MeshType::Pointer myMesh =
>>SphericalMesh->GetOutput();
>>
>>
>>For a description of the principles behind the use
>>of SmartPointers please look at the Tutorials
>>
>>     http://www.itk.org/HTML/Tutorials.htm
>>
>>in particular to
>>
>>
> 
> http://www.itk.org/CourseWare/Training/GettingStarted-I.pdf
> 
> http://www.itk.org/CourseWare/Training/GettingStartedI-WebPage/index.htm
> 
>>
>>
>>Regards,
>>
>>
>>    Luis
>>
>>
>>---------------------
>>Uday Kurkure wrote:
>>
>>>Hi All,
>>>
>>>I have some problem with visualization of mesh. I
>>>generate a spherical mesh and assign its output to
>>
>>a
>>
>>>mesh pointer:
>>>
>>>myMesh = SphericalMesh->GetOutput();
>>>(Rest of the code converts itkmesh to vtkpolydata)
>>>
>>>But the code is not compiling. I saw such
>>>implementation in some applications but I could
>>
>>not
>>
>>>find anything that I am missing. Can some one tell
>>
>>me
>>
>>>why there is a compiler error at this line.
>>>
>>>Error: 
>>>error C2679: binary '=' : no operator defined
>>
>>which
>>
>>>takes a right-hand operand of type 'class
>>>itk::Mesh<double,3,class
>>>itk::DefaultStaticMeshTraits<double,3,3,float,floa
>>>t,double> > *' (or there is no acceptable
>>
>>conversion)
>>
>>>Thanks,
>>>-Uday
>>>
>>>
>>>	
=========================================================



-------------- next part --------------
/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: DeformableModel2.cxx,v $
  Language:  C++
  Date:      $Date: 2004/04/03 13:00:22 $
  Version:   $Revision: 1.2 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
// Disable warning for long symbol names in this file only
#ifdef _MSC_VER
#pragma warning ( disable : 4786 )
#endif

// Software Guide : BeginLatex
//
// This example illustrates the use of the \doxygen{DeformableMesh3DFilter}.
// An initial mesh is created using the \doxygen{SphereMeshSource} filter.
//
// \index{Deformable Models}
// \index{DeformableMesh3DFilter}
// \index{SphereMeshSource}
//
// Software Guide : EndLatex



// Software Guide : BeginCodeSnippet
#include "itkMesh.h"

#include "itkDeformableMesh3DFilter.h"

#include "itkGradientRecursiveGaussianImageFilter.h"
#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"

#include "itkImage.h"
#include "itkCovariantVector.h"

#include "itkSphereMeshSource.h"

#include "itkImageFileReader.h"

#include "vtkPolyData.h"
#include "vtkPolyDataMapper.h"
#include "vtkPoints.h"
#include <itkTriangleCell.h>
#include <itkCellInterface.h>
#include <vtkCellArray.h>
#include "vtkRenderer.h"
#include "vtkRenderWindow.h"
#include "vtkActor.h"
#include "vtkRenderWindowInteractor.h"

#ifndef vtkDoubleType
#define vtkDoubleType double
#endif


  typedef itk::DefaultStaticMeshTraits<double, 3, 3,double,double> TriangleMeshTraits;

  typedef itk::Mesh<double,3, TriangleMeshTraits> TriangleMeshType;


	class VistVTKCellsClass
	{
  vtkCellArray* m_Cells;
  int* m_LastCell;
  int* m_TypeArray;


  public:



  //typedef the itk cells we are interested in
  typedef itk::CellInterface< TriangleMeshType::PixelType,
                              TriangleMeshType::CellTraits >  CellInterfaceType;

  typedef itk::VertexCell<CellInterfaceType>                 VertexCell;
  typedef itk::LineCell<CellInterfaceType>                   LineCell;

  typedef itk::TriangleCell<CellInterfaceType>      TriangleCell;

  /*!
    Set the vtkCellArray that will be constructed
   */
  void SetCellArray(vtkCellArray* a)
    {
      m_Cells = a;
    }

  /*!
    Set the cell counter pointer
  */
  void SetCellCounter(int* i)
    {
      m_LastCell = i;
    }

  /*!
    Set the type array for storing the vtk cell types
  */
  void SetTypeArray(int* i)
    {
      m_TypeArray = i;
    }

   void Visit(unsigned long, TriangleCell* t)
    {
      m_Cells->InsertNextCell(3,  (vtkIdType*)t->PointIdsBegin());
      m_TypeArray[*m_LastCell] = VTK_TRIANGLE;
      (*m_LastCell)++;
    }


};

	typedef itk::CellInterfaceVisitorImplementation< TriangleMeshType::PixelType,
                                                   TriangleMeshType::CellTraits,

  itk::TriangleCell< itk::CellInterface<TriangleMeshType::PixelType, TriangleMeshType::CellTraits > >,
  VistVTKCellsClass> TriangleVisitor;

//----------------------------------------------------------------

// DISPLAYS VTK MESH

void DisplayMesh(vtkPolyData* itkgrid)
{
// Create the renderer and window stuff
  vtkRenderer* ren1 = vtkRenderer::New();
  vtkRenderWindow* renWin = vtkRenderWindow::New();
  renWin->AddRenderer(ren1);
  vtkRenderWindowInteractor* inter = vtkRenderWindowInteractor::New();
  inter->SetRenderWindow(renWin);

  vtkPolyDataMapper* mapper = vtkPolyDataMapper::New();
  mapper->SetInput(itkgrid);
  vtkActor* actor = vtkActor::New();
  actor->SetMapper(mapper);
  ren1->SetViewport(0.0, 0.0, 0.5, 1.0);
  // add the actor and start the render loop
  ren1->AddActor(actor);
  renWin->Render();
  inter->Start();

  mapper->Delete();
  actor->Delete();
  ren1->Delete();
  renWin->Delete();
}

//---------------------------------------------------------------

// CONVERTS INTO VTK POLYDATA AND DISPLAYS IT

int MeshVisualization(TriangleMeshType *m_MeshToShow)
{

	int numPoints =  m_MeshToShow->GetNumberOfPoints();

/*  if (numPoints == 0)
    {
      fl_alert( "no points in Grid ");
      return;
    }
*/
 // fill in here the conversion between itkMesh versus vtkUnstructuredGrid

  vtkPolyData* vgrid = vtkPolyData::New();

  // Create the vtkPoints object and set the number of points
  vtkPoints* vpoints = vtkPoints::New();
  vpoints->SetNumberOfPoints(numPoints);

  // iterate over all the points in the itk mesh filling in
  // the vtkPoints object as we go
  TriangleMeshType::PointsContainer::Pointer points = m_MeshToShow->GetPoints();
  for(TriangleMeshType::PointsContainer::Iterator i = points->Begin();
      i != points->End(); ++i)
    {
      // Get the point index from the point container iterator
      int idx = i->Index();
      // Set the vtk point at the index with the the coord array from itk
      // itk returns a const pointer, but vtk is not const correct, so
      // we have to use a const cast to get rid of the const
      vtkDoubleType * pp = const_cast<vtkDoubleType*>(i->Value().GetDataPointer());
      vpoints->SetPoint(idx, pp);
    }

  // Set the points on the vtk grid
  vgrid->SetPoints(vpoints);

  TriangleMeshType::CellType::MultiVisitor::Pointer mv =
    TriangleMeshType::CellType::MultiVisitor::New();

  TriangleVisitor::Pointer tv = TriangleVisitor::New();

  //set up the visitors
  int vtkCellCount = 0; // running counter for current cell inserted into vtk
  int numCells = m_MeshToShow->GetNumberOfCells();
  int *types = new int[numCells]; //type array for vtk

  //create vtk cells and estimate the size
  vtkCellArray* cells = vtkCellArray::New();
  cells->EstimateSize(numCells, 4);

  // Set the TypeArray CellCount and CellArray for both visitors
  tv->SetTypeArray(types);
  tv->SetCellCounter(&vtkCellCount);
  tv->SetCellArray(cells);

  mv->AddVisitor(tv);

  m_MeshToShow->Accept(mv);


  vgrid->SetLines(cells);

  cells->Delete();
  vpoints->Delete();
	DisplayMesh(vgrid);
}

//---------------------------------------------------------------------

// MAIN

int main()
{

  const char *inputFileName      = "Z:\\ITK\\Insight\\data\\sphere.mhd";
  const float sigma              = 5;
  const int   numberOfIterations = 100;
  const float timeStep           = 0.1;
  const float externalForceScale = 10;
  const float stiffness          = 0.1;



//  typedef double                            MeshPixelType;
//  typedef itk::Mesh< MeshPixelType >        MeshType;


  unsigned const int Dimension = 3;

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


  typedef itk::GradientRecursiveGaussianImageFilter<
                                        ImageType
                                           > GradientFilterType;

  typedef itk::GradientMagnitudeRecursiveGaussianImageFilter<
                                        ImageType,ImageType>
                                                      GradientMagnitudeFilterType;


  typedef itk::DeformableMesh3DFilter<TriangleMeshType,TriangleMeshType>  DeformableFilterType;


  typedef itk::ImageFileReader< ImageType       >  ReaderType;

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

  imageReader->SetFileName( inputFileName );

  GradientMagnitudeFilterType::Pointer  gradientMagnitudeFilter
                      = GradientMagnitudeFilterType::New();

  try
    {
    imageReader->Update();
    }
  catch( itk::ExceptionObject & excep )
    {
    std::cerr << "Problem found while reading file " << inputFileName << std::endl;
    std::cerr << "Exception Caught !" << std::endl;
    std::cerr << excep << std::endl;
    }


  ImageType::ConstPointer inputImage = imageReader->GetOutput();
  std::cout << "read the image" << std::endl;

  gradientMagnitudeFilter->SetInput( inputImage );


  gradientMagnitudeFilter->SetSigma( sigma );


  GradientFilterType::Pointer gradientMapFilter = GradientFilterType::New();


  gradientMapFilter->SetInput( gradientMagnitudeFilter->GetOutput());
  gradientMapFilter->SetSigma( sigma );

  gradientMapFilter->Update();
  std::cout << "gradient computed" << std::endl;


  DeformableFilterType::Pointer deformableModelFilter =
                                     DeformableFilterType::New();



  typedef itk::SphereMeshSource< TriangleMeshType >        MeshSourceType;

  MeshSourceType::Pointer meshSource = MeshSourceType::New();


//  MESH POINTER
//  TriangleMeshType::Pointer meshToShow;// = TriangleMeshType::New();
  TriangleMeshType::Pointer meshToShow;


  // Set the initial sphere in the center of the image
  //
  const ImageType::SpacingType spacing = inputImage->GetSpacing();
  ImageType::PointType         origin  = inputImage->GetOrigin();
  ImageType::SizeType          size    = inputImage->GetBufferedRegion().GetSize();


  TriangleMeshType::PointType center;
  center[0] = origin[0] + spacing[0] * size[0] / 2.0;
  center[1] = origin[1] + spacing[1] * size[1] / 2.0;
  center[2] = origin[2] + spacing[2] * size[2] / 2.0;
  meshSource->SetCenter( center );

  TriangleMeshType::PointType radius;
  radius[0] = spacing[0] * size[0] / 4.0;
  radius[1] = spacing[1] * size[1] / 4.0;
  radius[2] = spacing[2] * size[2] / 4.0;
  meshSource->SetScale( radius );


  meshSource->SetResolutionX( 50 );
  meshSource->SetResolutionY( 50 );
  meshSource->Update();
  std::cout << " mesh generated" << std::endl;

  deformableModelFilter->SetInput(    meshSource->GetOutput()        );

  deformableModelFilter->SetGradient( gradientMapFilter->GetOutput() );



  typedef itk::CovariantVector<double, 2>           StiffnessType;

  StiffnessType stiffnessVector;
  stiffnessVector[0] = stiffness;
  stiffnessVector[1] = stiffness;



  deformableModelFilter->SetTimeStep( timeStep );
  deformableModelFilter->SetStiffness( stiffnessVector );
  deformableModelFilter->SetStepThreshold( numberOfIterations );
  deformableModelFilter->SetGradientMagnitude( externalForceScale );

  std::cout << "deformation in progress" << std::endl;
  try
    {
    deformableModelFilter->Update();
    }
  catch( itk::ExceptionObject & excep )
    {
    std::cerr << "Exception Caught !" << std::endl;
    std::cerr << excep << std::endl;
    }

  std::cout << "[TEST DONE]" << std::endl;


// MESH VISUALIZATION

	//	meshToShow = deformableModelFilter->GetOutput();
	meshToShow = meshSource->GetOutput();
	MeshVisualization(meshToShow);
	return 0;
}





More information about the Insight-users mailing list