diff --git a/cgalv2m.m b/cgalv2m.m
index 3a37499..428571e 100644
--- a/cgalv2m.m
+++ b/cgalv2m.m
@@ -10,7 +10,7 @@
 % author: Qianqian Fang (q.fang at neu.edu)
 %
 % input:
-%	 vol: a volumetric binary image
+%	 vol: a volumetric binary or labeled image
 %	 ix,iy,iz: subvolume selection indices in x,y,z directions
 %	 opt: parameters for CGAL mesher, if opt is a structure, then
 %	     opt.radbound: defines the maximum surface element size
@@ -37,10 +37,10 @@
 
 fprintf(1,'creating surface and tetrahedral mesh from a multi-domain volume ...\n');
 
-dtype=class(vol);
-if(~(islogical(vol) || strcmp(dtype,'uint8')))
-	error('cgalmesher can only handle uint8 volumes, you have to convert your image to unit8 first.');
-end
+%dtype=class(vol);
+%if(~(islogical(vol) || strcmp(dtype,'uint8')))
+%	error('cgalmesher can only handle uint8 volumes, you have to convert your image to unit8 first.');
+%end
 
 if(~any(vol))
         error('no labeled regions found in the input volume.');
diff --git a/tools/cgalmesh/initialise_triangulation_from_labeled_image.h b/tools/cgalmesh/initialise_triangulation_from_labeled_image.h
new file mode 100644
index 0000000..2d9b65b
--- /dev/null
+++ b/tools/cgalmesh/initialise_triangulation_from_labeled_image.h
@@ -0,0 +1,219 @@
+// Copyright (c) 2015,2016 GeometryFactory
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL: https://github.com/CGAL/cgal/blob/releases/CGAL-5.0-I-190/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h $
+// $Id: initialize_triangulation_from_labeled_image.h 254d60f 2019-10-19T15:23:19+02:00 Sébastien Loriot
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+//
+// Author(s)     : Laurent Rineau
+#ifndef CGAL_MESH_3_INITIALIZE_TRIANGULATION_FROM_LABELED_IMAGE_H
+#define CGAL_MESH_3_INITIALIZE_TRIANGULATION_FROM_LABELED_IMAGE_H
+#include <CGAL/license/Mesh_3.h>
+#include <CGAL/Mesh_3/search_for_connected_components_in_labeled_image.h>
+#include <CGAL/Mesh_3/squared_distance_Point_3_Triangle_3.h>
+#include <CGAL/make_mesh_3.h>
+#include <CGAL/enum.h>
+#include <CGAL/iterator.h>
+#include <CGAL/point_generators_3.h>
+#include <CGAL/Image_3.h>
+#include <iostream>
+#include <queue>
+template <typename Point>
+struct Get_point
+{
+	const double vx, vy, vz;
+	Get_point(const CGAL::Image_3* image)
+		: vx(image->vx())
+		, vy(image->vy())
+		, vz(image->vz())
+	{}
+	Point operator()(const std::size_t i,
+		const std::size_t j,
+		const std::size_t k) const
+	{
+		return Point(double(i) * vx, double(j) * vy, double(k) * vz);
+	}
+};
+template<class C3T3, class MeshDomain, class MeshCriteria>
+void init_tr_from_labeled_image_call_init_features(C3T3&,
+	const MeshDomain&,
+	const MeshCriteria&,
+	CGAL::Tag_false)
+{
+}
+template<class C3T3, class MeshDomain, class MeshCriteria>
+void init_tr_from_labeled_image_call_init_features(C3T3& c3t3,
+	const MeshDomain& domain,
+	const MeshCriteria& criteria,
+	CGAL::Tag_true)
+{
+	CGAL::Mesh_3::internal::init_c3t3_with_features(c3t3,
+		domain,
+		criteria);
+	std::cout << c3t3.triangulation().number_of_vertices()
+		<< " initial points on 1D-features" << std::endl;
+}
+template<class C3T3, class MeshDomain, class MeshCriteria,
+	typename Image_word_type>
+	void initialize_triangulation_from_labeled_image(C3T3& c3t3,
+		const MeshDomain&   domain,
+		const CGAL::Image_3& image,
+		const MeshCriteria& criteria,
+		Image_word_type,
+		bool protect_features = false
+	)
+{
+	typedef typename C3T3::Triangulation       Tr;
+	typedef typename Tr::Geom_traits           Gt;
+	typedef typename Gt::FT                    FT;
+	typedef typename Tr::Weighted_point        Weighted_point;
+	typedef typename Tr::Bare_point            Bare_point;
+	typedef typename Tr::Segment               Segment_3;
+	typedef typename Tr::Vertex_handle         Vertex_handle;
+	typedef typename Tr::Cell_handle           Cell_handle;
+	typedef typename Gt::Vector_3              Vector_3;
+	typedef MeshDomain                         Mesh_domain;
+	Tr& tr = c3t3.triangulation();
+	typename Gt::Compare_weighted_squared_radius_3 cwsr =
+		tr.geom_traits().compare_weighted_squared_radius_3_object();
+	typename Gt::Construct_point_3 cp =
+		tr.geom_traits().construct_point_3_object();
+	typename Gt::Construct_weighted_point_3 cwp =
+		tr.geom_traits().construct_weighted_point_3_object();
+	if (protect_features) {
+		init_tr_from_labeled_image_call_init_features
+		(c3t3, domain, criteria,
+			CGAL::Mesh_3::internal::Has_features<Mesh_domain>());
+	}
+	const double max_v = (std::max)((std::max)(image.vx(),
+		image.vy()),
+		image.vz());
+	typedef std::vector<std::pair<Bare_point, std::size_t> > Seeds;
+	Seeds seeds;
+	Get_point<Bare_point> get_point(&image);
+	std::cout << "Searching for connected components..." << std::endl;
+	CGAL::Identity<Image_word_type> no_transformation;
+	search_for_connected_components_in_labeled_image(image,
+		std::back_inserter(seeds),
+		CGAL::Emptyset_iterator(),
+		no_transformation,
+		get_point,
+		Image_word_type());
+	std::cout << "  " << seeds.size() << " components were found." << std::endl;
+	std::cout << "Construct initial points..." << std::endl;
+	for (typename Seeds::const_iterator it = seeds.begin(), end = seeds.end();
+		it != end; ++it)
+	{
+		const double radius = double(it->second + 1)* max_v;
+		CGAL::Random_points_on_sphere_3<Bare_point> points_on_sphere_3(radius);
+		typename Mesh_domain::Construct_intersection construct_intersection =
+			domain.construct_intersection_object();
+		std::vector<Vector_3> directions;
+		if (it->second < 2) {
+			// shoot in six directions
+			directions.push_back(Vector_3(-radius, 0, 0));
+			directions.push_back(Vector_3(+radius, 0, 0));
+			directions.push_back(Vector_3(0, -radius, 0));
+			directions.push_back(Vector_3(0, +radius, 0));
+			directions.push_back(Vector_3(0, 0, -radius));
+			directions.push_back(Vector_3(0, 0, +radius));
+		}
+		else {
+			for (int i = 0; i < 20; ++i)
+			{
+				// shoot 20 random directions
+				directions.push_back(*points_on_sphere_3++ - CGAL::ORIGIN);
+			}
+		}
+		for (const Vector_3& v : directions)
+		{
+			const Bare_point test = it->first + v;
+			const typename Mesh_domain::Intersection intersect =
+				construct_intersection(Segment_3(it->first, test));
+			if (std::get<2>(intersect) != 0)
+			{
+				const Bare_point& bpi = std::get<0>(intersect);
+				Weighted_point pi = cwp(bpi);
+				// This would cause trouble to optimizers
+				// check pi will not be hidden
+				typename Tr::Locate_type lt;
+				int li, lj;
+				Cell_handle pi_cell = tr.locate(pi, lt, li, lj);
+				if (lt != Tr::OUTSIDE_AFFINE_HULL) {
+					switch (tr.dimension())
+					{ //skip dimension 0
+					case 1:
+						if (tr.side_of_power_segment(pi_cell, pi, true) != CGAL::ON_BOUNDED_SIDE)
+							continue;
+						break;
+					case 2:
+						if (tr.side_of_power_circle(pi_cell, 3, pi, true) != CGAL::ON_BOUNDED_SIDE)
+							continue;
+						break;
+					case 3:
+						if (tr.side_of_power_sphere(pi_cell, pi, true) != CGAL::ON_BOUNDED_SIDE)
+							continue;
+					}
+				}
+				//check pi is not inside a protecting ball
+				std::vector<Vertex_handle> conflict_vertices;
+				if (tr.dimension() == 3)
+				{
+					tr.vertices_on_conflict_zone_boundary(pi, pi_cell
+						, std::back_inserter(conflict_vertices));
+				}
+				else
+				{
+					for (typename Tr::Finite_vertices_iterator vit = tr.finite_vertices_begin();
+						vit != tr.finite_vertices_end(); ++vit)
+					{
+						const Weighted_point& wp = tr.point(vit);
+						if (cwsr(wp, FT(0)) == CGAL::SMALLER) // 0 < wp's weight
+							conflict_vertices.push_back(vit);
+					}
+				}
+				bool pi_inside_protecting_sphere = false;
+				for (Vertex_handle cv : conflict_vertices)
+				{
+					if (tr.is_infinite(cv))
+						continue;
+					const Weighted_point& cv_wp = tr.point(cv);
+					if (cwsr(cv_wp, FT(0)) == CGAL::EQUAL) // 0 == wp's weight
+						continue;
+					// if the (squared) distance between bpi and cv is smaller or equal than cv's weight
+					if (cwsr(cv_wp, -tr.min_squared_distance(bpi, cp(cv_wp))) != CGAL::LARGER)
+					{
+						pi_inside_protecting_sphere = true;
+						break;
+					}
+				}
+				if (pi_inside_protecting_sphere)
+					continue;
+				const typename Mesh_domain::Index index = std::get<1>(intersect);
+				Vertex_handle v = tr.insert(pi);
+				// `v` could be null if `pi` is hidden by other vertices of `tr`.
+				CGAL_assertion(v != Vertex_handle());
+				c3t3.set_dimension(v, 2); // by construction, points are on surface
+				c3t3.set_index(v, index);
+			}
+			// else
+			// {
+			//   std::cerr <<
+			//     boost::format("Error. Segment (%1%, %2%) does not intersect the surface!\n")
+			//     % it->first % test;
+			// }
+		}
+	}
+	std::cout << "  " << tr.number_of_vertices() << " initial points." << std::endl;
+	if (c3t3.triangulation().dimension() != 3)
+	{
+		std::cout << "  not enough points: triangulation.dimension() == "
+			<< c3t3.triangulation().dimension() << std::endl;
+		CGAL::Mesh_3::internal::init_c3t3(c3t3, domain, criteria, 20);
+		std::cout << "  -> " << tr.number_of_vertices() << " initial points." << std::endl;
+	}
+}
+#endif // CGAL_MESH_3_INITIALIZE_TRIANGULATION_FROM_LABELED_IMAGE_H
\ No newline at end of file
diff --git a/tools/cgalmesh/mesh_3D_image.cpp b/tools/cgalmesh/mesh_3D_image.cpp
index d2e7349..9a8e6e4 100644
--- a/tools/cgalmesh/mesh_3D_image.cpp
+++ b/tools/cgalmesh/mesh_3D_image.cpp
@@ -8,141 +8,137 @@
 #include <CGAL/make_mesh_3.h>
 #include <CGAL/Image_3.h>
 
-// start PV //
 #include <string>
 #include <vector>
 #include <CGAL/Mesh_constant_domain_field_3.h>
-// end PV //
 
-// Domain
-struct K: public CGAL::Exact_predicates_inexact_constructions_kernel {};
+#include "initialise_triangulation_from_labeled_image.h"
+
+// Image
 typedef CGAL::Image_3 Image;
-typedef CGAL::Labeled_image_mesh_domain_3<Image,K> Mesh_domain;
+
+// Domain
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef CGAL::Labeled_mesh_domain_3<K> Mesh_domain;
 
 // Triangulation
-typedef CGAL::Mesh_triangulation_3<Mesh_domain>::type Tr;
+typedef CGAL::Parallel_tag Concurrency_tag;
+typedef CGAL::Mesh_triangulation_3<Mesh_domain, CGAL::Default, Concurrency_tag>::type Tr;
 typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
 
 // Mesh Criteria
 typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
-typedef Mesh_criteria::Facet_criteria    Facet_criteria;
-typedef Mesh_criteria::Cell_criteria     Cell_criteria;
+typedef Mesh_criteria::Facet_criteria Facet_criteria;
+typedef Mesh_criteria::Cell_criteria Cell_criteria;
 
-// start PV //
-typedef CGAL::Mesh_constant_domain_field_3<Mesh_domain::R,
-                                           Mesh_domain::Index> Sizing_field_cell;
-// end PV //
+typedef CGAL::Mesh_constant_domain_field_3<Mesh_domain::R, Mesh_domain::Index> Sizing_field_cell;
 
 
-void usage(char *exename){
-	printf("usage:\n\t%s input.inr output.mesh <angle|30> <surf-size|6> <approx|4> <rad-edge-ratio|3> \
+void usage(char *exename) {
+	printf(
+"usage:\n\t%s input.inr output.mesh <angle|30> <surf-size|6> <approx|4> <rad-edge-ratio|3> \
 <tetra-size|8> <randomseed|-1>\n\
-example:\n\t%s input.inr output.mesh 30 6 4 3 8 123456789\n",exename,exename);
+example:\n\t%s input.inr output.mesh 30 6 4 3 8 123456789\n", exename, exename);
 	exit(1);
 }
-int main(int argc,char *argv[])
+
+int main(int argc, char *argv[])
 {
-  // Loads image
-
-// start PV //
-  //float angle=30.f,ssize=6.f,approx=4.f,reratio=3.f,vsize=8.f;
-  float angle=30.f,ssize=6.f,approx=4.f,reratio=3.f,maxvol=0.f;
-// end PV //
-  int labelid=0, lid;
-
-  printf("Volume/Surface Mesh Generation Utility (Based on CGAL %s)\n\
-(modified for iso2mesh by Qianqian Fang and Peter Varga)\nhttp://iso2mesh.sf.net\n\n",CGAL_VERSION_STR);
-
-  if(argc!=3&&argc!=8&&argc!=9){
-  	usage(argv[0]);
-  }
-  if(argc>=8){
-  	angle=atof(argv[3]);
-  	ssize=atof(argv[4]);
-  	approx=atof(argv[5]);
-  	reratio=atof(argv[6]);
-// start PV //
-  	//vsize=atof(argv[7]);
-// end PV //
-  }
-  if(argc==9 && atoi(argv[8])>0){
-	printf("RNG seed=%d\n",atoi(argv[8]));
-        CGAL::Random rd(atoi(argv[8]));
-        CGAL::Random::State st;
-        rd.save_state(st);
-        CGAL::default_random.restore_state(st);
-  }
-  Image image;
-  image.read(argv[1]);
-  Mesh_domain domain(image);
-
-  // Mesh criteria
-  Facet_criteria facet_criteria(angle, ssize, approx); // angle, size, approximation
-
-// start PV //
-  //Cell_criteria cell_criteria(reratio, vsize); // radius-edge ratio, size
-
-  std::string vsize;
-  vsize=argv[7];
-
-  std::vector<double> vsize_vect;
-  std::vector<int> labels_vect;
-
-  std::stringstream stream_vsize(vsize);
-  std::string word_vsize;
-
-  while( std::getline(stream_vsize, word_vsize, ':') ){
-    int len=sscanf(word_vsize.c_str(), "%d=%f", &lid, &maxvol);
-    if(len==2){
-      labelid=lid;
-      if(maxvol<=0.f){
-         std::cerr << "cell volume must be positive" << std::endl;
-         exit(-2);
-      }
-      labels_vect.push_back(labelid);
-      vsize_vect.push_back(maxvol);
-    }else{
-      len=sscanf(word_vsize.c_str(), "%f", &maxvol);
-      if(len!=1 || maxvol<0 || word_vsize.find_first_of('=')!=std::string::npos){
-         std::cerr << "invalid sizing field label '" << word_vsize << "', please check your command" << std::endl;
-         exit(-1);
-      }
-      labels_vect.push_back(++labelid);
-      vsize_vect.push_back(maxvol);
-    }
-  }
-
-  float max_vsize = *max_element(vsize_vect.begin(),vsize_vect.end());
-
-  Sizing_field_cell vsize_cell(max_vsize);
-  int volume_dimension = 3;
-
-  for (int i=0; i<vsize_vect.size(); i++)
-    vsize_cell.set_size(vsize_vect[i], volume_dimension, 
-                domain.index_from_subdomain_index(labels_vect[i]));
-
-  std::cout << "Mesh sizes are (label=size) ";
-  for (int i=0; i<vsize_vect.size(); i++){
-    std::cout << "(" << labels_vect[i] <<  "=" << vsize_vect[i] << ") ";
-  }
-  std::cout << std::endl;
-
-  Cell_criteria cell_criteria(reratio, vsize_cell);
-// end PV //
-
-  Mesh_criteria criteria(facet_criteria, cell_criteria);
-
-  // Meshing
-  C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
-
-  // Output
-  std::ofstream medit_file;
-  if(argc>=8)
-	medit_file.open(argv[2]);
-  else
-	medit_file.open("output.mesh");
-  c3t3.output_to_medit(medit_file);
-  medit_file.close();
-
-  return 0;
+	//float angle=30.f,ssize=6.f,approx=4.f,reratio=3.f,vsize=8.f;
+	float angle = 30.f, ssize = 6.f, approx = 4.f, reratio = 3.f, maxvol = 0.f;
+	int labelid = 0, lid;
+
+	printf("Volume/Surface Mesh Generation Utility (Based on CGAL %s)\n\
+(modified for iso2mesh by Qianqian Fang and Peter Varga)\nhttp://iso2mesh.sf.net\n\n", CGAL_VERSION_STR);
+
+	if (argc != 3 && argc != 8 && argc != 9) {
+		usage(argv[0]);
+	}
+	if (argc >= 8) {
+		angle = (float)atof(argv[3]);
+		ssize = (float)atof(argv[4]);
+		approx = (float)atof(argv[5]);
+		reratio = (float)atof(argv[6]);
+	}
+	if (argc == 9 && atoi(argv[8]) > 0) {
+		printf("RNG seed=%d\n", atoi(argv[8]));
+		CGAL::Random rd(atoi(argv[8]));
+		CGAL::Random::State st;
+		rd.save_state(st);
+		CGAL::default_random.restore_state(st);
+	}
+	
+	Image image;
+	image.read(argv[1]);
+	Mesh_domain domain = Mesh_domain::create_labeled_image_mesh_domain(image);
+
+	// Mesh criteria
+	Facet_criteria facet_criteria(angle, ssize, approx); // angle, size, approximation
+
+	std::string vsize;
+	vsize = argv[7];
+
+	std::vector<double> vsize_vect;
+	std::vector<int> labels_vect;
+
+	std::stringstream stream_vsize(vsize);
+	std::string word_vsize;
+
+	while (std::getline(stream_vsize, word_vsize, ':')) {
+		int len = sscanf(word_vsize.c_str(), "%d=%f", &lid, &maxvol);
+		if (len == 2) {
+			labelid = lid;
+			if (maxvol <= 0.f) {
+				std::cerr << "cell volume must be positive" << std::endl;
+				exit(-2);
+			}
+			labels_vect.push_back(labelid);
+			vsize_vect.push_back(maxvol);
+		}
+		else {
+			len = sscanf(word_vsize.c_str(), "%f", &maxvol);
+			if (len != 1 || maxvol < 0 || word_vsize.find_first_of('=') != std::string::npos) {
+				std::cerr << "invalid sizing field label '" << word_vsize << "', please check your command" << std::endl;
+				exit(-1);
+			}
+			labels_vect.push_back(++labelid);
+			vsize_vect.push_back(maxvol);
+		}
+	}
+
+	float max_vsize = (float)*max_element(vsize_vect.begin(), vsize_vect.end());
+
+	Sizing_field_cell vsize_cell(max_vsize);
+	int volume_dimension = 3;
+
+	for (int i = 0; i < vsize_vect.size(); i++)
+		vsize_cell.set_size(vsize_vect[i], volume_dimension,
+			domain.index_from_subdomain_index(labels_vect[i]));
+
+	std::cout << "Mesh sizes are (label=size) ";
+	for (int i = 0; i < vsize_vect.size(); i++) {
+		std::cout << "(" << labels_vect[i] << "=" << vsize_vect[i] << ") ";
+	}
+	std::cout << std::endl;
+
+	Cell_criteria cell_criteria(reratio, vsize_cell);
+
+	Mesh_criteria criteria(facet_criteria, cell_criteria);
+
+	// Meshing
+	C3t3 c3t3;
+	CGAL_IMAGE_IO_CASE(image.image(), initialize_triangulation_from_labeled_image(c3t3, domain, image, criteria, (Word)0))
+
+	CGAL::refine_mesh_3(c3t3, domain, criteria);
+
+	// Output
+	std::ofstream medit_file;
+	if (argc >= 8)
+		medit_file.open(argv[2]);
+	else
+		medit_file.open("output.mesh");
+	c3t3.output_to_medit(medit_file);
+	medit_file.close();
+
+	return 0;
 }