CGAL users discussion list

[Sascha Zelzer <>]

Sure,

I have attached a file with my code. You can ignore the types in the 
"mitk" namespace, or replace them with CGAL types where appropriate.

In my case, I have a labelled image as the input for the CGAL mesher for 
which I can provide a mapping of labels to a sampling distance for the 
seed points (the SeedPointDistanceToLabelMap type in the attached code). 
So for example I know that my image contains three labels, with 
consecutive values from 1 to 3. For value 1 I need a sampling distance 
of 15mm (for example), for value 2 5mm, etc. I then simply iterate over 
all different sampling distances and for each distance, I iterate over 
the whole image (with stepping according to the current distance) and 
create a seed point if the current image value corresponds to a label 
value for which that sampling distance is required.

The created seed points are then inserted into the triangluation object 
of my domain (instead of relying on the default constructed seed points 
in make_mesh_3). I'm not sure if that approach has any disadvantages 
(apart from brute force iterating over the image) but it seems to work okay.

- Sascha

On 09/28/2012 05:47 PM, Hamid G wrote:
> Sascha,
>
> Would you mind elaborating on what exactly you did? I have same
> problem and using varying sizing field didn't help me either.
>
> thanks in advance,
>
> On Mon, Sep 17, 2012 at 3:37 PM, Sascha Zelzer
> <>
>  wrote:
> > Hi,
> >
> > I was finally able to construct a suitable set of seed points and run the
> > mesh generation step successfully.
> >
> > Many thanks for your hints!
> >
> > Before my first posting to this list, I have tried using a CGAL::SizingField
> > and setting different values for the tet radii for the different
> > sub-domains. That did obviously not work out and I now know it is because of
> > an inappropriate default seed point set.
> >
> > Thanks,
> >
> > Sascha
> >
> >
> > On 09/10/2012 09:57 AM, Mariette Yvinec wrote:
> >
> > Currently,
> > the best way to ensure that small component are mot missed
> > is to set initial points in them or on them boundary.
> > Then you have to use refine_mesh_3 instead of
> > make_mesh_3 for your initaial seeding points to be taken into account.
> >
> > Another way, if you know the area where the small components are supposed to
> > be,
> > is to define a  variable sizing field with a value locally similar to the
> > size of the components
> > to be meshed.
> >
> >
> >
> > Le 30/08/12 16:18, Sascha Zelzer a écrit :
> >
> > Hi,
> >
> > I have posted about my problem already here
> > http://cgal-discuss.949826.n4.nabble.com/3D-volume-mesher-ignoring-small-subdomains-td4655756.html
> > but thought I reformulate it:
> >
> > I have a labeled image as input to make_mesh_3 which contains a small
> > sub-domain (label 2) inside a larger domain (label 1). Small meaning the its
> > characteristic size is about two orders of magnitude smaller then the
> > embedding domain.
> >
> > So my question is, what are the criteria for refining cells in the large
> > domain to ensure that they do not "shadow" a complete small sub-domain? I
> > would like to be able to set a large maximum cell size for the large domain
> > and a suitable cell size for small sub-domain. When setting the cell and
> > facet criteria in such a way, my small sub-domain is ignored however.
> >
> > Any tips would be highly appreciated.
> >
> > Thanks,
> >
> > Sascha
> >
> >
> > --
> > Mariette Yvinec
> > Geometrica project team
> > INRIA  Sophia-Antipolis

template<typename MeshDomain, typename OutputIterator>
void mitk::LabeledImageToTetrahedralGridFilter::ConstructSeedPoints(const mitk::Image* image,
                                                                    const MeshDomain* domain,
                                                                    const SeedPointDistanceToLabelsMap& lengths,
                                                                    OutputIterator pts)
{
  typedef typename MeshDomain::Point_3 PointType;

  mitk::Point3D origin = image->GetGeometry()->GetOrigin();
  mitk::Point3D endPoint = image->GetGeometry()->GetCornerPoint(7);

  for (SeedPointDistanceToLabelsMap::const_iterator iter = lengths.begin();
       iter != lengths.end(); ++iter)
  {
    const std::set<mitk::LabeledImageToTetrahedralGridFilter::LabelType>& labels = iter->second;

    int xCount = (endPoint[0] - origin[0]) / iter->first;
    std::vector<std::vector<std::pair<mitk::Point3D,mitk::LabeledImageToTetrahedralGridFilter::LabelType> > > seedPoints(omp_get_max_threads());
    for (std::size_t i = 0; i < seedPoints.size(); ++i)
    {
      seedPoints[i].reserve(1000);
    }

    #pragma omp parallel for
    for (int i = 0; i < xCount; ++i)
    {
      mitk::Point3D seedPointCandidate = origin;
      seedPointCandidate[0] = origin[0] + i * iter->first;
      while (seedPointCandidate[1] < endPoint[1])
      {
        seedPointCandidate[2] = origin[2];
        while (seedPointCandidate[2] < endPoint[2])
        {
          mitk::LabeledImageToTetrahedralGridFilter::LabelType label = const_cast<mitk::Image*>(image)->GetPixelValueByWorldCoordinate(seedPointCandidate);
          if (label != 0 && labels.find(label) != labels.end())
          {
            seedPoints[omp_get_thread_num()].push_back(std::make_pair(seedPointCandidate, label));
          }
          seedPointCandidate[2] += iter->first;
        }
        seedPointCandidate[1] += iter->first;
      }
    }

    for (std::size_t i = 0; i < seedPoints.size(); ++i)
    {
      for (std::size_t j = 0; j < seedPoints[i].size(); ++j)
      {
        const mitk::Point3D& seedPoint = seedPoints[i][j].first;
        *pts++ = std::make_pair(PointType(seedPoint[0], seedPoint[1], seedPoint[2]), domain->index_from_subdomain_index(seedPoints[i][j].second));
      }
    }
  }
}

vtkUnstructuredGrid* mitk::LabeledImageToTetrahedralGridFilter::GenerateUnstructuredGridFromImage(const mitk::Image* mitkImage)
{
  typedef mitk::CGALLabeledImage_3 CGALImage;
  CGALImage image(mitkImage);

  typedef CGAL::Labeled_image_mesh_domain_3<CGALImage, CGAL::Exact_predicates_inexact_constructions_kernel> MeshDomain;
  MeshDomain domain(image);

  typedef CGAL::Mesh_triangulation_3<MeshDomain>::type Tr;
  typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
  typedef CGAL::Mesh_criteria_3<Tr> MeshCriteria;
  typedef CGAL::Mesh_constant_domain_field_3<MeshDomain::R, MeshDomain::Index> SizingField;
  typedef MeshCriteria::Facet_criteria FacetCriteria;
  typedef MeshCriteria::Cell_criteria CellCriteria;

  SizingField cellRadii(m_CellSize);
  for (std::map<LabelType,float>::iterator i = m_LabelToCellSize.begin();
       i != m_LabelToCellSize.end(); ++i)
  {
    cellRadii.set_size(i->second, 3, domain.index_from_subdomain_index(i->first));
  }

  SizingField facetRadii(m_FacetSize);
  for (std::map<LabelType,float>::iterator i = m_LabelToFacetSize.begin();
       i != m_LabelToFacetSize.end(); ++i)
  {
    facetRadii.set_size(i->second, 3, domain.index_from_subdomain_index(i->first));
  }

  FacetCriteria facet_crit(
        m_FacetAngle, // angle
        facetRadii, // size
        m_FacetDistance); // approximation
  CellCriteria cell_crit(
        m_CellRadiusEdgeRatio, // radius-edge ratio
        cellRadii); // size
  MeshCriteria mesh_crit(facet_crit, cell_crit);

  // Initialize c3t3
  C3t3 c3t3;

  typedef MeshDomain::Point_3 Point_3;
  typedef MeshDomain::Index Index;
  typedef std::vector<std::pair<Point_3, Index> > Initial_points_vector;
  typedef typename Initial_points_vector::iterator Ipv_iterator;
  typedef typename C3t3::Vertex_handle Vertex_handle;

  // Mesh initialization : get some points and add them to the mesh
  Initial_points_vector initial_points;
  int pointDimension = 3;
  this->ConstructSeedPoints(mitkImage, &domain, m_SeedPointDistanceToLabels, std::back_inserter(initial_points));

  if (initial_points.empty())
  {
    domain.construct_initial_points_object()(std::back_inserter(initial_points));
    pointDimension = 2;
  }

  // Insert points and set their index and dimension
  for ( Ipv_iterator it = initial_points.begin() ;
       it != initial_points.end() ;
       ++it )
  {
    Vertex_handle v = c3t3.triangulation().insert(it->first);

    // v could be null if point is hidden
    if ( v != Vertex_handle() )
    {
      c3t3.set_dimension(v, pointDimension);
      c3t3.set_index(v, it->second);
    }
  }

  CGAL_assertion( c3t3.triangulation().dimension() == 3 );

  // Build mesher and launch refinement process
  // Don't reset c3t3 as we just created it
  CGAL::refine_mesh_3(c3t3, domain, mesh_crit, CGAL::parameters::no_reset_c3t3());

  vtkUnstructuredGrid* vtkUG = mitk::CGALMeshUtil::CGALMeshToVtk(c3t3);
  return vtkUG;
}