CGAL users discussion list

[Michael Bieri <>]

Hi everybody,

I'm currently using CGAL to mesh a segmented image. I define the criteria as follows:

Image_sizing_field cell_size_fct(elem_size_image);

Mesh_criteria criteria(facet_angle=30, facet_size=cell_size_fct, facet_distance=cell_size_fct, cell_radius_edge_ratio=3, cell_size=cell_size_fct);

I'm not quite sure whether it's ok to use the cell_size_fct also for facet_size and facet_distance, but so far, it works for me. I'm open for advices on that.

Now, the Problem: When I add 1D features basically with "CGAL::polylines_to_protect" to protect the bounding box, I get very nasty results. The mesh along the boundary gets very coarse, it seems like no other points were allowed close to the boundary. I can fix the problem by specifying edge_size to a lower, constant value. BUT as soon as I try to use the Image_sizing_field cell_size_fct as edge_size, I get a compile error: 

In file included from /usr/local/include/CGAL/Mesh_criteria_3.h:32:
/usr/local/include/CGAL/Mesh_edge_criteria_3.h:63:38: error: no type named
'Point_3' in 'Image_sizing_field'
typename Sizing_field::Point_3,

Is it not possible to use a Image_sizing_field as edge_size criteria? Or how could I do this?

(You find the whole code attached.)

Best regards,

Michael

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/Labeled_image_mesh_domain_3.h>
#include <CGAL/Mesh_domain_with_polyline_features_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/Image_3.h>

#include <CGAL/Mesh_3/polylines_to_protect.h> // undocumented header

// Domain
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Labeled_image_mesh_domain_3<CGAL::Image_3,K> Image_domain;
typedef CGAL::Mesh_domain_with_polyline_features_3<Image_domain> Mesh_domain;
typedef K::FT FT;
typedef K::Point_3 Point;

#ifdef CGAL_CONCURRENT_MESH_3
typedef CGAL::Parallel_tag Concurrency_tag;
#else
typedef CGAL::Sequential_tag Concurrency_tag;
#endif
// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain,CGAL::Default,Concurrency_tag>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
// Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;

// To avoid verbose function and named parameters call
using namespace CGAL::parameters;

// For the 1D-features
bool add_1D_features(const CGAL::Image_3& image,
                     Mesh_domain& domain)
{
  typedef K::Point_3 Point_3;
  typedef Mesh_domain::Image_word_type Word_type; // that is `unsigned char`

  std::vector<std::vector<Point_3> > features_inside;
  std::vector<std::vector<Point_3> > polylines_on_bbox;
  CGAL::polylines_to_protect<Point_3, Word_type>(image, polylines_on_bbox,
                                                 features_inside.begin(),
                                                 features_inside.end());

  domain.add_features(polylines_on_bbox.begin(), polylines_on_bbox.end());
  return true;
}


// Sizing field from image
// VERY QUICK AND DIRTY
struct Image_sizing_field
{
  typedef ::FT FT;
  typedef Mesh_domain::Index Index;
 
  const CGAL::Image_3& img;

  std::size_t xdim;
  std::size_t ydim;
  std::size_t zdim;
  
  std::size_t xydim;
  std::size_t size;

  double vx;
  double vy;
  double vz;
  
  const double* data;
  
  Image_sizing_field(const CGAL::Image_3& img) : 
  	img(img),
	xdim(img.xdim()),
	ydim(img.ydim()),
	zdim(img.zdim()),
	size(img.size()),
	vx(img.vx()),
	vy(img.vy()),
	vz(img.vz()),
	data((double*)img.data())
  { 
    xydim = xdim * ydim;
  }
  
  FT operator()(const Point& p, const int, const Index&) const
  {
	FT x = p.x();
	FT y = p.y();
	FT z = p.z();

    //double v = img.trilinear_interpolation<double>(x, y, z, 0);
	//std::cout << x << " " << y << " " << z << std::endl;

	int xi = (int)(x / vx);
	int yi = (int)(y / vy);
	int zi = (int)(z / vz);
	
	xi = (xi < 0) ? 0 : xi;
	yi = (yi < 0) ? 0 : yi;
	zi = (zi < 0) ? 0 : zi;
	
	xi = (xi >= xdim) ? xdim - 1 : xi;
	yi = (yi >= ydim) ? ydim - 1 : yi;
	zi = (zi >= zdim) ? zdim - 1 : zi;  
	
	int idx = xydim * zi + xdim * yi + xi;
	if((idx < 0) || (idx >= size))
		std::cout << "Index for image access seems to be out of range.";
	double v = data[idx];
	
    return v;
  }
};

template <class T>
void print_image(const CGAL::Image_3& i){
	using namespace std;
	cout << "Image Info: " << endl;
	cout << "Size int: " << i.xdim() << " " << i.ydim() << " " << i.zdim() << endl;
	cout << "Voxel size: " << i.vx() << " " << i.vy() << " " << i.vz() << endl;
	size_t n = i.xdim() * i.ydim() * i.zdim();
	cout << "N voxels: " << n << endl;
	const T* data = (const T*)i.data();
	for(int j = 0; j < n; j ++){
		//cout << data[j] << " ";
	}
}

int main(int argc, char* argv[])
{
  std::cout << "CGAL meshing with fields for region and elem size (version 0.1)" << std::endl;
  const char* fname_out = argv[1];
  const char* fname_regions = argv[2];
  const char* fname_elem_size = argv[3];
  // Load region image
  CGAL::Image_3 region_image;
  if(!region_image.read(fname_regions)){
    std::cerr << "Error: Cannot read region file " <<  fname_regions << std::endl;
    return EXIT_FAILURE;
  }
  // Load sizing image
  CGAL::Image_3 elem_size_image;
  if(!elem_size_image.read(fname_elem_size)){
    std::cerr << "Error: Cannot read elem size file " <<  fname_elem_size << std::endl;
    return EXIT_FAILURE;
  }
  // Print sizing file for test
  print_image<double>(elem_size_image);
  
  // Domain
  Mesh_domain domain(region_image);
  
  if(!add_1D_features(region_image, domain)) {
    return EXIT_FAILURE;
  }
  
  // Mesh criteria
  Image_sizing_field cell_size_fct(elem_size_image);
  Mesh_criteria criteria(edge_size=cell_size_fct, facet_angle=30, facet_size=cell_size_fct, facet_distance=cell_size_fct,
                         cell_radius_edge_ratio=3, cell_size=cell_size_fct);
  C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
  // Output
  std::ofstream medit_file(fname_out);
  c3t3.output_to_medit(medit_file);
  std::cout << std::endl << "*** Meshing done ***" << std::endl;
  return 0;
}