inclusion.hpp 9.45 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
{
	typedef typename PFP::VEC3 VEC3 ;
	typedef typename PFP::REAL REAL;

	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
142
	for (face = visitedFaces.begin(); face != visitedFaces.end(); ++face)
Pierre Kraemer's avatar
Pierre Kraemer committed
143 144 145
	{
		if (!mark.isMarked(*face))
		{
146
			Geom::Plane3D<REAL> p = facePlane<PFP>(map, *face, position);
Pierre Kraemer's avatar
Pierre Kraemer committed
147 148 149 150
			Geom::Orientation3D o3d = p.orient(point);
			if(CCW)
			{
				if(o3d == Geom::OVER)
151
					return false;
Pierre Kraemer's avatar
Pierre Kraemer committed
152 153
			}
			else if(o3d == Geom::UNDER)
154
				return false;
Pierre Kraemer's avatar
Pierre Kraemer committed
155 156 157 158 159 160 161 162 163 164 165 166 167

			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) ;
		}
	}

168
	return true;
Pierre Kraemer's avatar
Pierre Kraemer committed
169 170 171
}

template <typename PFP>
172
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
173 174 175 176
{
	typedef typename PFP::VEC3 VEC3 ;
	typedef typename PFP::REAL REAL;

177
	Geom::Plane3D<REAL> pl = Geometry::facePlane<PFP>(map, d, position);
Pierre Kraemer's avatar
Pierre Kraemer committed
178 179 180
	Geom::Orientation3D o3d = pl.orient(point);
	if(o3d == Geom::ON)
	{
181 182
		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
183 184
		{
			VEC3 N = pl.normal();
185
			VEC3 v2(position[map.phi1(it)] - position[it]);
Pierre Kraemer's avatar
Pierre Kraemer committed
186
			VEC3 norm2 = N^v2;
187
			Geom::Plane3D<REAL> pl2(norm2, position[it]);
Pierre Kraemer's avatar
Pierre Kraemer committed
188 189 190 191
			o3d = pl2.orient(point);
			if(CCW)
			{
				if(o3d == Geom::UNDER)
192
					return false;
Pierre Kraemer's avatar
Pierre Kraemer committed
193 194
			}
			else if(o3d == Geom::OVER)
195
				return false;
196
		}
197
		return true;
Pierre Kraemer's avatar
Pierre Kraemer committed
198 199 200 201 202 203
	}

	return false;
}

template <typename PFP>
204
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
205 206 207 208
{
	typedef typename PFP::VEC3 VEC3 ;
	typedef typename PFP::REAL REAL;

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

	Geom::Orientation2D o2d;
212 213 214

	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
215
	{
216
		o2d = Geom::testOrientation2D(point, position[it], position[map.phi1(it)]);
Pierre Kraemer's avatar
Pierre Kraemer committed
217 218 219
		if(CCW)
		{
			if(o2d == Geom::RIGHT)
220
				return false;
Pierre Kraemer's avatar
Pierre Kraemer committed
221 222 223
		}
		else if(o2d == Geom::LEFT)
			return false;
224
	}
Pierre Kraemer's avatar
Pierre Kraemer committed
225

226
	return true;
Pierre Kraemer's avatar
Pierre Kraemer committed
227 228 229
}

template <typename PFP>
230
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
231 232 233 234
{
// 	typedef typename PFP::REAL REAL;
// 	typedef typename PFP::VEC3 VEC3 ;
// 
235
// 	VEC3 v1 = vectorOutOfDart<PFP>(map, d, positions);
Pierre Kraemer's avatar
Pierre Kraemer committed
236 237 238 239 240 241 242
// 	VEC3 v2(point - positions[d]);
// 
// 	v1.normalize();
// 	v2.normalize();
// 
// 	return fabs(REAL(1) - (v1*v2)) < std::numeric_limits<REAL>::min();

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

template <typename PFP>
258
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
259 260 261 262
{
	typedef typename PFP::REAL REAL;
	typedef typename PFP::VEC3 VEC3;

263
	VEC3 v1 = vectorOutOfDart<PFP>(map, d, position);
264
	VEC3 v2(point - position[d]);
Pierre Kraemer's avatar
Pierre Kraemer committed
265 266 267 268 269 270 271 272

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

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

template <typename PFP>
273
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
274
{
275
	return Geom::arePointsEquals(point, position[d]);
Pierre Kraemer's avatar
Pierre Kraemer committed
276 277 278
}

template <typename PFP>
279
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
280 281 282
{
	typedef typename PFP::VEC3 VEC3 ;

283 284
	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
285
	{
286
		if(Geom::isPointInTetrahedron(points, position[it], CCW))
Pierre Kraemer's avatar
Pierre Kraemer committed
287
			return true;
288
	}
Pierre Kraemer's avatar
Pierre Kraemer committed
289 290

	VEC3 inter;
291 292 293 294 295 296
	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
297 298 299 300 301 302
	)
		return true;

	return false;
}

303
} // namespace Geometry
Pierre Kraemer's avatar
Pierre Kraemer committed
304

305
} // namespace Surface
306

307
} // namspace Algo
Pierre Kraemer's avatar
Pierre Kraemer committed
308

309
} // namespace CGoGN