1 /*
2  * facet.cc -- Utility classes for shaded surface plotting
3  *
4  * This file is part of ePiX, a C++ library for creating high-quality
5  * figures in LaTeX
6  *
7  * Version 1.2.2
8  * Last Change: December 06, 2007
9  */
10 
11 /*
12  * Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007
13  * Andrew D. Hwang <rot 13 nujnat at zngupf dot ubylpebff dot rqh>
14  * Department of Mathematics and Computer Science
15  * College of the Holy Cross
16  * Worcester, MA, 01610-2395, USA
17  */
18 
19 /*
20  * ePiX is free software; you can redistribute it and/or modify it
21  * under the terms of the GNU General Public License as published by
22  * the Free Software Foundation; either version 2 of the License, or
23  * (at your option) any later version.
24  *
25  * ePiX is distributed in the hope that it will be useful, but WITHOUT
26  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
27  * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
28  * License for more details.
29  *
30  * You should have received a copy of the GNU General Public License
31  * along with ePiX; if not, write to the Free Software Foundation, Inc.,
32  * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
33  */
34 #include <cmath>
35 
36 #include "constants.h"
37 #include "errors.h"
38 #include "functions.h"
39 
40 #include "frame.h"
41 #include "camera.h"
42 
43 #include "Color.h"
44 #include "paint_style.h"
45 
46 #include "facet.h"
47 
48 namespace ePiX {
49 
facet(P f (double,double),double u0,double v0,double du,double dv,const unsigned int N1,const unsigned int N2)50   facet::facet(P f(double, double),
51 	       double u0, double v0,
52 	       double du, double dv,
53 	       const unsigned int N1, const unsigned int N2)
54     : m_tint(the_paint_style().fill_color()),
55       m_line(the_paint_style().line_pen()),
56       m_fill(the_paint_style().fill_flag()),
57       pt1(f(u0,         v0)),
58       pt2(f(u0 + N1*du, v0)),
59       pt3(f(u0 + N1*du, v0 + N2*dv)),
60       pt4(f(u0,         v0 + N2*dv)),
61       center(0.25*(pt1 + pt2 + pt3 + pt4)),
62       direction(center-cam().viewpt()),
63       distance(norm(direction))
64   {
65     perp = ((pt1 - center)*(pt2 - center));
66 
67     if (norm(perp) < EPIX_EPSILON)
68       perp = (pt3 - center)*(pt4 - center);
69 
70     perp *= recip(norm(perp));
71 
72     // bottom edge
73     for (unsigned int i=0; i<N1; ++i)
74       bd.pt(f(u0 + i*du, v0));
75 
76     // right edge
77     for (unsigned int j=0; j<N2; ++j)
78       bd.pt(f(u0 + N1*du, v0 + j*dv));
79 
80     // top edge (backward)
81     for (unsigned int i=0; i<N1; ++i)
82       bd.pt(f(u0 + (N1-i)*du, v0 + N2*dv));
83 
84     // left edge (downward)
85     for (unsigned int j=0; j<N2; ++j)
86       bd.pt(f(u0, v0 + (N2-j)*dv));
87 
88     bd.close().fill(!m_tint.is_unset());
89   }
90 
91 
92   // facet constructor for f(double, double, double)
facet(P f (double,double,double),double u0,double v0,double w0,double du,double dv,double dw,const unsigned int N1,const unsigned int N2)93   facet::facet(P f(double, double, double),
94 	       double u0, double v0, double w0,
95 	       double du, double dv, double dw,
96 	       const unsigned int N1, const unsigned int N2)
97     : m_tint(the_paint_style().fill_color()),
98       m_line(the_paint_style().line_pen()),
99       m_fill(the_paint_style().fill_flag()),
100       pt1(f(u0, v0, w0))
101   {
102     if (du == 0) // (y,z)
103       {
104 	// bottom edge
105 	for (unsigned int i=0; i<N1; ++i)
106 	  bd.pt(f(u0, v0 + i*dv, w0));
107 
108 	// right edge
109 	for (unsigned int j=0; j<N2; ++j)
110 	  bd.pt(f(u0, v0 + N1*dv, w0 + j*dw));
111 
112 	// top edge (backward)
113 	for (unsigned int i=0; i<N1; ++i)
114 	  bd.pt(f(u0, v0 + (N1-i)*dv, w0 + N2*dw));
115 
116 	// left edge (downward)
117 	for (unsigned int j=0; j<N2; ++j)
118 	  bd.pt(f(u0, v0, w0 + (N2-j)*dw));
119 
120 	// use corners to approximate distance to camera
121 	pt2 = f(u0, v0 + N1*dv, w0);
122 	pt3 = f(u0, v0 + N1*dv, w0 + N2*dw);
123 	pt4 = f(u0, v0,         w0 + N2*dw);
124       }
125 
126     else if (dv == 0) // (x,z)
127       {
128 	// bottom edge
129 	for (unsigned int i=0; i<N1; ++i)
130 	  bd.pt(f(u0 + i*du, v0, w0));
131 
132 	// right edge
133 	for (unsigned int j=0; j<N2; ++j)
134 	  bd.pt(f(u0 + N1*du, v0, w0 + j*dw));
135 
136 	// top edge (backward)
137 	for (unsigned int i=0; i<N1; ++i)
138 	  bd.pt(f(u0 + (N1-i)*du, v0, w0 + N2*dw));
139 
140 	// left edge (downward)
141 	for (unsigned int j=0; j<N2; ++j)
142 	  bd.pt(f(u0, v0, w0 + (N2-j)*dw));
143 
144 	// use corners to approximate distance to camera
145 	pt2 = f(u0 + N1*du, v0, w0);
146 	pt3 = f(u0 + N1*du, v0, w0 + N2*dw);
147 	pt4 = f(u0,         v0, w0 + N2*dw);
148       }
149 
150     else if (dw == 0) // (x,y)
151       {
152 	// bottom edge
153 	for (unsigned int i=0; i<N1; ++i)
154 	  bd.pt(f(u0 + i*du, v0, w0));
155 
156 	// right edge
157 	for (unsigned int j=0; j<N2; ++j)
158 	  bd.pt(f(u0 + N1*du, v0 + j*dv, w0));
159 
160 	// top edge (backward)
161 	for (unsigned int i=0; i<N1; ++i)
162 	  bd.pt(f(u0 + (N1-i)*du, v0 + N2*dv, w0));
163 
164 	// left edge (downward)
165 	for (unsigned int j=0; j<N2; ++j)
166 	  bd.pt(f(u0, v0 + (N2-j)*dv, w0));
167 
168 	// use corners to approximate distance to camera
169 	pt2 = f(u0 + N1*du, v0,         w0);
170 	pt3 = f(u0 + N1*du, v0 + N2*dv, w0);
171 	pt4 = f(u0,         v0 + N2*dv, w0);
172       }
173 
174     else
175       epix_error("Bad call to facet constructor"); // exits
176 
177     bd.close().fill(!m_tint.is_unset());
178 
179     center = 0.25*(pt1 + pt2 + pt3 + pt4);
180 
181     direction = center-cam().viewpt();
182     distance = norm(direction);
183 
184     perp = ((pt1 - center)*(pt2 - center));
185 
186     if (norm(perp) < EPIX_EPSILON)
187       perp = (pt3 - center)*(pt4 - center);
188 
189     perp *= recip(norm(perp));
190 
191   } // end of facet(f(double, double, double), ...)
192 
193 
194   // for surface of rotation
facet(double f (double),double g (double),double u0,double v0,double du,double dv,const unsigned int N1,const unsigned int N2,const frame & coords)195   facet::facet(double f(double), double g(double),
196 	       double u0, double v0,
197 	       double du, double dv,
198 	       const unsigned int N1, const unsigned int N2,
199 	       const frame& coords)
200     : m_tint(the_paint_style().fill_color()),
201       m_line(the_paint_style().line_pen()),
202       m_fill(the_paint_style().fill_flag()),
203       pt1(f(u0)*coords.sea() +
204 	  g(u0)*Cos(v0)*coords.sky() +
205 	  g(u0)*Sin(v0)*coords.eye()),
206       pt2(f(u0 + N1*du)*coords.sea() +
207 	  g(u0 + N1*du)*Cos(v0)*coords.sky() +
208 	  g(u0 + N1*du)*Sin(v0)*coords.eye()),
209       pt3(f(u0 + N1*du)*coords.sea() +
210 	  g(u0 + N1*du)*Cos(v0 + N2*dv)*coords.sky() +
211 	  g(u0 + N1*du)*Sin(v0 + N2*dv)*coords.eye()),
212       pt4(f(u0)*coords.sea() +
213 	  g(u0)*Cos(v0 + N2*dv)*coords.sky() +
214 	  g(u0)*Sin(v0 + N2*dv)*coords.eye()),
215       center(0.25*(pt1 + pt2 + pt3 + pt4)),
216       direction(center-cam().viewpt()), distance(norm(direction))
217   {
218     perp = ((pt1 - center)*(pt2 - center));
219 
220     if (norm(perp) < EPIX_EPSILON)
221       perp = (pt3 - center)*(pt4 - center);
222 
223     perp *= recip(norm(perp));
224 
225     // bottom edge
226     for (unsigned int i=0; i<N1; ++i)
227       bd.pt(f(u0 + i*du)*coords.sea() +
228 	    g(u0 + i*du)*Cos(v0)*coords.sky() +
229 	    g(u0 + i*du)*Sin(v0)*coords.eye());
230 
231     // right edge
232     for (unsigned int j=0; j<N2; ++j)
233       bd.pt(f(u0 + N1*du)*coords.sea() +
234 	    g(u0 + N1*du)*Cos(v0 + j*dv)*coords.sky() +
235 	    g(u0 + N1*du)*Sin(v0 + j*dv)*coords.eye());
236 
237     // top edge (backward)
238     for (unsigned int i=0; i<N1; ++i)
239       bd.pt(f(u0 + (N1-i)*du)*coords.sea() +
240 	    g(u0 + (N1-i)*du)*Cos(v0 + N2*dv)*coords.sky() +
241 	    g(u0 + (N1-i)*du)*Sin(v0 + N2*dv)*coords.eye());
242 
243     // left edge (downward)
244     for (unsigned int j=0; j<N2; ++j)
245       bd.pt(f(u0)*coords.sea() +
246 	    g(u0)*Cos(v0 + (N2-j)*dv)*coords.sky() +
247 	    g(u0)*Sin(v0 + (N2-j)*dv)*coords.eye());
248 
249     bd.close().fill(!m_tint.is_unset());
250   }
251 
252   //// Color-dependent constructors
facet(P f (double,double),double u0,double v0,double du,double dv,const unsigned int N1,const unsigned int N2,const Color & tint)253   facet::facet(P f(double, double),
254 	       double u0, double v0,
255 	       double du, double dv,
256 	       const unsigned int N1, const unsigned int N2,
257 	       const Color& tint)
258     : m_tint(tint), m_line(the_paint_style().line_pen()),
259       m_fill(the_paint_style().fill_flag()),
260       pt1(f(u0,         v0)),
261       pt2(f(u0 + N1*du, v0)),
262       pt3(f(u0 + N1*du, v0 + N2*dv)),
263       pt4(f(u0,         v0 + N2*dv)),
264       center(0.25*(pt1 + pt2 + pt3 + pt4)),
265       direction(center-cam().viewpt()),
266       distance(norm(direction))
267   {
268     perp = ((pt1 - center)*(pt2 - center));
269 
270     if (norm(perp) < EPIX_EPSILON)
271       perp = (pt3 - center)*(pt4 - center);
272 
273     perp *= recip(norm(perp));
274 
275     // bottom edge
276     for (unsigned int i=0; i<N1; ++i)
277       bd.pt(f(u0 + i*du, v0));
278 
279     // right edge
280     for (unsigned int j=0; j<N2; ++j)
281       bd.pt(f(u0 + N1*du, v0 + j*dv));
282 
283     // top edge (backward)
284     for (unsigned int i=0; i<N1; ++i)
285       bd.pt(f(u0 + (N1-i)*du, v0 + N2*dv));
286 
287     // left edge (downward)
288     for (unsigned int j=0; j<N2; ++j)
289       bd.pt(f(u0, v0 + (N2-j)*dv));
290 
291     bd.close().fill(!m_tint.is_unset());
292   }
293 
294 
295   // facet constructor for f(double, double, double)
facet(P f (double,double,double),double u0,double v0,double w0,double du,double dv,double dw,const unsigned int N1,const unsigned int N2,const Color & tint)296   facet::facet(P f(double, double, double),
297 	       double u0, double v0, double w0,
298 	       double du, double dv, double dw,
299 	       const unsigned int N1, const unsigned int N2,
300 	       const Color& tint)
301     : m_tint(tint), m_line(the_paint_style().line_pen()),
302       m_fill(the_paint_style().fill_flag()),
303       pt1(f(u0, v0, w0))
304   {
305     if (du == 0) // (y,z)
306       {
307 	// bottom edge
308 	for (unsigned int i=0; i<N1; ++i)
309 	  bd.pt(f(u0, v0 + i*dv, w0));
310 
311 	// right edge
312 	for (unsigned int j=0; j<N2; ++j)
313 	  bd.pt(f(u0, v0 + N1*dv, w0 + j*dw));
314 
315 	// top edge (backward)
316 	for (unsigned int i=0; i<N1; ++i)
317 	  bd.pt(f(u0, v0 + (N1-i)*dv, w0 + N2*dw));
318 
319 	// left edge (downward)
320 	for (unsigned int j=0; j<N2; ++j)
321 	  bd.pt(f(u0, v0, w0 + (N2-j)*dw));
322 
323 	// use corners to approximate distance to camera
324 	pt2 = f(u0, v0 + N1*dv, w0);
325 	pt3 = f(u0, v0 + N1*dv, w0 + N2*dw);
326 	pt4 = f(u0, v0,         w0 + N2*dw);
327       }
328 
329     else if (dv == 0) // (x,z)
330       {
331 	// bottom edge
332 	for (unsigned int i=0; i<N1; ++i)
333 	  bd.pt(f(u0 + i*du, v0, w0));
334 
335 	// right edge
336 	for (unsigned int j=0; j<N2; ++j)
337 	  bd.pt(f(u0 + N1*du, v0, w0 + j*dw));
338 
339 	// top edge (backward)
340 	for (unsigned int i=0; i<N1; ++i)
341 	  bd.pt(f(u0 + (N1-i)*du, v0, w0 + N2*dw));
342 
343 	// left edge (downward)
344 	for (unsigned int j=0; j<N2; ++j)
345 	  bd.pt(f(u0, v0, w0 + (N2-j)*dw));
346 
347 	// use corners to approximate distance to camera
348 	pt2 = f(u0 + N1*du, v0, w0);
349 	pt3 = f(u0 + N1*du, v0, w0 + N2*dw);
350 	pt4 = f(u0,         v0, w0 + N2*dw);
351       }
352 
353     else if (dw == 0) // (x,y)
354       {
355 	// bottom edge
356 	for (unsigned int i=0; i<N1; ++i)
357 	  bd.pt(f(u0 + i*du, v0, w0));
358 
359 	// right edge
360 	for (unsigned int j=0; j<N2; ++j)
361 	  bd.pt(f(u0 + N1*du, v0 + j*dv, w0));
362 
363 	// top edge (backward)
364 	for (unsigned int i=0; i<N1; ++i)
365 	  bd.pt(f(u0 + (N1-i)*du, v0 + N2*dv, w0));
366 
367 	// left edge (downward)
368 	for (unsigned int j=0; j<N2; ++j)
369 	  bd.pt(f(u0, v0 + (N2-j)*dv, w0));
370 
371 	// use corners to approximate distance to camera
372 	pt2 = f(u0 + N1*du, v0,         w0);
373 	pt3 = f(u0 + N1*du, v0 + N2*dv, w0);
374 	pt4 = f(u0,         v0 + N2*dv, w0);
375       }
376 
377     else
378       epix_error("Bad call to facet constructor"); // exits
379 
380     bd.close().fill(!m_tint.is_unset());
381 
382     center = 0.25*(pt1 + pt2 + pt3 + pt4);
383 
384     direction = center-cam().viewpt();
385     distance = norm(direction);
386 
387     perp = ((pt1 - center)*(pt2 - center));
388 
389     if (norm(perp) < EPIX_EPSILON)
390       perp = (pt3 - center)*(pt4 - center);
391 
392     perp *= recip(norm(perp));
393 
394   } // end of facet(f(double, double, double), ...)
395 
396 
397   // for surface of rotation
facet(double f (double),double g (double),double u0,double v0,double du,double dv,const unsigned int N1,const unsigned int N2,const Color & tint,const frame & coords)398   facet::facet(double f(double), double g(double),
399 	       double u0, double v0,
400 	       double du, double dv,
401 	       const unsigned int N1, const unsigned int N2,
402 	       const Color& tint, const frame& coords)
403     : m_tint(tint), m_line(the_paint_style().line_pen()),
404       m_fill(the_paint_style().fill_flag()),
405       pt1(f(u0)*coords.sea() +
406 	  g(u0)*Cos(v0)*coords.sky() +
407 	  g(u0)*Sin(v0)*coords.eye()),
408       pt2(f(u0 + N1*du)*coords.sea() +
409 	  g(u0 + N1*du)*Cos(v0)*coords.sky() +
410 	  g(u0 + N1*du)*Sin(v0)*coords.eye()),
411       pt3(f(u0 + N1*du)*coords.sea() +
412 	  g(u0 + N1*du)*Cos(v0 + N2*dv)*coords.sky() +
413 	  g(u0 + N1*du)*Sin(v0 + N2*dv)*coords.eye()),
414       pt4(f(u0)*coords.sea() +
415 	  g(u0)*Cos(v0 + N2*dv)*coords.sky() +
416 	  g(u0)*Sin(v0 + N2*dv)*coords.eye()),
417       center(0.25*(pt1 + pt2 + pt3 + pt4)),
418       direction(center-cam().viewpt()), distance(norm(direction))
419   {
420     perp = ((pt1 - center)*(pt2 - center));
421 
422     if (norm(perp) < EPIX_EPSILON)
423       perp = (pt3 - center)*(pt4 - center);
424 
425     perp *= recip(norm(perp));
426 
427     // bottom edge
428     for (unsigned int i=0; i<N1; ++i)
429       bd.pt(f(u0 + i*du)*coords.sea() +
430 	    g(u0 + i*du)*Cos(v0)*coords.sky() +
431 	    g(u0 + i*du)*Sin(v0)*coords.eye());
432 
433     // right edge
434     for (unsigned int j=0; j<N2; ++j)
435       bd.pt(f(u0 + N1*du)*coords.sea() +
436 	    g(u0 + N1*du)*Cos(v0 + j*dv)*coords.sky() +
437 	    g(u0 + N1*du)*Sin(v0 + j*dv)*coords.eye());
438 
439     // top edge (backward)
440     for (unsigned int i=0; i<N1; ++i)
441       bd.pt(f(u0 + (N1-i)*du)*coords.sea() +
442 	    g(u0 + (N1-i)*du)*Cos(v0 + N2*dv)*coords.sky() +
443 	    g(u0 + (N1-i)*du)*Sin(v0 + N2*dv)*coords.eye());
444 
445     // left edge (downward)
446     for (unsigned int j=0; j<N2; ++j)
447       bd.pt(f(u0)*coords.sea() +
448 	    g(u0)*Cos(v0 + (N2-j)*dv)*coords.sky() +
449 	    g(u0)*Sin(v0 + (N2-j)*dv)*coords.eye());
450 
451     bd.close().fill(!m_tint.is_unset());
452   }
453 
454 
clone() const455   facet* facet::clone() const
456   {
457     return new facet(*this);
458   }
459 
how_far() const460   double facet::how_far() const { return distance; }
461 
front_facing() const462   bool facet::front_facing() const
463   {
464     return (-direction|perp) > -EPIX_EPSILON;
465   }
466 
467   // N.B. We assume the fill state is stored by our caller
draw(int cull) const468   void facet::draw(int cull) const
469   {
470     if (( cull ==  1 && front_facing() ) || ( cull == -1 && !front_facing() ))
471       return;
472 
473     // else
474     Color paint(m_tint);
475     if (paint.is_unset())
476       paint = White();
477 
478     Color ink(m_line.color());
479     if (ink.is_unset())
480       ink = paint;
481 
482     if (m_fill)
483       {
484 	// calculate cosine^2 of normal angle
485 	// Magic formula (simulated ambient lighting)
486 	const double dens(0.5*(1+pow(perp|(recip(distance)*direction), 2)));
487 
488 	paint *= dens;
489 	ink   *= dens;
490       }
491 
492     bd.draw(paint, pen_data(ink, m_line.width()));
493   } // end of facet::draw(bool, int)
494 
495 
operator ()(const facet & arg1,const facet & arg2)496   bool by_distance::operator() (const facet& arg1, const facet& arg2)
497   {
498     return arg1.how_far() > arg2.how_far();
499   }
500 
operator ()(const facet * arg1,const facet * arg2)501   bool by_distance::operator() (const facet* arg1, const facet* arg2)
502   {
503     return arg1->how_far() > arg2->how_far();
504   }
505 } // end of namespace
506