inclusion.hpp 9.69 KB
Newer Older
Pierre Kraemer's avatar
Pierre Kraemer committed
1 2 3
/*******************************************************************************
 * CGoGN: Combinatorial and Geometric modeling with Generic N-dimensional Maps  *
 * version 0.1                                                                  *
4
 * Copyright (C) 2009-2012, IGG Team, LSIIT, University of Strasbourg           *
Pierre Kraemer's avatar
Pierre Kraemer committed
5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
 *                                                                              *
 * This library is free software; you can redistribute it and/or modify it      *
 * under the terms of the GNU Lesser General Public License as published by the *
 * Free Software Foundation; either version 2.1 of the License, or (at your     *
 * option) any later version.                                                   *
 *                                                                              *
 * This library is distributed in the hope that it will be useful, but WITHOUT  *
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or        *
 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License  *
 * for more details.                                                            *
 *                                                                              *
 * You should have received a copy of the GNU Lesser General Public License     *
 * along with this library; if not, write to the Free Software Foundation,      *
 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301 USA.           *
 *                                                                              *
20
 * Web site: http://cgogn.unistra.fr/                                           *
Pierre Kraemer's avatar
Pierre Kraemer committed
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
 * Contact information: cgogn@unistra.fr                                        *
 *                                                                              *
 *******************************************************************************/

#include "Algo/Geometry/intersection.h"
#include "Algo/Geometry/orientation.h"
#include "Algo/Geometry/basic.h"
#include "Algo/Geometry/plane.h"

#include <limits>

namespace CGoGN
{

namespace Algo
{

38 39 40
namespace Surface
{

Pierre Kraemer's avatar
Pierre Kraemer committed
41 42 43 44
namespace Geometry
{

template <typename PFP>
Sylvain Thery's avatar
Sylvain Thery committed
45
bool isConvex(typename PFP::MAP& map, Vol v, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, bool CCW)
Pierre Kraemer's avatar
Pierre Kraemer committed
46 47 48
{
	//get all the dart of the volume
	std::vector<Dart> vStore;
Sylvain Thery's avatar
Sylvain Thery committed
49
	map.foreach_dart_of_orbit(v, [&] (Dart d) { vStore.push_back(d); });
Pierre Kraemer's avatar
Pierre Kraemer committed
50 51 52

	bool convex = true;

Sylvain Thery's avatar
Sylvain Thery committed
53
	DartMarkerStore<typename PFP::MAP> m(map);
Pierre Kraemer's avatar
Pierre Kraemer committed
54 55 56 57 58
	for (std::vector<Dart>::iterator it = vStore.begin() ; it != vStore.end() && convex ; ++it)
	{
		Dart e = *it;
		if (!m.isMarked(e))
		{
Sylvain Thery's avatar
Sylvain Thery committed
59
			m.template markOrbit<EDGE>(e) ;
60
			convex = isTetrahedronWellOriented<PFP>(map, e, position, CCW) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
61 62 63 64 65 66
		}
	}

	return convex;
}

Sylvain Thery's avatar
Sylvain Thery committed
67

Pierre Kraemer's avatar
Pierre Kraemer committed
68
template <typename PFP>
69
bool isPointInVolume(typename PFP::MAP& map, Vol v, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, const typename PFP::VEC3& point)
Pierre Kraemer's avatar
Pierre Kraemer committed
70 71
{
	typedef typename PFP::VEC3 VEC3;
Sylvain Thery's avatar
Sylvain Thery committed
72
	typedef typename PFP::REAL REAL;
Pierre Kraemer's avatar
Pierre Kraemer committed
73 74 75

	//number of intersection between a ray and the volume must be odd
	int countInter = 0;
76
	int countInter2 = 0;
Pierre Kraemer's avatar
Pierre Kraemer committed
77

78
	VEC3 dir(0.9f,1.1f,1.3f);
79

Pierre Kraemer's avatar
Pierre Kraemer committed
80
	std::vector<VEC3> interPrec;
Sylvain Thery's avatar
Sylvain Thery committed
81
	interPrec.reserve(16);
82 83

	std::vector<Face> visitedFaces;			// Faces that are traversed
84
	visitedFaces.reserve(64);
85 86 87
	Face f(v.dart);
	visitedFaces.push_back(f);				// Start with the first face of v

Pierre Kraemer's avatar
Pierre Kraemer committed
88
	DartMarkerStore<typename PFP::MAP> mark(map);
89 90 91
	mark.markOrbit(f) ;

	for(unsigned int iface = 0; iface != visitedFaces.size(); ++iface)
Pierre Kraemer's avatar
Pierre Kraemer committed
92
	{
93
		f = visitedFaces[iface];
94
		VEC3 inter;
95
		bool interRes = intersectionLineConvexFace<PFP>(map, f, position, point, dir, inter);
96
		if (interRes)
Pierre Kraemer's avatar
Pierre Kraemer committed
97
		{
Sylvain Thery's avatar
Sylvain Thery committed
98
			// check if already intersect on same point (a vertex certainly)
Pierre Kraemer's avatar
Pierre Kraemer committed
99
			bool alreadyfound = false;
Sylvain Thery's avatar
Sylvain Thery committed
100 101
			for(typename std::vector<VEC3>::iterator it = interPrec.begin(); !alreadyfound && it != interPrec.end(); ++it)
			{
102 103
				if (Geom::arePointsEquals(*it, inter))
					alreadyfound = true;
Sylvain Thery's avatar
Sylvain Thery committed
104 105 106
			}

			if (!alreadyfound)
Pierre Kraemer's avatar
Pierre Kraemer committed
107
			{
Sylvain Thery's avatar
Sylvain Thery committed
108
				REAL v = dir * (inter - point);
109
				if (v > 0)
110
					++countInter;
111
				if (v < 0)
112
					++countInter2;
Sylvain Thery's avatar
Sylvain Thery committed
113
				interPrec.push_back(inter);
Pierre Kraemer's avatar
Pierre Kraemer committed
114 115
			}
		}
116
		// add all face neighbours to the table
117
		foreach_adjacent2<EDGE>(map, f, [&] (Face ff)
118
		{
119
			if(!mark.isMarked(ff)) // not already marked
120
			{
121 122
				visitedFaces.push_back(ff) ;
				mark.markOrbit(ff) ;
123
			}
124
		});
Pierre Kraemer's avatar
Pierre Kraemer committed
125
	}
126

Pierre Kraemer's avatar
Pierre Kraemer committed
127
	//if the point is in the volume there is an odd number of intersection with all faces with any direction
128
	return ((countInter % 2) != 0) && ((countInter2 % 2) != 0); //	return (countInter % 2) == 1;
Pierre Kraemer's avatar
Pierre Kraemer committed
129 130 131
}

template <typename PFP>
132
bool isPointInConvexVolume(typename PFP::MAP& map, Vol v, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, const typename PFP::VEC3& point, bool CCW)
Pierre Kraemer's avatar
Pierre Kraemer committed
133 134 135 136
{
	typedef typename PFP::VEC3 VEC3 ;
	typedef typename PFP::REAL REAL;

137 138 139 140
	std::vector<Face> visitedFaces;			// Faces that are traversed
	visitedFaces.reserve(64);
	Face f(v.dart);
	visitedFaces.push_back(f);				// Start with the first face of v
Pierre Kraemer's avatar
Pierre Kraemer committed
141

Pierre Kraemer's avatar
Pierre Kraemer committed
142
	DartMarkerStore<typename PFP::MAP> mark(map);		// Lock a marker
143 144

	for (std::vector<Face>::iterator face = visitedFaces.begin(); face != visitedFaces.end(); ++face)
Pierre Kraemer's avatar
Pierre Kraemer committed
145
	{
146 147
		f = *face;
		if (!mark.isMarked(f))
Pierre Kraemer's avatar
Pierre Kraemer committed
148
		{
149 150 151
			mark.markOrbit(f);

			Geom::Plane3D<REAL> p = facePlane<PFP>(map, f, position);
Pierre Kraemer's avatar
Pierre Kraemer committed
152 153 154 155
			Geom::Orientation3D o3d = p.orient(point);
			if(CCW)
			{
				if(o3d == Geom::OVER)
156
					return false;
Pierre Kraemer's avatar
Pierre Kraemer committed
157 158
			}
			else if(o3d == Geom::UNDER)
159
				return false;
Pierre Kraemer's avatar
Pierre Kraemer committed
160

161 162
			// add all face neighbours to the table
			foreach_adjacent2<EDGE>(map, f, [&] (Face ff)
Pierre Kraemer's avatar
Pierre Kraemer committed
163
			{
164 165 166
				if(!mark.isMarked(ff)) // not already marked
					visitedFaces.push_back(ff) ;
			});
Pierre Kraemer's avatar
Pierre Kraemer committed
167 168 169
		}
	}

170
	return true;
Pierre Kraemer's avatar
Pierre Kraemer committed
171 172 173
}

template <typename PFP>
174
bool isPointInConvexFace(typename PFP::MAP& map, Face f, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, const typename PFP::VEC3& point, bool CCW)
Pierre Kraemer's avatar
Pierre Kraemer committed
175 176 177 178
{
	typedef typename PFP::VEC3 VEC3 ;
	typedef typename PFP::REAL REAL;

179
	Geom::Plane3D<REAL> pl = Geometry::facePlane<PFP>(map, f, position);
Pierre Kraemer's avatar
Pierre Kraemer committed
180 181 182
	Geom::Orientation3D o3d = pl.orient(point);
	if(o3d == Geom::ON)
	{
183
		Traversor2FV<typename PFP::MAP> tfv(map, f) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
184
		for(Vertex v = tfv.begin(); v.dart != tfv.end(); v = tfv.next())
Pierre Kraemer's avatar
Pierre Kraemer committed
185 186
		{
			VEC3 N = pl.normal();
187 188 189
			VEC3 v2(position[map.phi1(v.dart)] - position[v]);
			VEC3 norm2 = N ^ v2;
			Geom::Plane3D<REAL> pl2(norm2, position[v]);
Pierre Kraemer's avatar
Pierre Kraemer committed
190 191 192 193
			o3d = pl2.orient(point);
			if(CCW)
			{
				if(o3d == Geom::UNDER)
194
					return false;
Pierre Kraemer's avatar
Pierre Kraemer committed
195 196
			}
			else if(o3d == Geom::OVER)
197
				return false;
198
		}
199
		return true;
Pierre Kraemer's avatar
Pierre Kraemer committed
200 201 202 203 204 205
	}

	return false;
}

template <typename PFP>
206
bool isPointInConvexFace2D(typename PFP::MAP& map, Face f, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, const typename PFP::VEC3& point, bool CCW )
Pierre Kraemer's avatar
Pierre Kraemer committed
207 208 209 210
{
	typedef typename PFP::VEC3 VEC3 ;
	typedef typename PFP::REAL REAL;

211
// 	CGoGNout << "point " << point << "d " << d << "faceDeg " << map.faceDegree(d) << CGoGNendl;
Pierre Kraemer's avatar
Pierre Kraemer committed
212 213

	Geom::Orientation2D o2d;
214

215
	Traversor2FV<typename PFP::MAP> tfv(map, f) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
216
	for(Vertex v = tfv.begin(); v.dart != tfv.end(); v = tfv.next())
Pierre Kraemer's avatar
Pierre Kraemer committed
217
	{
218
		o2d = Geom::testOrientation2D(point, position[v], position[map.phi1(v.dart)]);
Pierre Kraemer's avatar
Pierre Kraemer committed
219 220 221
		if(CCW)
		{
			if(o2d == Geom::RIGHT)
222
				return false;
Pierre Kraemer's avatar
Pierre Kraemer committed
223 224 225
		}
		else if(o2d == Geom::LEFT)
			return false;
226
	}
Pierre Kraemer's avatar
Pierre Kraemer committed
227

228
	return true;
Pierre Kraemer's avatar
Pierre Kraemer committed
229 230 231
}

template <typename PFP>
232
bool isPointOnEdge(typename PFP::MAP& map, Edge e, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, const typename PFP::VEC3& point)
Pierre Kraemer's avatar
Pierre Kraemer committed
233 234 235 236
{
// 	typedef typename PFP::REAL REAL;
// 	typedef typename PFP::VEC3 VEC3 ;
// 
237
// 	VEC3 v1 = vectorOutOfDart<PFP>(map, d, positions);
Pierre Kraemer's avatar
Pierre Kraemer committed
238 239 240 241 242 243 244
// 	VEC3 v2(point - positions[d]);
// 
// 	v1.normalize();
// 	v2.normalize();
// 
// 	return fabs(REAL(1) - (v1*v2)) < std::numeric_limits<REAL>::min();

245 246
	Dart d = e.dart;

247
	if(
248 249 250
		( isPointOnHalfEdge<PFP>(map,d,position,point) && isPointOnHalfEdge<PFP>(map,map.phi2(d),position,point) ) ||
		isPointOnVertex<PFP>(map,d,position,point) ||
		isPointOnVertex<PFP>(map,map.phi1(d),position,point)
251
	)
Pierre Kraemer's avatar
Pierre Kraemer committed
252
		return true;
253 254
	else
	{
255
		CGoGNout << " point " << point << CGoGNendl;
256
		CGoGNout << " d1 " << position[d] << " d2 " << position[map.phi2(d)] << CGoGNendl;
Pierre Kraemer's avatar
Pierre Kraemer committed
257 258 259 260 261
		return false;
	}
}

template <typename PFP>
262
bool isPointOnHalfEdge(typename PFP::MAP& map, Dart d, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, const typename PFP::VEC3& point)
Pierre Kraemer's avatar
Pierre Kraemer committed
263 264 265 266
{
	typedef typename PFP::REAL REAL;
	typedef typename PFP::VEC3 VEC3;

Sylvain Thery's avatar
Sylvain Thery committed
267
	VEC3 v1 = Algo::Geometry::vectorOutOfDart<PFP>(map, d, position);
268
	VEC3 v2(point - position[d]);
Pierre Kraemer's avatar
Pierre Kraemer committed
269 270 271 272

	v1.normalize();
	v2.normalize();

Sylvain Thery's avatar
Sylvain Thery committed
273
	return std::abs(v1*v2) <= REAL(0.00001);
Pierre Kraemer's avatar
Pierre Kraemer committed
274 275 276
}

template <typename PFP>
277
bool isPointOnVertex(typename PFP::MAP& map, Vertex v, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, const typename PFP::VEC3& point)
Pierre Kraemer's avatar
Pierre Kraemer committed
278
{
279
	return Geom::arePointsEquals(point, position[v]);
Pierre Kraemer's avatar
Pierre Kraemer committed
280 281 282
}

template <typename PFP>
283
bool isConvexFaceInOrIntersectingTetrahedron(typename PFP::MAP& map, Face f, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, const typename PFP::VEC3 points[4], bool CCW)
Pierre Kraemer's avatar
Pierre Kraemer committed
284 285 286
{
	typedef typename PFP::VEC3 VEC3 ;

287
	Traversor2FV<typename PFP::MAP> tfv(map, f) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
288
	for(Vertex v = tfv.begin(); v.dart != tfv.end(); v = tfv.next())
Pierre Kraemer's avatar
Pierre Kraemer committed
289
	{
290
		if(Geom::isPointInTetrahedron(points, position[v], CCW))
Pierre Kraemer's avatar
Pierre Kraemer committed
291
			return true;
292
	}
Pierre Kraemer's avatar
Pierre Kraemer committed
293 294

	VEC3 inter;
Sylvain Thery's avatar
Sylvain Thery committed
295 296 297 298 299 300
	if( intersectionSegmentConvexFace<PFP>(map, f, position, points[0], points[1], inter)
		|| intersectionSegmentConvexFace<PFP>(map, f, position, points[1], points[2], inter)
		|| intersectionSegmentConvexFace<PFP>(map, f, position, points[2], points[0], inter)
		|| intersectionSegmentConvexFace<PFP>(map, f, position, points[0], points[3], inter)
		|| intersectionSegmentConvexFace<PFP>(map, f, position, points[1], points[3], inter)
		|| intersectionSegmentConvexFace<PFP>(map, f, position, points[2], points[3], inter)
Pierre Kraemer's avatar
Pierre Kraemer committed
301 302 303 304 305 306
	)
		return true;

	return false;
}

307
} // namespace Geometry
Pierre Kraemer's avatar
Pierre Kraemer committed
308

309
} // namespace Surface
310

311
} // namspace Algo
Pierre Kraemer's avatar
Pierre Kraemer committed
312

313
} // namespace CGoGN