inclusion.hpp 9.62 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>
45
bool isConvex(typename PFP::MAP& map, Dart d, const VertexAttribute<typename PFP::VEC3>& position, bool CCW, unsigned int thread)
Pierre Kraemer's avatar
Pierre Kraemer committed
46 47 48 49
{
	//get all the dart of the volume
	std::vector<Dart> vStore;
	FunctorStore fs(vStore);
50
	map.foreach_dart_of_volume(d, fs, thread);
Pierre Kraemer's avatar
Pierre Kraemer committed
51 52 53

	bool convex = true;

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

	return convex;
}

Sylvain Thery's avatar
Sylvain Thery committed
68
// TODO add thread Pameter
Pierre Kraemer's avatar
Pierre Kraemer committed
69
template <typename PFP>
70
bool isPointInVolume(typename PFP::MAP& map, Dart d, const VertexAttribute<typename PFP::VEC3>& position, const typename PFP::VEC3& point)
Pierre Kraemer's avatar
Pierre Kraemer committed
71 72 73 74 75
{
	typedef typename PFP::VEC3 VEC3;

	//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);
Pierre Kraemer's avatar
Pierre Kraemer committed
79
	std::vector<VEC3> interPrec;
Sylvain Thery's avatar
Sylvain Thery committed
80
	interPrec.reserve(16);
81 82 83 84 85 86
	std::vector<Dart> visitedFaces;			// Faces that are traversed
	visitedFaces.reserve(64);
	visitedFaces.push_back(d);				// Start with the face of d
	DartMarkerStore mark(map);
	mark.markOrbit<FACE>(d) ;
	for(unsigned int  iface = 0; iface != visitedFaces.size(); ++iface)
Pierre Kraemer's avatar
Pierre Kraemer committed
87
	{
88 89
		Dart e = visitedFaces[iface];
		VEC3 inter;
90
		bool interRes = intersectionLineConvexFace<PFP>(map, e, position, point, dir, inter);
91
		if (interRes)
Pierre Kraemer's avatar
Pierre Kraemer committed
92
		{
Sylvain Thery's avatar
Sylvain Thery committed
93
			// check if already intersect on same point (a vertex certainly)
Pierre Kraemer's avatar
Pierre Kraemer committed
94
			bool alreadyfound = false;
Sylvain Thery's avatar
Sylvain Thery committed
95 96 97 98 99 100 101
			for(typename std::vector<VEC3>::iterator it = interPrec.begin(); !alreadyfound && it != interPrec.end(); ++it)
			{
				if (Geom::arePointsEquals(*it,inter))
						alreadyfound = true;
			}

			if (!alreadyfound)
Pierre Kraemer's avatar
Pierre Kraemer committed
102
			{
103 104 105 106 107
				float v = dir * (inter-point);
				if (v>0)
					++countInter;
				if (v<0)
					++countInter2;
Sylvain Thery's avatar
Sylvain Thery committed
108
				interPrec.push_back(inter);
Pierre Kraemer's avatar
Pierre Kraemer committed
109 110
			}
		}
111 112 113 114 115 116 117 118 119 120 121 122
		// add all face neighbours to the table
		Dart currentFace = e;
		do
		{
			Dart ee = map.phi2(e) ;
			if(!mark.isMarked(ee)) // not already marked
			{
				visitedFaces.push_back(ee) ;
				mark.markOrbit<FACE>(ee) ;
			}
			e = map.phi1(e) ;
		} while(e != currentFace) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
123
	}
124

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

Pierre Kraemer's avatar
Pierre Kraemer committed
128 129 130
}

template <typename PFP>
131
bool isPointInConvexVolume(typename PFP::MAP& map, Dart d, const VertexAttribute<typename PFP::VEC3>& position, const typename PFP::VEC3& point, bool CCW)
Pierre Kraemer's avatar
Pierre Kraemer committed
132 133 134 135 136 137 138 139 140 141 142 143 144 145 146
{
	typedef typename PFP::VEC3 VEC3 ;
	typedef typename PFP::REAL REAL;

	bool inside = true;
	std::list<Dart> visitedFaces;			// Faces that are traversed
	visitedFaces.push_back(d);				// Start with the face of d
	std::list<Dart>::iterator face;
	VEC3 N;

	DartMarkerStore mark(map);					// Lock a marker
	for (face = visitedFaces.begin(); inside && face != visitedFaces.end(); ++face)
	{
		if (!mark.isMarked(*face))
		{
147
			Geom::Plane3D<REAL> p = facePlane<PFP>(map, *face, position);
Pierre Kraemer's avatar
Pierre Kraemer committed
148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173
			Geom::Orientation3D o3d = p.orient(point);
			if(CCW)
			{
				if(o3d == Geom::OVER)
					inside = false;
			}
			else if(o3d == Geom::UNDER)
				inside = false;

			Dart dNext = *face ;
			do
			{
				mark.mark(dNext);					// Mark
				Dart adj = map.phi2(dNext);			// Get adjacent face
				if (adj != dNext && !mark.isMarked(adj))
					visitedFaces.push_back(adj);	// Add it
				dNext = map.phi1(dNext) ;
			} while(dNext != *face) ;
		}
	}

	//if the point is in the volume there is an odd number of intersection with all faces with any direction
	return inside;
}

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

	bool inside = true;

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

	return false;
}

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

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

	bool convex=true;
	Geom::Orientation2D o2d;
217 218 219

	Traversor2FV<typename PFP::MAP> tfv(map, d) ;
	for(Dart it = tfv.begin(); it != tfv.end(); it = tfv.next())
Pierre Kraemer's avatar
Pierre Kraemer committed
220
	{
221
		o2d = Geom::testOrientation2D(point, position[it], position[map.phi1(it)]);
Pierre Kraemer's avatar
Pierre Kraemer committed
222 223 224 225 226 227 228
		if(CCW)
		{
			if(o2d == Geom::RIGHT)
				convex= false;
		}
		else if(o2d == Geom::LEFT)
			return false;
229
	}
Pierre Kraemer's avatar
Pierre Kraemer committed
230 231 232 233 234

	return convex;
}

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

248
	if(
249 250 251
		( 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)
252
	)
Pierre Kraemer's avatar
Pierre Kraemer committed
253
		return true;
254 255
	else
	{
256
		CGoGNout << " point " << point << CGoGNendl;
257
		CGoGNout << " d1 " << position[d] << " d2 " << position[map.phi2(d)] << CGoGNendl;
Pierre Kraemer's avatar
Pierre Kraemer committed
258 259 260 261 262
		return false;
	}
}

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

268
	VEC3 v1 = vectorOutOfDart<PFP>(map, d, position);
269
	VEC3 v2(point - position[d]);
Pierre Kraemer's avatar
Pierre Kraemer committed
270 271 272 273 274 275 276 277

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

	return (v1*v2) >= -REAL(0.00001);
}

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

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

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

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

	return false;
}

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

310 311
}

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

314
} // namespace CGoGN