Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

source: trunk/source/graphics_reps/src/BooleanProcessor.src@ 927

Last change on this file since 927 was 830, checked in by garnier, 17 years ago
import all except CVS
File size: 69.2 KB

Line
1	/***********************************************************************
2	* *
3	* Name: BooleanProcessor Date: 10.12.99 *
4	* Author: E.Chernyaev Revised: *
5	* *
6	* Function: Internal class for executing boolean operations *
7	* on Polyhedra *
8	* *
9	***********************************************************************/
10
11	#include <CLHEP/Geometry/Plane3D.h>
12
13	using namespace HepGeom;
14
15	#define INITIAL_SIZE 200
16	#define CRAZY_POINT Point3D<double>(-10.e+6, -10.e+6, -10.e+6)
17	#define GRANULARITY 10.e+5;
18
19	#define SWAP(A,B) w = A; A = B; B = w
20
21	#define OP_UNION 0 // Operations
22	#define OP_INTERSECTION 1
23	#define OP_SUBTRACTION 2
24
25	#define OUT_OF_PLANE 0 // Face vs face statuses
26	#define ON_PLANE 1
27	#define INTERSECTION 2
28	#define EDGE 3
29	#define NON_PLANAR_FACE 4
30
31	#define UNKNOWN_FACE 0 // Face statuses
32	#define ORIGINAL_FACE -1
33	#define NEW_FACE -2
34	#define UNSUITABLE_FACE -3
35	#define DEFECTIVE_FACE -4
36
37	// -------------------------------------------- Simplified STL vector ---
38	template<class T>
39	class vector {
40	private:
41	int cur_size, max_size;
42	T * v;
43
44	public:
45	vector(): cur_size(0), max_size(INITIAL_SIZE), v(new T[INITIAL_SIZE]) {}
46	~vector() { delete [] v; }
47
48	void clear() { cur_size = 0; }
49	int size() const { return cur_size; }
50	T & operator[](int i) { return v[i]; }
51	const T & operator[](int i) const { return v[i]; }
52	T & front() { return v[0]; }
53	const T & front() const { return v[0]; }
54	T & back() { return v[cur_size-1]; }
55	const T & back() const { return v[cur_size-1]; }
56	void pop_back() { cur_size--; }
57	void push_back(const T & a) {
58	if (cur_size == max_size) {
59	T * w = v;
60	max_size *= 2;
61	v = new T[max_size];
62	for (int i=0; i<cur_size; i++) v[i] = w[i];
63	v[cur_size++] = a;
64	delete [] w;
65	}else{
66	v[cur_size++] = a;
67	}
68	}
69	};
70
71	// ---------------------------------------------------- Extended node ---
72	class ExtNode {
73	public:
74	Point3D<double> v;
75	int s;
76
77	public:
78	ExtNode(Point3D<double> vertex=Point3D<double>(), int status=0)
79	: v(vertex), s(status) {}
80	~ExtNode() {}
81
82	ExtNode(const ExtNode & node) : v(node.v), s(node.s) {}
83
84	ExtNode & operator=(const ExtNode & node) {
85	v = node.v;
86	s = node.s;
87	return *this;
88	}
89	};
90
91	// ---------------------------------------------------- Extended edge ---
92	class ExtEdge {
93	public:
94	int i1, i2; // end points
95	int iface1; // native face
96	int iface2; // neighbouring face
97	int ivis; // visibility: +1 (visible), -1 (invisible)
98	int inext; // index of next edge
99
100	public:
101	ExtEdge(int k1=0, int k2=0, int kface1=0, int kface2=0, int kvis=0) :
102	i1(k1), i2(k2), iface1(kface1), iface2(kface2), ivis(kvis), inext(0) {}
103
104	~ExtEdge() {};
105
106	ExtEdge(const ExtEdge & edge) :
107	i1(edge.i1), i2(edge.i2), iface1(edge.iface1), iface2(edge.iface2),
108	ivis(edge.ivis), inext(edge.inext) {}
109
110	ExtEdge & operator=(const ExtEdge & edge) {
111	i1 = edge.i1;
112	i2 = edge.i2;
113	iface1 = edge.iface1;
114	iface2 = edge.iface2;
115	ivis = edge.ivis;
116	inext = edge.inext;
117	return *this;
118	}
119
120	void invert() {
121	int w;
122	SWAP(i1, i2);
123	}
124	};
125
126	// ---------------------------------------------------- Extended face ---
127	class ExtFace {
128	public:
129	int iedges[4]; // indices of original edges
130	Plane3D<double> plane; // face plane
131	double rmin[3], rmax[3]; // bounding box
132	int iold; // head of the list of the original edges
133	int inew; // head of the list of the new edges
134	int iprev; // index of previous face
135	int inext; // index of next face
136
137	public:
138	ExtFace(int iedge=0) : iold(iedge), inew(0), iprev(iprev), inext(0) {}
139	~ExtFace() {}
140
141	ExtFace(const ExtFace & face) :
142	plane(face.plane), iold(face.iold), inew(face.inew),
143	iprev(face.iprev), inext(face.inext)
144	{
145	int i;
146	for (i=0; i<4; i++) { iedges[i] = face.iedges[i]; }
147	for (i=0; i<3; i++) { rmin[i] = face.rmin[i]; rmax[i] = face.rmax[i]; }
148	}
149
150	ExtFace & operator=(const ExtFace & face) {
151	int i;
152	for (i=0; i<4; i++) { iedges[i] = face.iedges[i]; }
153	plane = face.plane;
154	for (i=0; i<3; i++) { rmin[i] = face.rmin[i]; rmax[i] = face.rmax[i]; }
155	iold = face.iold;
156	inew = face.inew;
157	iprev = face.iprev;
158	inext = face.inext;
159	return *this;
160	}
161
162	void invert();
163	};
164
165	// ---------------------------------------------------- Global arrays ---
166	static vector<ExtNode> nodes; // vector of nodes
167	static vector<ExtEdge> edges; // vector of edges
168	static vector<ExtFace> faces; // vector of faces
169
170	// ---------------------------------------------------- List of faces ---
171	class FaceList {
172	private:
173	int ihead;
174	int ilast;
175
176	public:
177	FaceList() : ihead(0), ilast(0) {}
178	~FaceList() {}
179
180	void clean() { ihead = 0; ilast = 0; }
181	int front() { return ihead; }
182
183	void push_back(int i) {
184	if (ilast == 0) { ihead = i; } else { faces[ilast].inext = i; }
185	faces[i].iprev = ilast;
186	faces[i].inext = 0;
187	ilast = i;
188	}
189
190	void remove(int i) {
191	if (ihead == i) {
192	ihead = faces[i].inext;
193	}else{
194	faces[faces[i].iprev].inext = faces[i].inext;
195	}
196	if (ilast == i) {
197	ilast = faces[i].iprev;
198	}else{
199	faces[faces[i].inext].iprev = faces[i].iprev;
200	}
201	faces[i].iprev = 0;
202	faces[i].inext = 0;
203	}
204	};
205
206	// --------------------- Polyhedron with extended access to
207	// its members from the BooleanProcessor class ---
208	class ExtPolyhedron : public HepPolyhedron {
209	friend class BooleanProcessor;
210	virtual HepPolyhedron& operator = (const HepPolyhedron& from) {
211	return HepPolyhedron::operator = (from);
212	}
213	};
214
215	// ----------------------------------------- Boolean processor class ---
216	class BooleanProcessor {
217	private:
218	int processor_error; // is set in case of error
219	int operation; // 0 (union), 1 (intersection), 2 (subtraction)
220	int ifaces1, ifaces2; // lists of faces
221	int iout1, iout2; // lists of faces with status "out"
222	int iunk1, iunk2; // lists of faces with status "unknown"
223	double rmin[3], rmax[3]; // intersection of bounding boxes
224	double del; // precision (tolerance)
225
226	FaceList result_faces; // list of accepted faces
227	FaceList suitable_faces; // list of suitable faces
228	FaceList unsuitable_faces; // list of unsuitable faces
229	FaceList unknown_faces; // list of unknown faces
230
231	vector<int> external_contours; // heads of external contours
232	vector<int> internal_contours; // heads of internal contours
233
234	private:
235	void takePolyhedron(const HepPolyhedron & p, double, double, double);
236	double findMinMax();
237	void selectOutsideFaces(int & ifaces, int & iout);
238	int testFaceVsPlane(ExtEdge & edge);
239	void renumberNodes(int & i1, int & i2, int & i3, int & i4);
240	int testEdgeVsEdge(ExtEdge & edge1, ExtEdge & edge2);
241	void removeJunkNodes() { while(nodes.back().s != 0) nodes.pop_back(); }
242	void divideEdge(int & i1, int & i2);
243	void insertEdge(const ExtEdge & edge);
244	void caseII(ExtEdge & edge1, ExtEdge & edge2);
245	void caseIE(ExtEdge & edge1, ExtEdge & edge2);
246	void caseEE(ExtEdge & edge1, ExtEdge & edge2);
247	void testFaceVsFace(int iface1, int iface2);
248	void invertNewEdges(int iface);
249	void checkDoubleEdges(int iface);
250	void assembleFace(int what, int iface);
251	void assembleNewFaces(int what, int ihead);
252	void initiateLists();
253	void assemblePolyhedra();
254	void findABC(double x1, double y1, double x2, double y2,
255	double &a, double &b, double &c) const;
256	int checkDirection(double x, double y) const;
257	int checkIntersection(int ix, int iy, int i1, int i2) const;
258	void mergeContours(int ix, int iy, int kext, int kint);
259	int checkTriangle(int iedge1, int iedge2, int ix, int iy) const;
260	void triangulateContour(int ix, int iy, int ihead);
261	void modifyReference(int iface, int i1, int i2, int iref);
262	void triangulateFace(int iface);
263	HepPolyhedron createPolyhedron();
264
265	public:
266	BooleanProcessor() {}
267	~BooleanProcessor() {}
268
269	HepPolyhedron execute(int op,
270	const HepPolyhedron &a,
271	const HepPolyhedron &b);
272
273	void draw();
274	void draw_edge(int, int);
275	void draw_contour(int, int, int);
276	void draw_faces(int, int, int);
277	void print_face(int);
278	void print_edge(int);
279	int get_processor_error() const {return processor_error;}
280	};
281
282	inline void ExtFace::invert()
283	/***********************************************************************
284	* *
285	* Name: ExtFace::invert() Date: 28.02.00 *
286	* Author: E.Chernyaev Revised: *
287	* *
288	* Function: Invert face *
289	* *
290	***********************************************************************/
291	{
292	int iEprev, iEcur, iEnext;
293
294	iEprev = 0; iEcur = iold;
295	while (iEcur > 0) {
296	edges[iEcur].invert();
297	iEnext = edges[iEcur].inext;
298	edges[iEcur].inext = iEprev;
299	iEprev = iEcur;
300	iEcur = iEnext;
301	}
302	if (iold > 0) iold = iEprev;
303
304	iEprev = 0; iEcur = inew;
305	while (iEcur > 0) {
306	edges[iEcur].invert();
307	iEnext = edges[iEcur].inext;
308	edges[iEcur].inext = iEprev;
309	iEprev = iEcur;
310	iEcur = iEnext;
311	}
312	if (inew > 0) inew = iEprev;
313
314	plane = Plane3D<double>(-plane.a(), -plane.b(), -plane.c(), -plane.d());
315	}
316
317	void BooleanProcessor::takePolyhedron(const HepPolyhedron & p,
318	double dx, double dy, double dz)
319	/***********************************************************************
320	* *
321	* Name: BooleanProcessor::takePolyhedron Date: 16.12.99 *
322	* Author: E.Chernyaev Revised: *
323	* *
324	* Function: Transfer Polyhedron to internal representation *
325	* *
326	***********************************************************************/
327	{
328	int i, k, nnode, iNodes[5], iVis[4], iFaces[4];
329	int dnode = nodes.size() - 1;
330	int dface = faces.size() - 1;
331
332	// S E T N O D E S
333
334	// for (i=1; i <= p.GetNoVertices(); i++) {
335	// nodes.push_back(ExtNode(p.GetVertex(i)));
336	// }
337
338	Point3D<double> ppp;
339	for (i=1; i <= p.GetNoVertices(); i++) {
340	ppp = p.GetVertex(i);
341	ppp.setX(ppp.x()+dx);
342	ppp.setY(ppp.y()+dy);
343	ppp.setZ(ppp.z()+dz);
344	nodes.push_back(ExtNode(ppp));
345	}
346
347	// S E T F A C E S
348
349	for (int iface=1; iface <= p.GetNoFacets(); iface++) {
350	faces.push_back(ExtFace(edges.size()));
351
352	// S E T F A C E N O D E S
353
354	p.GetFacet(iface, nnode, iNodes, iVis, iFaces);
355	for (i=0; i<nnode; i++) {
356	if (iNodes[i] < 1 \|\| iNodes[i] > p.GetNoVertices()) processor_error = 1;
357	if (iFaces[i] < 1 \|\| iFaces[i] > p.GetNoFacets()) processor_error = 1;
358	iNodes[i] += dnode;
359	iFaces[i] += dface;
360	}
361
362	// S E T E D G E S
363
364	iNodes[nnode] = iNodes[0];
365	faces.back().iedges[3] = 0;
366	for (i=0; i<nnode; i++) {
367	faces.back().iedges[i] = edges.size();
368	edges.push_back(ExtEdge(iNodes[i], iNodes[i+1],
369	iface+dface, iFaces[i], iVis[i]));
370	edges.back().inext = edges.size();
371	}
372	edges.back().inext = 0;
373
374	// S E T F A C E M I N - M A X
375
376	for (i=0; i<3; i++) {
377	faces.back().rmin[i] = nodes[iNodes[0]].v[i];
378	faces.back().rmax[i] = nodes[iNodes[0]].v[i];
379	}
380	for (i=1; i<nnode; i++) {
381	for (k=0; k<3; k++) {
382	if (faces.back().rmin[k] > nodes[iNodes[i]].v[k])
383	faces.back().rmin[k] = nodes[iNodes[i]].v[k];
384	if (faces.back().rmax[k] < nodes[iNodes[i]].v[k])
385	faces.back().rmax[k] = nodes[iNodes[i]].v[k];
386	}
387	}
388
389	// S E T F A C E P L A N E
390
391	Normal3D<double> n = (nodes[iNodes[2]].v-nodes[iNodes[0]].v).cross
392	(nodes[iNodes[3]].v-nodes[iNodes[1]].v);
393	Point3D<double> p(0,0,0);
394
395	for (i=0; i<nnode; i++) { p += nodes[iNodes[i]].v; }
396	p *= (1./nnode);
397	faces.back().plane = Plane3D<double>(n.unit(), p);
398
399	// S E T R E F E R E N C E T O T H E N E X T F A C E
400
401	faces.back().inext = faces.size();
402	}
403	faces.back().inext = 0;
404	}
405
406	double BooleanProcessor::findMinMax()
407	/***********************************************************************
408	* *
409	* Name: BooleanProcessor::findMinMax Date: 16.12.99 *
410	* Author: E.Chernyaev Revised: *
411	* *
412	* Function: Find min-max (bounding) boxes for polyhedra *
413	* *
414	***********************************************************************/
415	{
416	if (ifaces1 == 0 \|\| ifaces2 == 0) return 0;
417
418	int i, iface;
419	double rmin1[3], rmax1[3];
420	double rmin2[3], rmax2[3];
421
422	// F I N D B O U N D I N G B O X E S
423
424	for (i=0; i<3; i++) {
425	rmin1[i] = faces[ifaces1].rmin[i];
426	rmax1[i] = faces[ifaces1].rmax[i];
427	rmin2[i] = faces[ifaces2].rmin[i];
428	rmax2[i] = faces[ifaces2].rmax[i];
429	}
430
431	iface = faces[ifaces1].inext;
432	while(iface > 0) {
433	for (i=0; i<3; i++) {
434	if (rmin1[i] > faces[iface].rmin[i]) rmin1[i] = faces[iface].rmin[i];
435	if (rmax1[i] < faces[iface].rmax[i]) rmax1[i] = faces[iface].rmax[i];
436	}
437	iface = faces[iface].inext;
438	}
439
440	iface = faces[ifaces2].inext;
441	while(iface > 0) {
442	for (i=0; i<3; i++) {
443	if (rmin2[i] > faces[iface].rmin[i]) rmin2[i] = faces[iface].rmin[i];
444	if (rmax2[i] < faces[iface].rmax[i]) rmax2[i] = faces[iface].rmax[i];
445	}
446	iface = faces[iface].inext;
447	}
448
449	// F I N D I N T E R S E C T I O N O F B O U N D I N G B O X E S
450
451	for (i=0; i<3; i++) {
452	rmin[i] = (rmin1[i] > rmin2[i]) ? rmin1[i] : rmin2[i];
453	rmax[i] = (rmax1[i] < rmax2[i]) ? rmax1[i] : rmax2[i];
454	}
455
456	// F I N D T O L E R A N C E
457
458	double del1 = 0;
459	double del2 = 0;
460	for (i=0; i<3; i++) {
461	if ((rmax1[i]-rmin1[i]) > del1) del1 = rmax1[i]-rmin1[i];
462	if ((rmax2[i]-rmin2[i]) > del2) del2 = rmax2[i]-rmin2[i];
463	}
464	return ((del1 < del2) ? del1 : del2) / GRANULARITY;
465	}
466
467	void BooleanProcessor::selectOutsideFaces(int & ifaces, int & iout)
468	/***********************************************************************
469	* *
470	* Name: BooleanProcessor::selectOutsideFaces Date: 10.01.00 *
471	* Author: E.Chernyaev Revised: *
472	* *
473	* Function: Preselection of outside faces *
474	* *
475	***********************************************************************/
476	{
477	int i, outflag, iface = ifaces, *prev;
478	Point3D<double> mmbox[8] = { Point3D<double>(rmin[0],rmin[1],rmin[2]),
479	Point3D<double>(rmax[0],rmin[1],rmin[2]),
480	Point3D<double>(rmin[0],rmax[1],rmin[2]),
481	Point3D<double>(rmax[0],rmax[1],rmin[2]),
482	Point3D<double>(rmin[0],rmin[1],rmax[2]),
483	Point3D<double>(rmax[0],rmin[1],rmax[2]),
484	Point3D<double>(rmin[0],rmax[1],rmax[2]),
485	Point3D<double>(rmax[0],rmax[1],rmax[2]) };
486	prev = &ifaces;
487	while (iface > 0) {
488
489	// B O U N D I N G B O X vs B O U N D I N G B O X
490
491	outflag = 0;
492	for (i=0; i<3; i++) {
493	if (faces[iface].rmin[i] > rmax[i] + del) { outflag = 1; break; }
494	if (faces[iface].rmax[i] < rmin[i] - del) { outflag = 1; break; }
495	}
496
497	// B O U N D I N G B O X vs P L A N E
498
499	if (outflag == 0) {
500	int npos = 0, nneg = 0;
501	double d;
502	for (i=0; i<8; i++) {
503	d = faces[iface].plane.distance(mmbox[i]);
504	if (d > +del) npos++;
505	if (d < -del) nneg++;
506	}
507	if (npos == 8 \|\| nneg == 8) outflag = 1;
508	}
509
510	// U P D A T E L I S T S
511
512	if (outflag == 1) {
513	*prev = faces[iface].inext;
514	faces[iface].inext = iout;
515	iout = iface;
516	}else{
517	prev = &faces[iface].inext;
518	}
519	iface = *prev;
520	}
521	}
522
523	int BooleanProcessor::testFaceVsPlane(ExtEdge & edge)
524	/***********************************************************************
525	* *
526	* Name: BooleanProcessor::testFaceVsPlane Date: 19.01.00 *
527	* Author: E.Chernyaev Revised: *
528	* *
529	* Function: Find result of intersection of face by plane *
530	* *
531	***********************************************************************/
532	{
533	int iface = edge.iface1;
534	Plane3D<double> plane = faces[edge.iface2].plane;
535	int i, nnode, npos = 0, nneg = 0, nzer = 0;
536	double dd[5];
537
538	// F I N D D I S T A N C E S
539
540	nnode = (faces[iface].iedges[3] == 0) ? 3 : 4;
541	for (i=0; i<nnode; i++) {
542	dd[i] = plane.distance(nodes[edges[faces[iface].iedges[i]].i1].v);
543	if (dd[i] > del) {
544	npos++;
545	}else if (dd[i] < -del) {
546	nneg++;
547	}else{
548	nzer++; dd[i] = 0;
549	}
550	}
551
552	// S O M E S I M P L E C A S E S ( N O I N T E R S E C T I O N )
553
554	if (npos == nnode \|\| nneg == nnode) return OUT_OF_PLANE;
555	if (nzer == 1 && nneg == 0) return OUT_OF_PLANE;
556	if (nzer == 1 && npos == 0) return OUT_OF_PLANE;
557	if (nzer == nnode) return ON_PLANE;
558	if (nzer == 3) return NON_PLANAR_FACE;
559
560	// F I N D I N T E R S E C T I O N
561
562	int ie1 = 0, ie2 = 0, s1 = 0, s2 = 0, status, nint = 0;
563	enum { PLUS_MINUS, MINUS_PLUS, ZERO_ZERO, ZERO_PLUS, ZERO_MINUS };
564
565	dd[nnode] = dd[0];
566	for (i=0; i<nnode; i++) {
567	if (dd[i] > 0) {
568	if (dd[i+1] >= 0) continue;
569	status = PLUS_MINUS;
570	}else if (dd[i] < 0) {
571	if (dd[i+1] <= 0) continue;
572	status = MINUS_PLUS;
573	}else{
574	status = ZERO_ZERO;
575	if (dd[i+1] > 0) status = ZERO_PLUS;
576	if (dd[i+1] < 0) status = ZERO_MINUS;
577	}
578	switch (nint) {
579	case 0:
580	ie1 = i; s1 = status; nint++; break;
581	case 1:
582	ie2 = i; s2 = status; nint++; break;
583	default:
584	return NON_PLANAR_FACE;
585	}
586	}
587	if (nint != 2) return NON_PLANAR_FACE;
588
589	// F O R M I N T E R S E C T I O N S E G M E N T
590
591	if (s1 != ZERO_ZERO && s2 != ZERO_ZERO) {
592	if (s1 == s2) return NON_PLANAR_FACE;
593	int iedge, i1 = 0, i2 = 0, ii[2];
594	double d1 = 0., d2 = 0., dd ;
595	ii[0] = ie1; ii[1] = ie2;
596	for (i=0; i<2; i++) {
597	iedge = faces[iface].iedges[ii[i]];
598	while (iedge > 0) {
599	i1 = edges[iedge].i1;
600	i2 = edges[iedge].i2;
601	d1 = plane.distance(nodes[i1].v);
602	d2 = plane.distance(nodes[i2].v);
603	if (d1 > del) {
604	if (d2 < -del) { ii[i] = nodes.size(); break; } // +-
605	}else if (d1 < -del) {
606	if (d2 > del) { ii[i] = nodes.size(); break; } // -+
607	}else{
608	ii[i] = i1; break; // 0+ or 0-
609	}
610	iedge = edges[iedge].inext;
611	}
612	if (ii[i] == nodes.size()) {
613	dd = d2-d1; d1 = d1/dd; d2 = d2/dd;
614	nodes.push_back(ExtNode(d2nodes[i1].v-d1nodes[i2].v, iedge));
615	}
616	}
617	edge.inext = 0;
618	if (s1 == MINUS_PLUS \|\| s1 == ZERO_PLUS) {
619	edge.i1 = ii[1];
620	edge.i2 = ii[0];
621	}else{
622	edge.i1 = ii[0];
623	edge.i2 = ii[1];
624	}
625	return INTERSECTION;
626	}else{
627	if (npos == nneg) return NON_PLANAR_FACE;
628	edge.inext = (s1 == ZERO_ZERO) ? ie1+1 : ie2+1;
629	if (s1 == ZERO_PLUS \|\| s2 == ZERO_MINUS) {
630	edge.i1 = edges[faces[iface].iedges[ie2]].i1;
631	edge.i2 = edges[faces[iface].iedges[ie1]].i1;
632	}else{
633	edge.i1 = edges[faces[iface].iedges[ie1]].i1;
634	edge.i2 = edges[faces[iface].iedges[ie2]].i1;
635	}
636	return EDGE;
637	}
638	}
639
640	void BooleanProcessor::renumberNodes(int & i1, int & i2, int & i3, int & i4)
641	/***********************************************************************
642	* *
643	* Name: BooleanProcessor::renumberNodes Date: 19.01.00 *
644	* Author: E.Chernyaev Revised: *
645	* *
646	* Function: Renumber nodes and remove last temporary node. *
647	* Remark: In principal this routine can be replaced just *
648	* with i1 = i2; *
649	* Removal of temporary nodes provides additional control *
650	* on number of nodes, that is very useful for debugging. *
651	* *
652	***********************************************************************/
653	{
654	if (i1 == i2) return;
655	if (nodes[i1].s == 0 \|\| nodes.back().s == 0) { i1 = i2; return; }
656
657	int ilast = nodes.size()-1;
658	if (i1 == ilast) { i1 = i2; nodes.pop_back(); return; }
659	if (i2 == ilast) { i2 = i1; }
660	if (i3 == ilast) { i3 = i1; }
661	if (i4 == ilast) { i4 = i1; }
662	nodes[i1] = nodes.back(); i1 = i2; nodes.pop_back();
663	}
664
665	int BooleanProcessor::testEdgeVsEdge(ExtEdge & edge1, ExtEdge & edge2)
666	/***********************************************************************
667	* *
668	* Name: BooleanProcessor::testEdgeVsEdge Date: 19.01.00 *
669	* Author: E.Chernyaev Revised: *
670	* *
671	* Function: Find common part of two edges *
672	* *
673	***********************************************************************/
674	{
675	int i, ii = 0;
676	double d, dd = 0.;
677
678	for (i=0; i<3; i++) {
679	d = nodes[edge1.i1].v[i]-nodes[edge1.i2].v[i];
680	if (d < 0.) d = -d;
681	if (d > dd) { dd = d; ii = i; }
682	}
683	double t1 = nodes[edge1.i1].v[ii];
684	double t2 = nodes[edge1.i2].v[ii];
685	double t3 = nodes[edge2.i1].v[ii];
686	double t4 = nodes[edge2.i2].v[ii];
687	if (t2-t1 < 0.) { t1 = -t1; t2 = -t2; t3 = -t3; t4 = -t4; }
688
689	if (t3 <= t1+del \|\| t4 >= t2-del) return 0;
690	if (t3 > t2+del) {
691	renumberNodes(edge2.i1, edge1.i2, edge1.i1, edge2.i2);
692	}else if (t3 < t2-del) {
693	renumberNodes(edge1.i2, edge2.i1, edge1.i1, edge2.i2);
694	}
695	if (t4 < t1-del) {
696	renumberNodes(edge2.i2, edge1.i1, edge1.i2, edge2.i1);
697	}else if (t4 > t1+del) {
698	renumberNodes(edge1.i1, edge2.i2, edge1.i2, edge2.i1);
699	}
700	return 1;
701	}
702
703	void BooleanProcessor::divideEdge(int & i1, int & i2)
704	/***********************************************************************
705	* *
706	* Name: BooleanProcessor::divideEdge Date: 24.01.00 *
707	* Author: E.Chernyaev Revised: *
708	* *
709	* Function: Unify the nodes and divide edge on two parts by the node. *
710	* *
711	***********************************************************************/
712	{
713	int iedges[2];
714	iedges[0] = nodes[i1].s;
715	iedges[1] = nodes[i2].s;
716
717	// U N I F Y N O D E S
718
719	if (i1 < i2) { i2 = i1; }
720	else if (i1 > i2) { i1 = i2; }
721	else { iedges[1] = 0; }
722	if (iedges[0] == iedges[1]) return;
723
724	int ie1, ie2, inode = i1;
725	nodes[inode].s = 0;
726	for (int i=0; i<2; i++) {
727
728	// F I N D C O R R E S P O N D I N G E D G E
729
730	if ((ie1 = iedges[i]) == 0) continue;
731	ie2 = faces[edges[ie1].iface2].iedges[0];
732	while (ie2 > 0) {
733	if (edges[ie2].i1 == edges[ie1].i2 &&
734	edges[ie2].i2 == edges[ie1].i1) break;
735	ie2 = edges[ie2].inext;
736	}
737
738	// D I V I D E E D G E S
739
740	edges.push_back(edges[ie1]);
741	edges[ie1].inext = edges.size() - 1;
742	edges[ie1].i2 = inode;
743	edges.back().i1 = inode;
744
745	edges.push_back(edges[ie2]);
746	edges[ie2].inext = edges.size() - 1;
747	edges[ie2].i2 = inode;
748	edges.back().i1 = inode;
749	}
750	}
751
752	void BooleanProcessor::insertEdge(const ExtEdge & edge)
753	/***********************************************************************
754	* *
755	* Name: BooleanProcessor::insertEdge Date: 24.01.00 *
756	* Author: E.Chernyaev Revised: *
757	* *
758	* Function: Insert edge to the list of new edges *
759	* *
760	***********************************************************************/
761	{
762	int iface = edge.iface1;
763	edges.push_back(edge);
764	edges.back().inext = faces[iface].inew;
765	faces[iface].inew = edges.size() - 1;
766	}
767
768	void BooleanProcessor::caseII(ExtEdge & edge1, ExtEdge & edge2)
769	/***********************************************************************
770	* *
771	* Name: BooleanProcessor::caseII Date: 19.01.00 *
772	* Author: E.Chernyaev Revised: *
773	* *
774	* Function: Intersection/Intersection case *
775	* *
776	***********************************************************************/
777	{
778	divideEdge(edge1.i1, edge2.i2);
779	divideEdge(edge1.i2, edge2.i1);
780	edge1.ivis = 1;
781	edge2.ivis = 1;
782	insertEdge(edge1);
783	insertEdge(edge2);
784	}
785
786	void BooleanProcessor::caseIE(ExtEdge &, ExtEdge &)
787	/***********************************************************************
788	* *
789	* Name: BooleanProcessor::caseIE Date: 19.01.00 *
790	* Author: E.Chernyaev Revised: *
791	* *
792	* Function: Intersection/Edge-touch case *
793	* *
794	***********************************************************************/
795	{
796	processor_error = 1;
797	std::cout
798	<< "BooleanProcessor::caseIE : unimplemented case"
799	<< std::endl;
800	}
801
802	void BooleanProcessor::caseEE(ExtEdge &, ExtEdge &)
803	/***********************************************************************
804	* *
805	* Name: BooleanProcessor::caseEE Date: 19.01.00 *
806	* Author: E.Chernyaev Revised: *
807	* *
808	* Function: Edge-touch/Edge-touch case *
809	* *
810	***********************************************************************/
811	{
812	processor_error = 1;
813	std::cout
814	<< "BooleanProcessor::caseEE : unimplemented case"
815	<< std::endl;
816	}
817
818	void BooleanProcessor::testFaceVsFace(int iface1, int iface2)
819	/***********************************************************************
820	* *
821	* Name: BooleanProcessor::testFaceVsFace Date: 11.01.00 *
822	* Author: E.Chernyaev Revised: *
823	* *
824	* Function: Find result (an edge) of intersection of two faces *
825	* *
826	***********************************************************************/
827	{
828	ExtEdge edge1, edge2;
829	int irep1, irep2;
830
831	// M I N - M A X
832
833	for (int i=0; i<3; i++) {
834	if (faces[iface1].rmin[i] > faces[iface2].rmax[i] + del) return;
835	if (faces[iface1].rmax[i] < faces[iface2].rmin[i] - del) return;
836	}
837
838	// F A C E - 1 vs P L A N E - 2
839
840	edge1.iface1 = iface1;
841	edge1.iface2 = iface2;
842	irep1 = testFaceVsPlane(edge1);
843	if (irep1 == OUT_OF_PLANE \|\| irep1 == ON_PLANE) {
844	removeJunkNodes();
845	return;
846	}
847
848	// F A C E - 2 vs P L A N E - 1
849
850	edge2.iface1 = iface2;
851	edge2.iface2 = iface1;
852	irep2 = testFaceVsPlane(edge2);
853	if (irep2 == OUT_OF_PLANE \|\| irep2 == ON_PLANE) {
854	removeJunkNodes();
855	return;
856	}
857
858	// C H E C K F O R N O N P L A N A R F A C E
859
860	if (irep1 == NON_PLANAR_FACE \|\| irep2 == NON_PLANAR_FACE) {
861	removeJunkNodes();
862	return;
863	}
864
865	// F I N D I N T E R S E C T I O N P A R T
866
867	if (testEdgeVsEdge(edge1, edge2) == 0) return;
868
869	// C O N S I D E R D I F F E R E N T C A S E S
870
871	if (irep1 == INTERSECTION && irep2 == INTERSECTION) caseII(edge1, edge2);
872	if (irep1 == INTERSECTION && irep2 == EDGE) caseIE(edge1, edge2);
873	if (irep1 == EDGE && irep2 == INTERSECTION) caseIE(edge2, edge1);
874	if (irep1 == EDGE && irep2 == EDGE) caseEE(edge1, edge2);
875	removeJunkNodes();
876
877	}
878
879	void BooleanProcessor::invertNewEdges(int iface)
880	/***********************************************************************
881	* *
882	* Name: BooleanProcessor::invertNewEdges Date: 04.02.00 *
883	* Author: E.Chernyaev Revised: *
884	* *
885	* Function: Invert direction of new edges *
886	* *
887	***********************************************************************/
888	{
889	int iedge = faces[iface].inew;
890	while (iedge > 0) {
891	edges[iedge].invert();
892	iedge = edges[iedge].inext;
893	}
894	}
895
896	void BooleanProcessor::checkDoubleEdges(int)
897	/***********************************************************************
898	* *
899	* Name: BooleanProcessor::checkDoubleEdges Date: 04.02.00 *
900	* Author: E.Chernyaev Revised: *
901	* *
902	* Function: Eliminate duplication of edges *
903	* *
904	***********************************************************************/
905	{
906
907	}
908
909	void BooleanProcessor::assembleFace(int what, int iface)
910	/***********************************************************************
911	* *
912	* Name: BooleanProcessor::assembleFace Date: 19.02.00 *
913	* Author: E.Chernyaev Revised: *
914	* *
915	* Function: Assemble face *
916	* *
917	***********************************************************************/
918	{
919	// A S S E M B L E N E W F A C E
920
921	int ihead; // head of the list of edges for new face
922	int icur; // current edge in the list - last edge inserted to the list
923	int *ilink; // pointer to the current link
924	int ifirst; // first node of a contour
925	int *i; // pointer to the index of the current edge in a loop
926	int ioldflag=0; // is set if an edge from iold has been taken
927
928	#define INSERT_EDGE_TO_THE_LIST(A) \
929	ilink = A; ilink = &edges[A].inext; ilink = 0
930
931	ilink = &ihead;
932	for(;;) {
933	if (faces[iface].inew == 0) break;
934
935	// S T A R T N E W C O N T O U R
936
937	icur = faces[iface].inew;
938	faces[iface].inew = edges[icur].inext;
939	INSERT_EDGE_TO_THE_LIST(icur);
940	ifirst = edges[icur].i1;
941
942	// C O N S T R U C T T H E C O N T O U R
943
944	for (;;) {
945	i = &faces[iface].inew;
946	while(*i > 0) {
947	if (edges[*i].i1 == edges[icur].i2) break;
948	i = &edges[*i].inext;
949	}
950	if (*i == 0) {
951	i = &faces[iface].iold;
952	while(*i > 0) {
953	if (edges[*i].i1 == edges[icur].i2) ioldflag = 1;
954	if (edges[*i].i1 == edges[icur].i2) break;
955	i = &edges[*i].inext;
956	}
957	}
958	if (*i > 0) {
959	icur = *i;
960	*i = edges[icur].inext;
961	INSERT_EDGE_TO_THE_LIST(icur);
962	if (edges[icur].i2 == ifirst) { break; } else { continue; }
963	}else{
964	processor_error = 1;
965	std::cerr
966	<< "BooleanProcessor::assembleFace(" << iface << ") : "
967	<< "could not find next edge of the contour"
968	<< std::endl;
969	faces[iface].inew = DEFECTIVE_FACE;
970	return;
971	}
972	}
973	}
974
975	// C H E C K O R I G I N A L C O N T O U R
976
977	int iedge;
978	iedge = faces[iface].iold;
979	if (what == 0 && ioldflag == 0 && iedge > 0) {
980	for (;;) {
981	if (edges[iedge].inext > 0) {
982	if (edges[iedge].i2 == edges[edges[iedge].inext].i1) {
983	iedge = edges[iedge].inext;
984	}else{
985	break;
986	}
987	}else{
988	if (edges[iedge].i2 == edges[faces[iface].iold].i1) {
989	edges[iedge].inext = ihead; // set new face
990	return;
991	}else{
992	break;
993	}
994	}
995	}
996	}
997
998	// M A R K U N S U I T A B L E N E I G H B O U R I N G F A C E S
999
1000	int iface2;
1001	iedge = faces[iface].iold;
1002	while(iedge > 0) {
1003	iface2 = edges[iedge].iface2;
1004	if (faces[iface2].inew == 0) faces[iface2].inew = UNSUITABLE_FACE;
1005	iedge = edges[iedge].inext;
1006	}
1007	faces[iface].iold = ihead; // set new face
1008	}
1009
1010	void BooleanProcessor::assembleNewFaces(int what, int ihead)
1011	/***********************************************************************
1012	* *
1013	* Name: BooleanProcessor::assembleNewFaces Date: 30.01.00 *
1014	* Author: E.Chernyaev Revised: *
1015	* *
1016	* Function: Assemble internal or external parts of faces *
1017	* *
1018	***********************************************************************/
1019	{
1020	int iface = ihead;
1021	while(iface > 0) {
1022	if (faces[iface].inew > 0) {
1023	if (what != 0) invertNewEdges(iface);
1024	checkDoubleEdges(iface);
1025	assembleFace(what, iface);
1026	faces[iface].inew =
1027	(faces[iface].iold == 0) ? UNSUITABLE_FACE : NEW_FACE;
1028	}
1029	iface = faces[iface].inext;
1030	}
1031	}
1032
1033	void BooleanProcessor::initiateLists()
1034	/***********************************************************************
1035	* *
1036	* Name: BooleanProcessor::initiateLists Date: 28.02.00 *
1037	* Author: E.Chernyaev Revised: *
1038	* *
1039	* Function: Initiate lists of faces. *
1040	* *
1041	***********************************************************************/
1042	{
1043	int i, iface;
1044
1045	// R E S E T L I S T S O F F A C E S
1046
1047	result_faces.clean();
1048	suitable_faces.clean();
1049	unsuitable_faces.clean();
1050	unknown_faces.clean();
1051
1052	// I N I T I A T E T H E L I S T S
1053
1054	iface = iout1;
1055	while (iface > 0) {
1056	i = iface;
1057	iface = faces[i].inext;
1058	if (operation == OP_INTERSECTION) {
1059	unsuitable_faces.push_back(i);
1060	faces[i].inew = UNSUITABLE_FACE;
1061	}else{
1062	suitable_faces.push_back(i);
1063	faces[i].inew = ORIGINAL_FACE;
1064	}
1065	}
1066	iface = iout2;
1067	while (iface > 0) {
1068	i = iface;
1069	iface = faces[i].inext;
1070	if (operation == OP_UNION) {
1071	suitable_faces.push_back(i);
1072	faces[i].inew = ORIGINAL_FACE;
1073	}else{
1074	unsuitable_faces.push_back(i);
1075	faces[i].inew = UNSUITABLE_FACE;
1076	}
1077	}
1078
1079	iface = iunk1;
1080	while (iface > 0) {
1081	i = iface;
1082	iface = faces[i].inext;
1083	unknown_faces.push_back(i);
1084	}
1085	iface = iunk2;
1086	while (iface > 0) {
1087	i = iface;
1088	iface = faces[i].inext;
1089	if (operation == OP_SUBTRACTION) faces[i].invert();
1090	unknown_faces.push_back(i);
1091	}
1092
1093	iface = ifaces1;
1094	while (iface > 0) {
1095	i = iface;
1096	iface = faces[i].inext;
1097	switch(faces[i].inew) {
1098	case UNKNOWN_FACE:
1099	unknown_faces.push_back(i);
1100	break;
1101	case ORIGINAL_FACE: case NEW_FACE:
1102	suitable_faces.push_back(i);
1103	break;
1104	case UNSUITABLE_FACE:
1105	unsuitable_faces.push_back(i);
1106	break;
1107	default:
1108	faces[i].iprev = 0;
1109	faces[i].inext = 0;
1110	break;
1111	}
1112	}
1113	iface = ifaces2;
1114	while (iface > 0) {
1115	i = iface;
1116	iface = faces[i].inext;
1117	if (operation == OP_SUBTRACTION) faces[i].invert();
1118	switch(faces[i].inew) {
1119	case UNKNOWN_FACE:
1120	unknown_faces.push_back(i);
1121	break;
1122	case ORIGINAL_FACE: case NEW_FACE:
1123	suitable_faces.push_back(i);
1124	break;
1125	case UNSUITABLE_FACE:
1126	unsuitable_faces.push_back(i);
1127	break;
1128	default:
1129	faces[i].iprev = 0;
1130	faces[i].inext = 0;
1131	break;
1132	}
1133	}
1134	ifaces1 = ifaces2 = iout1 = iout2 = iunk1 = iunk2 = 0;
1135	}
1136
1137	void BooleanProcessor::assemblePolyhedra()
1138	/***********************************************************************
1139	* *
1140	* Name: BooleanProcessor::assemblePolyhedra() Date: 10.12.99 *
1141	* Author: E.Chernyaev Revised: *
1142	* *
1143	* Function: Collect suitable faces and remove unsuitable ones. *
1144	* *
1145	***********************************************************************/
1146	{
1147	int i, iedge, iface;
1148
1149	// L O O P A L O N G S U I T A B L E F A C E S
1150
1151	iface = suitable_faces.front();
1152	while(iface > 0) {
1153	i = iface;
1154	iedge = faces[i].iold;
1155	while(iedge > 0) {
1156	iface = edges[iedge].iface2;
1157	if (faces[iface].inew == UNKNOWN_FACE) {
1158	unknown_faces.remove(iface);
1159	suitable_faces.push_back(iface);
1160	faces[iface].inew = ORIGINAL_FACE;
1161	}
1162	iedge = edges[iedge].inext;
1163	}
1164	iface = faces[i].inext;
1165	suitable_faces.remove(i);
1166	result_faces.push_back(i);
1167	}
1168	if (unknown_faces.front() == 0) return;
1169
1170	// L O O P A L O N G U N S U I T A B L E F A C E S
1171
1172	iface = unsuitable_faces.front();
1173	while(iface > 0) {
1174	i = iface;
1175	iedge = faces[i].iold;
1176	while(iedge > 0) {
1177	iface = edges[iedge].iface2;
1178	if (faces[iface].inew == UNKNOWN_FACE) {
1179	unknown_faces.remove(iface);
1180	unsuitable_faces.push_back(iface);
1181	faces[iface].inew = UNSUITABLE_FACE;
1182	}
1183	iedge = edges[iedge].inext;
1184	}
1185	iface = faces[i].inext;
1186	unsuitable_faces.remove(i);
1187	}
1188	}
1189
1190	inline void
1191	BooleanProcessor::findABC(double x1, double y1, double x2, double y2,
1192	double &a, double &b, double &c) const
1193	/***********************************************************************
1194	* *
1195	* Name: BooleanProcessor::findABC Date: 07.03.00 *
1196	* Author: E.Chernyaev Revised: *
1197	* *
1198	* Function: Find line equation Ax+By+C=0 *
1199	* *
1200	***********************************************************************/
1201	{
1202	double w;
1203	a = y1 - y2;
1204	b = x2 - x1;
1205	w = std::abs(a)+std::abs(b);
1206	a /= w;
1207	b /= w;
1208	c = -(ax2 + by2);
1209	}
1210
1211	int BooleanProcessor::checkDirection(double x, double y) const
1212	/***********************************************************************
1213	* *
1214	* Name: BooleanProcessor::checkDirection Date: 06.03.00 *
1215	* Author: E.Chernyaev Revised: *
1216	* *
1217	* Function: Check direction of line 1-4 *
1218	* *
1219	***********************************************************************/
1220	{
1221	double a1, b1, c1, a2, b2, c2, d1, d2;
1222
1223	// T E S T L I N E 1 - 4 V S E X T E R N A L C O N T O U R
1224
1225	findABC(x[0], y[0], x[1], y[1], a1, b1, c1);
1226	findABC(x[1], y[1], x[2], y[2], a2, b2, c2);
1227	d1 = a1x[4] + b1y[4] + c1;
1228	d2 = a2x[4] + b2y[4] + c2;
1229	if (d1 <= del && d2 <= del) return 1;
1230	if (! (d1 > del && d2 > del)) {
1231	if ( a1x[2] + b1y[2] + c1 >= -del) return 1;
1232	}
1233
1234	// T E S T L I N E 1 - 4 V S I N T E R N A L C O N T O U R
1235
1236	findABC(x[3], y[3], x[4], y[4], a1, b1, c1);
1237	findABC(x[4], y[4], x[5], y[5], a2, b2, c2);
1238	d1 = a1x[1] + b1y[1] + c1;
1239	d2 = a2x[1] + b2y[1] + c2;
1240	if (d1 <= del && d2 <= del) return 1;
1241	if (!(d1 > del && d2 > del)) {
1242	if ( a1x[5] + b1y[5] + c1 >= -del) return 1;
1243	}
1244	return 0;
1245	}
1246
1247	int BooleanProcessor::checkIntersection(int ix, int iy, int i1, int i2) const
1248	/***********************************************************************
1249	* *
1250	* Name: BooleanProcessor::checkDirection Date: 06.03.00 *
1251	* Author: E.Chernyaev Revised: *
1252	* *
1253	* Function: Check line i1-i2 on intersection with contours *
1254	* *
1255	***********************************************************************/
1256	{
1257	// F I N D L I N E E Q U A T I O N
1258
1259	double x1, y1, x2, y2, a1, b1, c1;
1260	x1 = nodes[i1].v[ix];
1261	y1 = nodes[i1].v[iy];
1262	x2 = nodes[i2].v[ix];
1263	y2 = nodes[i2].v[iy];
1264	findABC(x1, y1, x2, y2, a1, b1, c1);
1265
1266	// L O O P A L O N G E X T E R N A L C O N T O U R S
1267
1268	int icontour, iedge, k1, k2;
1269	double x3, y3, x4, y4, a2, b2, c2, d1, d2;
1270	for(icontour=0; icontour<external_contours.size(); icontour++) {
1271	iedge = external_contours[icontour];
1272	while(iedge > 0) {
1273	k1 = edges[iedge].i1;
1274	k2 = edges[iedge].i2;
1275	iedge = edges[iedge].inext;
1276	if (k1 == i1 \|\| k2 == i1) continue;
1277	if (k1 == i2 \|\| k2 == i2) continue;
1278	x3 = nodes[k1].v[ix];
1279	y3 = nodes[k1].v[iy];
1280	x4 = nodes[k2].v[ix];
1281	y4 = nodes[k2].v[iy];
1282	d1 = a1x3 + b1y3 + c1;
1283	d2 = a1x4 + b1y4 + c1;
1284	if (d1 > del && d2 > del) continue;
1285	if (d1 < -del && d2 < -del) continue;
1286
1287	findABC(x3, y3, x4, y4, a2, b2, c2);
1288	d1 = a2x1 + b2y1 + c2;
1289	d2 = a2x2 + b2y2 + c2;
1290	if (d1 > del && d2 > del) continue;
1291	if (d1 < -del && d2 < -del) continue;
1292	return 1;
1293	}
1294	}
1295
1296	// L O O P A L O N G E X T E R N A L C O N T O U R S
1297
1298	for(icontour=0; icontour<internal_contours.size(); icontour++) {
1299	iedge = internal_contours[icontour];
1300	while(iedge > 0) {
1301	k1 = edges[iedge].i1;
1302	k2 = edges[iedge].i2;
1303	iedge = edges[iedge].inext;
1304	if (k1 == i1 \|\| k2 == i1) continue;
1305	if (k1 == i2 \|\| k2 == i2) continue;
1306	x3 = nodes[k1].v[ix];
1307	y3 = nodes[k1].v[iy];
1308	x4 = nodes[k2].v[ix];
1309	y4 = nodes[k2].v[iy];
1310	d1 = a1x3 + b1y3 + c1;
1311	d2 = a1x4 + b1y4 + c1;
1312	if (d1 > del && d2 > del) continue;
1313	if (d1 < -del && d2 < -del) continue;
1314
1315	findABC(x3, y3, x4, y4, a2, b2, c2);
1316	d1 = a2x1 + b2y1 + c2;
1317	d2 = a2x2 + b2y2 + c2;
1318	if (d1 > del && d2 > del) continue;
1319	if (d1 < -del && d2 < -del) continue;
1320	return 1;
1321	}
1322	}
1323	return 0;
1324	}
1325
1326	void BooleanProcessor::mergeContours(int ix, int iy, int kext, int kint)
1327	/***********************************************************************
1328	* *
1329	* Name: BooleanProcessor::mergeContours Date: 06.03.00 *
1330	* Author: E.Chernyaev Revised: *
1331	* *
1332	* Function: Attemp to merge internal contour with external one *
1333	* *
1334	***********************************************************************/
1335	{
1336	int i1ext, i2ext, i1int, i2int, i, k[6];
1337	double x[6], y[6];
1338
1339	// L O O P A L O N G E X T E R N A L C O N T O U R
1340
1341	i1ext = external_contours[kext];
1342	while (i1ext > 0) {
1343	i2ext = edges[i1ext].inext;
1344	if (i2ext == 0) i2ext = external_contours[kext];
1345	k[0] = edges[i1ext].i1;
1346	k[1] = edges[i1ext].i2;
1347	k[2] = edges[i2ext].i2;
1348	for (i=0; i<3; i++) {
1349	x[i] = nodes[k[i]].v[ix];
1350	y[i] = nodes[k[i]].v[iy];
1351	}
1352
1353	// L O O P A L O N G I N T E R N A L C O N T O U R
1354
1355	i1int = internal_contours[kint];
1356	while (i1int > 0) {
1357	i2int = edges[i1int].inext;
1358	if (i2int == 0) i2int = internal_contours[kint];
1359	k[3] = edges[i1int].i1;
1360	k[4] = edges[i1int].i2;
1361	k[5] = edges[i2int].i2;
1362	for (i=3; i<6; i++) {
1363	x[i] = nodes[k[i]].v[ix];
1364	y[i] = nodes[k[i]].v[iy];
1365	}
1366
1367	// T E S T L I N E K1 - K4
1368	// I F O K T H E N M E R G E C O N T O U R S
1369
1370	if (checkDirection(x, y) == 0) {
1371	if (checkIntersection(ix, iy, k[1], k[4]) == 0) {
1372	i = i1int;
1373	for(;;) {
1374	if (edges[i].inext == 0) {
1375	edges[i].inext = internal_contours[kint];
1376	internal_contours[kint] = 0;
1377	break;
1378	}else{
1379	i = edges[i].inext;
1380	}
1381	}
1382	i = edges[i1int].iface1;
1383	edges.push_back(ExtEdge(k[1], k[4], i, -(edges.size()+1), -1));
1384	edges.back().inext = i2int;
1385	edges.push_back(ExtEdge(k[4], k[1], i, -(edges.size()-1), -1));
1386	edges.back().inext = edges[i1ext].inext;
1387	edges[i1ext].inext = edges.size()-2;
1388	edges[i1int].inext = edges.size()-1;
1389	return;
1390	}
1391	}
1392	i1int = edges[i1int].inext;
1393	}
1394	i1ext = edges[i1ext].inext;
1395	}
1396	}
1397
1398	int
1399	BooleanProcessor::checkTriangle(int iedge1, int iedge2, int ix, int iy) const
1400	/***********************************************************************
1401	* *
1402	* Name: BooleanProcessor::checkTriangle Date: 08.03.00 *
1403	* Author: E.Chernyaev Revised: *
1404	* *
1405	* Function: Check triangle for correctness *
1406	* *
1407	***********************************************************************/
1408	{
1409	int k[3];
1410	double x[3], y[3];
1411	double a1, b1, c1;
1412
1413	k[0] = edges[iedge1].i1;
1414	k[1] = edges[iedge1].i2;
1415	k[2] = edges[iedge2].i2;
1416	for (int i=0; i<3; i++) {
1417	x[i] = nodes[k[i]].v[ix];
1418	y[i] = nodes[k[i]].v[iy];
1419	}
1420
1421	// C H E C K P R I N C I P A L C O R R E C T N E S S
1422
1423	findABC(x[2], y[2], x[0], y[0], a1, b1, c1);
1424	if (a1x[1]+b1y[1]+c1 <= 0.1*del) return 1;
1425
1426	// C H E C K T H A T T H E R E I S N O P O I N T S I N S I D E
1427
1428	int inode, iedge;
1429	double a2, b2, c2, a3, b3, c3;
1430
1431	findABC(x[0], y[0], x[1], y[1], a2, b2, c2);
1432	findABC(x[1], y[1], x[2], y[2], a3, b3, c3);
1433	iedge = iedge2;
1434	for (;;) {
1435	iedge = edges[iedge].inext;
1436	if (edges[iedge].inext == iedge1) return 0;
1437	inode = edges[iedge].i2;
1438	if (inode == k[0]) continue;
1439	if (inode == k[1]) continue;
1440	if (inode == k[2]) continue;
1441	x[1] = nodes[inode].v[ix];
1442	y[1] = nodes[inode].v[iy];
1443	if (a1x[1]+b1y[1]+c1 < -0.1*del) continue;
1444	if (a2x[1]+b2y[1]+c2 < -0.1*del) continue;
1445	if (a3x[1]+b3y[1]+c3 < -0.1*del) continue;
1446	return 1;
1447	}
1448	}
1449
1450	void BooleanProcessor::triangulateContour(int ix, int iy, int ihead)
1451	/***********************************************************************
1452	* *
1453	* Name: BooleanProcessor::triangulateContour Date: 06.03.00 *
1454	* Author: E.Chernyaev Revised: *
1455	* *
1456	* Function: Triangulate external contour *
1457	* *
1458	***********************************************************************/
1459	{
1460
1461	//std::cout << "Next Countour" << std::endl;
1462	//int draw_flag = 0;
1463	//if (draw_flag) draw_contour(5, 3, ihead);
1464
1465	// C L O S E C O N T O U R
1466
1467	int ipnext = ihead, nnode = 1;
1468	for (;;) {
1469	if (edges[ipnext].inext > 0) {
1470	ipnext = edges[ipnext].inext;
1471	nnode++;
1472	}else{
1473	edges[ipnext].inext = ihead;
1474	break;
1475	}
1476	}
1477
1478	// L O O P A L O N G C O N T O U R
1479
1480	int iedge1, iedge2, iedge3, istart = 0;
1481	for (;;) {
1482	iedge1 = edges[ipnext].inext;
1483	iedge2 = edges[iedge1].inext;
1484	if (istart == 0) {
1485	istart = iedge1;
1486	if (nnode <= 3) {
1487	iedge3 = edges[iedge2].inext;
1488	edges[iedge1].iface1 = faces.size();
1489	edges[iedge2].iface1 = faces.size();
1490	edges[iedge3].iface1 = faces.size();
1491	edges[iedge3].inext = 0;
1492	faces.push_back(ExtFace());
1493	faces.back().iold = iedge1;
1494	faces.back().inew = ORIGINAL_FACE;
1495
1496	//if (draw_flag) draw_contour(4, 2, iedge1);
1497
1498	break;
1499	}
1500	}else if (istart == iedge1) {
1501	processor_error = 1;
1502	std::cerr
1503	<< "BooleanProcessor::triangulateContour : "
1504	<< "could not generate a triangle (infinite loop)"
1505	<< std::endl;
1506	break;
1507	}
1508
1509	// C H E C K C O R E C T N E S S O F T H E T R I A N G L E
1510
1511	if (checkTriangle(iedge1,iedge2,ix,iy) != 0) {
1512	ipnext = edges[ipnext].inext;
1513	continue;
1514	}
1515
1516	// M O D I F Y C O N T O U R
1517
1518	int i1 = edges[iedge1].i1;
1519	int i3 = edges[iedge2].i2;
1520	int iface1 = edges[iedge1].iface1;
1521	int iface2 = faces.size();
1522
1523	edges[ipnext].inext = edges.size();
1524	edges.push_back(ExtEdge(i1, i3, iface1, -(edges.size()+1), -1));
1525	edges.back().inext = edges[iedge2].inext;
1526
1527	// A D D N E W T R I A N G L E T O T H E L I S T
1528
1529	edges[iedge2].inext = edges.size();
1530	edges.push_back(ExtEdge(i3, i1, iface2, -(edges.size()-1), -1));
1531	faces.push_back(ExtFace());
1532	faces.back().iold = iedge1;
1533	faces.back().inew = ORIGINAL_FACE;
1534	edges[iedge1].iface1 = iface2;
1535	edges[iedge2].iface1 = iface2;
1536	ipnext = edges[ipnext].inext;
1537	istart = 0;
1538	nnode--;
1539
1540	//if (draw_flag) draw_contour(4, 2, iedge1);
1541
1542	}
1543	}
1544
1545	void BooleanProcessor::modifyReference(int iface, int i1, int i2, int iref)
1546	/***********************************************************************
1547	* *
1548	* Name: BooleanProcessor::modifyReference Date: 13.03.00 *
1549	* Author: E.Chernyaev Revised: *
1550	* *
1551	* Function: Modify reference to the neighbouring face *
1552	* *
1553	***********************************************************************/
1554	{
1555	int iedge = faces[iface].iold;
1556	while (iedge > 0) {
1557	if (edges[iedge].i1 == i2 && edges[iedge].i2 == i1) {
1558	edges[iedge].iface2 = iref;
1559	return;
1560	}
1561	iedge = edges[iedge].inext;
1562	}
1563	processor_error = 1;
1564	std::cerr
1565	<< "BooleanProcessor::modifyReference : could not find the edge, "
1566	<< "iface=" << iface << ", i1,i2=" << i1 << "," << i2 << ", iref=" << iref
1567	<< std::endl;
1568	}
1569
1570	void BooleanProcessor::triangulateFace(int iface)
1571	/***********************************************************************
1572	* *
1573	* Name: BooleanProcessor::triangulateFace Date: 02.03.00 *
1574	* Author: E.Chernyaev Revised: *
1575	* *
1576	* Function: Triangulation of an extended face *
1577	* *
1578	***********************************************************************/
1579	{
1580
1581	// F I N D M A X C O M P O N E N T O F T H E N O R M A L
1582	// S E T IX, IY, IZ
1583
1584	Normal3D<double> normal = faces[iface].plane.normal();
1585	int ix, iy, iz = 0;
1586	if (std::abs(normal[1]) > std::abs(normal[iz])) iz = 1;
1587	if (std::abs(normal[2]) > std::abs(normal[iz])) iz = 2;
1588	if (normal[iz] > 0) {
1589	ix = (iz+1)%3; iy = (ix+1)%3;
1590	}else{
1591	iy = (iz+1)%3; ix = (iy+1)%3;
1592	}
1593
1594	// F I L L L I S T S O F C O N T O U R S
1595
1596	external_contours.clear();
1597	internal_contours.clear();
1598	double z;
1599	int i1, i2, ifirst, iedge, icontour = faces[iface].iold;
1600	while (icontour > 0) {
1601	iedge = icontour;
1602	ifirst = edges[iedge].i1;
1603	z = 0.0;
1604	for(;;) {
1605	if (iedge > 0) {
1606	i1 = edges[iedge].i1;
1607	i2 = edges[iedge].i2;
1608	z += nodes[i1].v[ix]nodes[i2].v[iy]-nodes[i2].v[ix]nodes[i1].v[iy];
1609	if (ifirst != i2) {
1610	iedge = edges[iedge].inext;
1611	continue;
1612	}else{
1613	if (z > del*del) {
1614	external_contours.push_back(icontour);
1615	}else if (z < -del*del) {
1616	internal_contours.push_back(icontour);
1617	}else{
1618	processor_error = 1;
1619	std::cerr
1620	<< "BooleanProcessor::triangulateFace : too small contour"
1621	<< std::endl;
1622	}
1623	icontour = edges[iedge].inext;
1624	edges[iedge].inext = 0;
1625	break;
1626	}
1627	}else{
1628	processor_error = 1;
1629	std::cerr
1630	<< "BooleanProcessor::triangulateFace : broken contour"
1631	<< std::endl;
1632	icontour = 0;
1633	break;
1634	}
1635	}
1636	}
1637
1638	// G E T R I D O F I N T E R N A L C O N T O U R S
1639
1640	int kint, kext;
1641	for (kint=0; kint < internal_contours.size(); kint++) {
1642	for (kext=0; kext < external_contours.size(); kext++) {
1643	mergeContours(ix, iy, kext, kint);
1644	if (internal_contours[kint] == 0) break;
1645	}
1646	if (kext == external_contours.size()) {
1647	processor_error = 1;
1648	std::cerr
1649	<< "BooleanProcessor::triangulateFace : "
1650	<< "could not merge internal contour " << kint
1651	<< std::endl;
1652	}
1653	}
1654
1655	// T R I A N G U L A T E C O N T O U R S
1656
1657	int nface = faces.size();
1658	for (kext=0; kext < external_contours.size(); kext++) {
1659	triangulateContour(ix, iy, external_contours[kext]);
1660	}
1661	faces[iface].inew = UNSUITABLE_FACE;
1662
1663	// M O D I F Y R E F E R E N C E S
1664
1665	for (int ifa=nface; ifa<faces.size(); ifa++) {
1666	iedge = faces[ifa].iold;
1667	while (iedge > 0) {
1668	if (edges[iedge].iface1 != ifa) {
1669	processor_error = 1;
1670	std::cerr
1671	<< "BooleanProcessor::triangulateFace : wrong reference to itself, "
1672	<< "iface=" << ifa << ", iface1=" << edges[iedge].iface1
1673	<< std::endl;
1674	}else if (edges[iedge].iface2 > 0) {
1675	modifyReference(edges[iedge].iface2,
1676	edges[iedge].i1, edges[iedge].i2, ifa);
1677	}else if (edges[iedge].iface2 < 0) {
1678	edges[iedge].iface2 = edges[-edges[iedge].iface2].iface1;
1679	}
1680	iedge = edges[iedge].inext;
1681	}
1682	}
1683	}
1684
1685	HepPolyhedron BooleanProcessor::createPolyhedron()
1686	/***********************************************************************
1687	* *
1688	* Name: BooleanProcessor::createPolyhedron() Date: 14.03.00 *
1689	* Author: E.Chernyaev Revised: *
1690	* *
1691	* Function: Create HepPolyhedron. *
1692	* *
1693	***********************************************************************/
1694	{
1695	int i, iedge, nnode = 0, nface = 0;
1696
1697	// R E N U M E R A T E N O D E S A N D F A C E S
1698
1699	for (i=1; i<nodes.size(); i++) nodes[i].s = 0;
1700
1701	for (i=1; i<faces.size(); i++) {
1702	if (faces[i].inew == ORIGINAL_FACE) {
1703	faces[i].inew = ++nface;
1704	iedge = faces[i].iold;
1705	while (iedge > 0) {
1706	nodes[edges[iedge].i1].s = 1;
1707	iedge = edges[iedge].inext;
1708	}
1709	}else{
1710	faces[i].inew = 0;
1711	}
1712	}
1713
1714	for (i=1; i<nodes.size(); i++) {
1715	if (nodes[i].s == 1) nodes[i].s = ++nnode;
1716	}
1717
1718	// A L L O C A T E M E M O R Y
1719
1720	ExtPolyhedron polyhedron;
1721	if (nface == 0) return polyhedron;
1722	polyhedron.AllocateMemory(nnode, nface);
1723
1724	// S E T N O D E S
1725
1726	for (i=1; i<nodes.size(); i++) {
1727	if (nodes[i].s != 0) polyhedron.pV[nodes[i].s] = nodes[i].v;
1728	}
1729
1730	// S E T F A C E S
1731
1732	int k, v[4], f[4];
1733	for (i=1; i<faces.size(); i++) {
1734	if (faces[i].inew == 0) continue;
1735	v[3] = f[3] = k = 0;
1736	iedge = faces[i].iold;
1737	while (iedge > 0) {
1738	if (k > 3) {
1739	std::cerr
1740	<< "BooleanProcessor::createPolyhedron : too many edges"
1741	<< std::endl;
1742	break;
1743	}
1744	v[k] = nodes[edges[iedge].i1].s;
1745	if (edges[iedge].ivis < 0) v[k] = -v[k];
1746	f[k] = faces[edges[iedge].iface2].inew;
1747	iedge = edges[iedge].inext;
1748	k++;
1749	}
1750	if (k < 3) {
1751	std::cerr
1752	<< "BooleanProcessor::createPolyhedron : "
1753	<< "face has only " << k << " edges"
1754	<< std::endl;
1755	}
1756	polyhedron.pF[faces[i].inew] =
1757	G4Facet(v[0],f[0], v[1],f[1], v[2],f[2], v[3],f[3]);
1758	}
1759	return polyhedron;
1760	}
1761
1762	HepPolyhedron BooleanProcessor::execute(int op,
1763	const HepPolyhedron & a,
1764	const HepPolyhedron & b)
1765	/***********************************************************************
1766	* *
1767	* Name: BooleanProcessor::execute Date: 10.12.99 *
1768	* Author: E.Chernyaev Revised: *
1769	* *
1770	* Function: Execute boolean operation. *
1771	* *
1772	***********************************************************************/
1773	{
1774	static int ishift = 0;
1775	static double shift[8][3] = {
1776	{ 31, 23, 17},
1777	{ -31, -23, -17},
1778	{ -23, 17, 31},
1779	{ 23, -17, -31},
1780	{ -17, -31, 23},
1781	{ 17, 31, -23},
1782	{ 31, -23, 17},
1783	{ -31, 23, -17}
1784	};
1785
1786	// I N I T I A T E P R O C E S S O R
1787
1788	processor_error = 0;
1789	operation = op;
1790	nodes.clear(); nodes.push_back(CRAZY_POINT);
1791	edges.clear(); edges.push_back(ExtEdge());
1792	faces.clear(); faces.push_back(ExtFace());
1793
1794	// T A K E P O L Y H E D R A
1795
1796	ifaces1 = faces.size(); takePolyhedron(a,0,0,0);
1797	ifaces2 = faces.size(); takePolyhedron(b,0,0,0);
1798
1799	if (processor_error) { // corrapted polyhedron
1800	std::cerr
1801	<< "BooleanProcessor: corrapted input polyhedron"
1802	<< std::endl;
1803	return HepPolyhedron();
1804	}
1805	if (ifaces1 == ifaces2) { // a is empty
1806	switch (operation) {
1807	case OP_UNION:
1808	return b;
1809	case OP_INTERSECTION:
1810	std::cerr
1811	<< "BooleanProcessor: intersection with empty polyhedron"
1812	<< std::endl;
1813	return HepPolyhedron();
1814	case OP_SUBTRACTION:
1815	std::cerr
1816	<< "BooleanProcessor: subtraction from empty polyhedron"
1817	<< std::endl;
1818	return HepPolyhedron();
1819	}
1820	}
1821	if (ifaces2 == faces.size()) { // b is empty
1822	switch (operation) {
1823	case OP_UNION:
1824	return a;
1825	case OP_INTERSECTION:
1826	std::cerr
1827	<< "BooleanProcessor: intersection with empty polyhedron"
1828	<< std::endl;
1829	return HepPolyhedron();
1830	case OP_SUBTRACTION:
1831	return a;
1832	}
1833	}
1834
1835	// S E T I N I T I A L M I N - M A X A N D T O L E R A N C E
1836
1837	del = findMinMax();
1838
1839	// W O R K A R O U N D T O A V O I D I E A N D E E
1840
1841	double ddxx = del*shift[ishift][0];
1842	double ddyy = del*shift[ishift][1];
1843	double ddzz = del*shift[ishift][2];
1844	ishift++; if (ishift == 8) ishift = 0;
1845
1846	operation = op;
1847	nodes.clear(); nodes.push_back(CRAZY_POINT);
1848	edges.clear(); edges.push_back(ExtEdge());
1849	faces.clear(); faces.push_back(ExtFace());
1850
1851	ifaces1 = faces.size(); takePolyhedron(a,0,0,0);
1852	ifaces2 = faces.size(); takePolyhedron(b,ddxx,ddyy,ddzz);
1853
1854	del = findMinMax();
1855
1856	// P R E S E L E C T O U T S I D E F A C E S
1857
1858	iout1 = iout2 = 0;
1859	selectOutsideFaces(ifaces1, iout1);
1860	selectOutsideFaces(ifaces2, iout2);
1861
1862	// P R E S E L E C T N O I N T E R S E C T I O N F A C E S
1863
1864	int ifa1, ifa2;
1865	iunk1 = iunk2 = 0;
1866	if (iout1 != 0 \|\| iout2 != 0) {
1867	for(;;) {
1868	ifa1 = iunk1;
1869	ifa2 = iunk2;
1870	selectOutsideFaces(ifaces1, iunk1);
1871	selectOutsideFaces(ifaces2, iunk2);
1872	if (iunk1 == ifa1 && iunk2 == ifa2) break;
1873	findMinMax();
1874	}
1875	}
1876
1877	#define PROCESSOR_ERROR \
1878	std::cerr << "BooleanProcessor: boolean operation failed" << std::endl;\
1879	return a;
1880
1881	// F I N D N E W E D G E S
1882
1883	if (ifaces1 != 0 && ifaces2 != 0 ) {
1884	ifa1 = ifaces1;
1885	while (ifa1 > 0) {
1886	ifa2 = ifaces2;
1887	while (ifa2 > 0) {
1888	testFaceVsFace(ifa1, ifa2);
1889	ifa2 = faces[ifa2].inext;
1890	}
1891	ifa1 = faces[ifa1].inext;
1892	}
1893	}
1894	if (processor_error) { PROCESSOR_ERROR }
1895
1896	// C O N S T R U C T N E W F A C E S
1897
1898	assembleNewFaces((operation == OP_INTERSECTION) ? 1 : 0, ifaces1);
1899	if (processor_error) { PROCESSOR_ERROR }
1900	assembleNewFaces((operation == OP_UNION) ? 0 : 1, ifaces2);
1901	if (processor_error) { PROCESSOR_ERROR }
1902
1903	// A S S E M B L E S U I T A B L E F A C E S
1904
1905	initiateLists();
1906	for (;;) {
1907	assemblePolyhedra();
1908	if (unknown_faces.front() != 0) {
1909	processor_error = 1;
1910	std::cerr
1911	<< "BooleanProcessor::execute : unknown faces !!!"
1912	<< std::endl;
1913	}
1914	break;
1915	}
1916	if (processor_error) { PROCESSOR_ERROR }
1917
1918	// T R I A N G U L A T E A C C E P T E D F A C E S
1919
1920	ifa1 = result_faces.front();
1921	while (ifa1 > 0) {
1922	ifa2 = ifa1;
1923	ifa1 = faces[ifa2].inext;
1924	if (faces[ifa2].inew == NEW_FACE) triangulateFace(ifa2);
1925	if (processor_error) { PROCESSOR_ERROR }
1926	}
1927
1928	// C R E A T E P O L Y H E D R O N
1929
1930	return createPolyhedron();
1931	}
1932
1933
1934	//#include <cfortran.h>
1935	//#include <higz.h>
1936	//#include "zbuf.h"
1937	//void BooleanProcessor::draw()
1938	/***********************************************************************
1939	* *
1940	* Name: BooleanProcessor::draw Date: 10.12.99 *
1941	* Author: E.Chernyaev Revised: *
1942	* *
1943	* Function: Draw *
1944	* *
1945	***********************************************************************/
1946	/*
1947	{
1948	int II;
1949	int icol, i1, i2, iedge, iface, ilist[4];
1950	float p1[3], p2[3];
1951
1952	ilist[0] = ifaces1;
1953	ilist[1] = ifaces2;
1954	ilist[2] = iout1;
1955	ilist[3] = iout2;
1956
1957	for (int i=0; i<4; i++) {
1958
1959	if (i == 0) cout << "========= Ifaces_1" << endl;
1960	if (i == 1) cout << "========= Ifaces_2" << endl;
1961	if (i == 2) cout << "========= Iout_1" << endl;
1962	if (i == 3) cout << "========= Iout_2" << endl;
1963
1964	icol = i+1;
1965	iface = ilist[i];
1966	while (iface > 0) {
1967
1968	cout << "iface = " << iface << endl;
1969	cout << "--- iold" << endl;
1970
1971	iedge = faces[iface].iold;
1972	icol = 2;
1973
1974	while (iedge > 0) {
1975
1976	cout << " iegde = " << iedge
1977	<< " i1,i2 =" << edges[iedge].i1 << "," << edges[iedge].i2
1978	<< " iface1,iface2 = "
1979	<< edges[iedge].iface1 << "," << edges[iedge].iface2
1980	<< endl;
1981
1982	i1 = edges[iedge].i1;
1983	p1[0] = nodes[i1].v.x();
1984	p1[1] = nodes[i1].v.y();
1985	p1[2] = nodes[i1].v.z();
1986	IHWTON(p1,p1);
1987	i2 = edges[iedge].i2;
1988	p2[0] = nodes[i2].v.x();
1989	p2[1] = nodes[i2].v.y();
1990	p2[2] = nodes[i2].v.z();
1991	IHWTON(p2,p2);
1992	// icol = (edges[iedge].ivis > 0) ? 1 : 2;
1993	IHZLIN(icol,p1[0],p1[1],p1[2], p2[0],p2[1],p2[2]);
1994	iedge = edges[iedge].inext;
1995	}
1996
1997	cout << "--- inew" << endl;
1998
1999	iedge = faces[iface].inew;
2000	icol = 3;
2001
2002	while (iedge > 0) {
2003
2004	cout << " iegde = " << iedge
2005	<< " i1,i2 =" << edges[iedge].i1 << "," << edges[iedge].i2
2006	<< " iface1,iface2 = "
2007	<< edges[iedge].iface1 << "," << edges[iedge].iface2
2008	<< endl;
2009
2010	i1 = edges[iedge].i1;
2011	p1[0] = nodes[i1].v.x();
2012	p1[1] = nodes[i1].v.y();
2013	p1[2] = nodes[i1].v.z();
2014	IHWTON(p1,p1);
2015	i2 = edges[iedge].i2;
2016	p2[0] = nodes[i2].v.x();
2017	p2[1] = nodes[i2].v.y();
2018	p2[2] = nodes[i2].v.z();
2019	IHWTON(p2,p2);
2020	// icol = (edges[iedge].ivis > 0) ? 1 : 2;
2021	IHZLIN(icol,p1[0],p1[1],p1[2], p2[0],p2[1],p2[2]);
2022	iedge = edges[iedge].inext;
2023	}
2024	iface = faces[iface].inext;
2025
2026	IHZTOX(0,100,100);
2027	ixupdwi(0);
2028	cin >> II;
2029	ixclrwi();
2030	IHZCLE(0);
2031	}
2032	}
2033	}
2034	*/
2035
2036	/*
2037	//--------------------------------------------------------------------
2038	void
2039	BooleanProcessor::draw_edge(int icol, int iedge) {
2040	int i1, i2;
2041	float p1[3], p2[3];
2042
2043	i1 = edges[iedge].i1;
2044	p1[0] = nodes[i1].v.x();
2045	p1[1] = nodes[i1].v.y();
2046	p1[2] = nodes[i1].v.z();
2047	IHWTON(p1,p1);
2048	i2 = edges[iedge].i2;
2049	p2[0] = nodes[i2].v.x();
2050	p2[1] = nodes[i2].v.y();
2051	p2[2] = nodes[i2].v.z();
2052	IHWTON(p2,p2);
2053	IHZLIN(icol,p1[0],p1[1],p1[2], p2[0],p2[1],p2[2]);
2054	}
2055
2056	//--------------------------------------------------------------------
2057	void
2058	BooleanProcessor::draw_contour(int i1col, int i2col, int ihead) {
2059	int iedge, icol;
2060	iedge = ihead;
2061	while (iedge > 0) {
2062	icol = (edges[iedge].ivis > 0) ? i1col : i2col;
2063	draw_edge(icol, iedge);
2064	iedge = edges[iedge].inext;
2065	}
2066
2067	IHZTOX(0,100,100);
2068	ixupdwi(0);
2069
2070	int i;
2071	std::cin >> i;
2072	}
2073
2074	//--------------------------------------------------------------------
2075	void
2076	BooleanProcessor::print_face(int iface) {
2077	cout.precision(3);
2078	cout << "\n====== Face N " << iface << endl;
2079	cout << "iedges[4] = "
2080	<< faces[iface].iedges[0] << ", "
2081	<< faces[iface].iedges[1] << ", "
2082	<< faces[iface].iedges[2] << ", "
2083	<< faces[iface].iedges[3] << endl;
2084	cout << "rmin[3] = "
2085	<< faces[iface].rmin[0] << ", "
2086	<< faces[iface].rmin[1] << ", "
2087	<< faces[iface].rmin[2] << endl;
2088	cout << "rmax[3] = "
2089	<< faces[iface].rmax[0] << ", "
2090	<< faces[iface].rmax[1] << ", "
2091	<< faces[iface].rmax[2] << endl;
2092	cout << "iprev,inext = "
2093	<< faces[iface].iprev << ", "
2094	<< faces[iface].inext << endl;
2095	cout << "iold = " << faces[iface].iold << endl;
2096	for(int i = faces[iface].iold; i != 0;) {
2097	print_edge(i);
2098	i = edges[abs(i)].inext;
2099	}
2100
2101	cout << "inew = ";
2102	switch (faces[iface].inew) {
2103	case UNKNOWN_FACE:
2104	cout << "UNKNOWN_FACE" << endl;
2105	break;
2106	case ORIGINAL_FACE:
2107	cout << "ORIGINAL_FACE" << endl;
2108	break;
2109	case NEW_FACE:
2110	cout << "NEW_FACE" << endl;
2111	break;
2112	case UNSUITABLE_FACE:
2113	cout << "UNSUITABLE_FACE" << endl;
2114	break;
2115	case DEFECTIVE_FACE:
2116	cout << "DEFECTIVE_FACE" << endl;
2117	break;
2118	default:
2119	cout << faces[iface].inew << endl;
2120	for(int k = faces[iface].inew; k != 0;) {
2121	print_edge(k);
2122	k = edges[abs(k)].inext;
2123	}
2124	}
2125	}
2126
2127	//--------------------------------------------------------------------
2128	void
2129	BooleanProcessor::print_edge(int iedge) {
2130	cout << "==== Edge N " << iedge << endl;
2131	int i = std::abs(iedge);
2132	int i1 = edges[i].i1;
2133	int i2 = edges[i].i2;
2134	cout << "node[" << i1 << "] = "
2135	<< nodes[i1].v.x() << ", "
2136	<< nodes[i1].v.y() << ", "
2137	<< nodes[i1].v.z() << endl;
2138
2139	cout << "node[" << i2 << "] = "
2140	<< nodes[i2].v.x() << ", "
2141	<< nodes[i2].v.y() << ", "
2142	<< nodes[i2].v.z() << endl;
2143
2144	cout << "iface1,iface2,ivis,inext = "
2145	<< edges[i].iface1 << ", "
2146	<< edges[i].iface2 << ", "
2147	<< edges[i].ivis << ", "
2148	<< edges[i].inext << endl;
2149	}
2150	*/

Note: See TracBrowser for help on using the repository browser.

Download in other formats: