songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

WIP: implementation for scalar image 

... tested only on segmented images!
This commit is contained in:
Laurent Rineau 2018-01-15 16:14:26 +01:00
parent 306d3031ab
commit 1d5e9e9150
1 changed files with 229 additions and 83 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

312
Mesh_3/include/CGAL/Mesh_3/polylines_to_protect.h
View File

@ -37,6 +37,7 @@
// CGAL::Null_subdomain_index
#include <boost/utility.hpp> // for boost::prior
#include <boost/foreach.hpp>
#include <boost/optional.hpp>
#include <CGAL/Search_traits_3.h>
#include <CGAL/Orthogonal_incremental_neighbor_search.h>
@ -58,6 +59,116 @@ struct Returns_midpoint {
= K().construct_midpoint_3_object();
return a < b ? midpt(a, b) : midpt(b, a);
}
double operator()(const NT,
const NT) const
{
return 0.5;
}
};
class Isoline_equation {
double a;
double b;
double c;
double d;
public:
Isoline_equation(double a, double b, double c, double d)
: a(a)
, b(b)
, c(c)
, d(d)
{}
double y(double x) const {
return ( x*(a-b)-a ) /
( c-a+x*(a-b-c+d) );
// If (a+d) == (b+c), then the curve is a straight line.
}
double x(double y) const {
return ( y*(a-c)-a ) /
( b-a+y*(a-b-c+d) );
// If (a+d) == (b+c), then the curve is a straight line.
}
}; // end of class Isoline_equation
template <typename G_manip, typename vertex_descriptor, typename Geom_traits>
class Insert_curve_in_graph {
typedef typename Geom_traits::Point_3 Point;
typedef typename Geom_traits::Vector_3 Vector;
G_manip& g_manip;
const Geom_traits& gt;
const typename Geom_traits::Construct_translated_point_3 translate;
const typename Geom_traits::Construct_scaled_vector_3 scale;
const int prec;
public:
Insert_curve_in_graph(G_manip& g_manip,
const Geom_traits& gt,
const int prec = 10)
: g_manip(g_manip)
, gt(gt)
, translate(gt.construct_translated_point_3_object())
, scale(gt.construct_scaled_vector_3_object())
, prec(prec)
{}
struct Coords {
double x;
double y;
};
void insert_curve(Isoline_equation equation,
vertex_descriptor start_v,
vertex_descriptor end_v,
Coords start,
Coords end,
Point p00,
Vector vx,
Vector vy)
{
if(CGAL::abs(start.x - end.x) >= CGAL::abs(start.y - end.y)) {
insert_curve(equation, &Isoline_equation::y,
start_v, end_v,
start.x, end.x, p00, vx, vy);
} else {
insert_curve(equation, &Isoline_equation::x,
start_v, end_v,
start.y, end.y, p00, vy, vx);
}
}
private:
void insert_curve(Isoline_equation equation,
double (Isoline_equation::* f)(double) const,
vertex_descriptor start_v,
vertex_descriptor end_v,
double start,
double end,
Point p00,
Vector vx,
Vector vy)
{
const double step = (end - start)/prec;
const double stop = end-step/2;
const bool step_is_positive = (step > 0);
vertex_descriptor old = start_v;
vertex_descriptor v_int;
for(double x = start + step;
(step_is_positive ? x < stop : x > stop);
x += step)
{
const double y = (equation.*f)(x);
const Point inter_p =
translate(translate(p00,
scale(vx, x)),
scale(vy, y));
v_int = g_manip.get_vertex(inter_p);
g_manip.try_add_edge(old, v_int);
old = v_int;
}
g_manip.try_add_edge(v_int, end_v);
}
};
template <typename Pixel,
@ -225,15 +336,18 @@ template <typename P,
typename InterpolationFunctor,
typename PolylineInputIterator>
void
polylines_to_protect(const CGAL::Image_3& cgal_image,
const double vx, const double vy, const double vz,
std::vector<std::vector<P> >& polylines,
Image_word_type*,
Null_subdomain_index null,
DomainFunctor domain_fct,
InterpolationFunctor interpolate,
PolylineInputIterator existing_polylines_begin,
PolylineInputIterator existing_polylines_end)
polylines_to_protect
(const CGAL::Image_3& cgal_image,
const double vx, const double vy, const double vz,
std::vector<std::vector<P> >& polylines,
Image_word_type*,
Null_subdomain_index null,
DomainFunctor domain_fct,
InterpolationFunctor interpolate,
PolylineInputIterator existing_polylines_begin,
PolylineInputIterator existing_polylines_end,
boost::optional<Image_word_type> scalar_interpolation_value = boost::none,
int prec = 10)
{
typedef typename DomainFunctor::result_type Domain_type;
typedef typename Kernel_traits<P>::Kernel K;
@ -249,10 +363,18 @@ polylines_to_protect(const CGAL::Image_3& cgal_image,
const int image_dims[3] = { xdim, ydim, zdim };
Graph graph;
internal::Mesh_3::Graph_manipulations<Graph,
Point_3,
Image_word_type,
InterpolationFunctor> g_manip(graph, interpolate);
using namespace CGAL::internal::Mesh_3;
typedef internal::Mesh_3::Graph_manipulations<Graph,
Point_3,
Image_word_type,
InterpolationFunctor> G_manip;
G_manip g_manip(graph, interpolate);
Insert_curve_in_graph<G_manip,
vertex_descriptor,
K> insert_curve_in_graph(g_manip, K(), prec);
std::size_t
case4 = 0, // 4 colors
@ -306,9 +428,8 @@ polylines_to_protect(const CGAL::Image_3& cgal_image,
{ pix11, Point_3(), Domain_type(), 0 } }} }};
std::map<Image_word_type, int> pixel_values_set;
for(int ii = 0; ii < 2; ++ii)
for(int jj = 0; jj < 2; ++jj)
{
for(int ii = 0; ii < 2; ++ii) {
for(int jj = 0; jj < 2; ++jj) {
const Pixel& pixel = square[ii][jj].pixel;
double x = pixel[0] * vx;
double y = pixel[1] * vy;
@ -319,8 +440,15 @@ polylines_to_protect(const CGAL::Image_3& cgal_image,
pixel[1],
pixel[2]));
square[ii][jj].domain = domain_fct(square[ii][jj].word);
if(scalar_interpolation_value != boost::none) {
square[ii][jj].word =
Image_word_type(square[ii][jj].word -
(*scalar_interpolation_value));
}
++pixel_values_set[square[ii][jj].domain];
}
}
const Point_3& p00 = square[0][0].point;
const Point_3& p10 = square[1][0].point;
const Point_3& p01 = square[0][1].point;
@ -375,6 +503,7 @@ polylines_to_protect(const CGAL::Image_3& cgal_image,
CGAL_assertion(square[1][0].domain != square[1][1].domain);
CGAL_assertion(square[0][1].domain != square[1][1].domain);
CGAL_assertion(square[0][1].domain != square[1][0].domain);
case_4:
// case 4 or case 2-2
//
@ -393,7 +522,9 @@ case_4:
vertex_descriptor top = g_manip.split(p01, p11, v01, v11, out01, out11);
vertex_descriptor bottom = g_manip.split(p00, p10, v00, v10, out00, out10);
vertex_descriptor vmid = g_manip.get_vertex(midpoint(p00, p11));
vertex_descriptor vmid = g_manip.split(graph[left], graph[right],
v00, v10,
false, false);
g_manip.try_add_edge(left , vmid);
g_manip.try_add_edge(right , vmid);
g_manip.try_add_edge(top , vmid);
@ -432,37 +563,34 @@ case_4:
CGAL_assertion(square[1][1].domain != square[0][0].domain);
CGAL_assertion(square[0][1].domain != square[0][0].domain);
CGAL_assertion(square[0][1].domain != square[1][1].domain);
case_1_2_1:
++case121;
vertex_descriptor left = g_manip.split(p00, p01, v00, v01, out00, out01);
vertex_descriptor right = g_manip.split(p10, p11, v10, v11, out10, out11);
vertex_descriptor top = g_manip.split(p01, p11, v01, v11, out01, out11);
vertex_descriptor bottom = g_manip.split(p00, p10, v00, v10, out00, out10);
vertex_descriptor old_left = left;
vertex_descriptor old_right = right;
vertex_descriptor v_int_left, v_int_right;
// approximate the arcs by 10 segments
// -> 9 intermediate vertices
for(double x = 0.05; x < 0.5; x+= 0.05)
{
const Point_3 inter_left =
translate(p00
, x * vector(p00, p10) // x
+ ((1.-2.*x)/(2.-3.*x)) * vector(p00, p01)); // y
const Point_3 inter_right =
translate(p11
, x * vector(p11, p01) // x
+ ((1.-2.*x)/(2.-3.*x)) * vector(p11, p10)); // y
v_int_left = g_manip.get_vertex(inter_left);
v_int_right = g_manip.get_vertex(inter_right);
g_manip.try_add_edge(old_left, v_int_left);
g_manip.try_add_edge(old_right, v_int_right);
old_left = v_int_left;
old_right = v_int_right;
Isoline_equation equation =
(scalar_interpolation_value == boost::none) ?
Isoline_equation(1, -1, -1, 0) :
Isoline_equation(v00, v10, v01, v11);
insert_curve_in_graph.insert_curve(equation,
left, bottom,
{0., interpolate(v00, v01) },
{interpolate(v00, v10), 0. },
p00,
p10 - p00,
p01 - p00);
if(scalar_interpolation_value == boost::none) {
equation = Isoline_equation(0, -1, -1, 1);
}
g_manip.try_add_edge(v_int_left, bottom);
g_manip.try_add_edge(v_int_right, top);
insert_curve_in_graph.insert_curve(equation,
top, right,
{interpolate(v01, v11), 1. },
{1., interpolate(v10, v11) },
p00,
p10 - p00,
p01 - p00);
} else {
// case 2-1-1
if(square[0][0].domain == square[1][0].domain) {
@ -483,6 +611,8 @@ case_4:
CGAL_assertion(square[0][0].domain != square[1][1].domain);
CGAL_assertion(square[1][0].domain != square[1][1].domain);
++case211;
// Note: this case cannot occur with non-segmented scalar
// images, because it needs three domains.
//
// A-------------B
// | | |
@ -514,32 +644,20 @@ case_4:
vertex_descriptor top = g_manip.split(p01, p11, v01, v11, out01, out11);
vertex_descriptor bottom = g_manip.split(p00, p10, v00, v10, out00, out10);
vertex_descriptor old_top = top;
vertex_descriptor old_bottom = bottom;
vertex_descriptor v_int_top, v_int_bottom;
// approximate the arcs by 10 segments
// -> 9 intermediate vertices
for(double x = 0.51666; x < 0.66; x+= 0.016666)
{
const Point_3 inter_top =
translate(p00
, x * vector(p00, p10) // x
+ ((1./x) - 1.) * vector(p00, p01)); // y
const Point_3 inter_bottom =
translate(p00
, x * vector(p00, p10) // x
+ (2.-(1./x)) * vector(p00, p01)); // y
v_int_top = g_manip.get_vertex(inter_top);
v_int_bottom = g_manip.get_vertex(inter_bottom);
g_manip.try_add_edge(old_top, v_int_top);
g_manip.try_add_edge(old_bottom, v_int_bottom);
old_top = v_int_top;
old_bottom = v_int_bottom;
}
g_manip.try_add_edge(v_int_bottom, v_inter);
g_manip.try_add_edge(v_int_top, v_inter);
insert_curve_in_graph.insert_curve(Isoline_equation(1, -1, 1, 0),
bottom, v_inter,
{ 0.5 , 0. },
{ 2./3., 0.5 },
p00,
p10 - p00,
p01 - p00);
insert_curve_in_graph.insert_curve(Isoline_equation(1, 0, 1, -1),
top, v_inter,
{ 0.5 , 1. },
{ 2./3., 0.5 },
p00,
p10 - p00,
p01 - p00);
g_manip.try_add_edge(right, v_inter);
} // end case 2-1-1
} // end `case 3:`
@ -585,13 +703,43 @@ case_4:
CGAL_assertion(square[1][0].domain!=square[0][1].domain);
vertex_descriptor top = g_manip.split(p01, p11, v01, v11, out01, out11);
vertex_descriptor bottom = g_manip.split(p00, p10, v00, v10, out00, out10);
g_manip.try_add_edge(top, bottom);
if(scalar_interpolation_value == boost::none) {
g_manip.try_add_edge(top, bottom);
} else {
insert_curve_in_graph.insert_curve(Isoline_equation(v00, v10,
v01, v11),
top, bottom,
{interpolate(v01,v11), 1.},
{interpolate(v00,v10), 0.},
p00,
p10 - p00,
p01 - p00);
}
} else {
// Else diagonal case case 2-2
// Same as the case with 4 colors
CGAL_assertion(square[0][0].domain==square[1][1].domain);
CGAL_assertion(square[1][0].domain==square[0][1].domain);
CGAL_assertion(square[0][0].domain!=square[0][1].domain);
if(scalar_interpolation_value != boost::none) {
// Compute the squared distance between the two branches of
// the hyperbola.
const double discrimant = double(v00) * v11 - double(v01) * v10;
const double squared_distance =
8 * CGAL::abs(discrimant) / CGAL::square(double(v00) - v10 - v01 + v11);
if(CGAL::square(prec) * squared_distance > 1)
{
// In that case, the case 1-2-1 will be applied
if(discrimant > 0.) {
// Vertical swap
std::swap(square[0][1], square[0][0]); std::swap(out01, out00);
std::swap(square[1][1], square[1][0]); std::swap(out11, out10);
}
goto case_1_2_1;
}
} // scalar interpolation
goto case_4;
}
}
@ -639,22 +787,20 @@ case_4:
CGAL_assertion(square[1][1].domain == square[0][0].domain);
CGAL_assertion(square[0][1].domain == square[0][0].domain);
vertex_descriptor bottom = g_manip.split(p00, p10, v00, v10, out00, out10);
vertex_descriptor old = bottom;
vertex_descriptor v_int;
for(double x = 0.55; x < 1.; x+= 0.05)
{
const Point_3 inter =
translate(p00
, x * vector(p00, p10) // x
+ (1.-1./(2.*x)) * vector(p00, p01)); // y
v_int = g_manip.get_vertex(inter);
g_manip.try_add_edge(old, v_int);
old = v_int;
}
vertex_descriptor right = g_manip.split(p10, p11, v10, v11, out10, out11);
g_manip.try_add_edge(v_int, right);
Isoline_equation equation =
(scalar_interpolation_value == boost::none) ?
Isoline_equation(1, -1, 1, 1) :
Isoline_equation(v00, v10, v01, v11);
insert_curve_in_graph.insert_curve(equation,
bottom, right,
{interpolate(v00, v10), 0. },
{1., interpolate(v10, v11) },
p00,
p10 - p00,
p01 - p00);
}
}
break;