home *** CD-ROM | disk | FTP | other *** search
/ Amiga ISO Collection / AmigaUtilCD2.iso / Programming / GCC / GERLIB_DEV08B.LHA / gerlib / libg++ / etc / PlotFile3D / PlotFile3D.cc < prev    next >
Encoding:
C/C++ Source or Header  |  1993-12-12  |  14.3 KB  |  502 lines
  1. // Source file for PlotFile3D class -*- C++ -*-
  2. /* 
  3. Copyright (C) 1990 Free Software Foundation
  4.     written by J. Thomas Ngo, Harvard University
  5.  
  6. This file is part of the GNU C++ Library.
  7.  
  8. The GNU C++ Library is distributed in the hope that it will be useful,
  9. but WITHOUT ANY WARRANTY.  No author or distributor accepts
  10. responsibility to anyone for the consequences of using it or for
  11. whether it serves any particular purpose or works at all, unless he
  12. says so in writing.  Refer to the GNU General Public License for full
  13. details.
  14.  
  15. Everyone is granted permission to copy, modify and redistribute The
  16. GNU C++ Library, but only under the conditions described in the GNU
  17. General Public License.  A copy of this license is supposed to have
  18. been given to you along with The GNU C++ Library so you can know your
  19. rights and responsibilities.  It should be in a file named COPYING.
  20. Among other things, the copyright notice and this notice must be
  21. preserved on all copies.  
  22. */
  23.  
  24. // An analog of the PlotFile class, but written for 3D graphics.  In
  25. // the initializer you specify a camera location, and decide whether
  26. // you want stereo viewing.  Most other methods are analogous to the
  27. // unix plot(3) commands, except that they take 3D positions.
  28.  
  29. #ifdef __GNUG__
  30. #pragma implementation
  31. #endif
  32. #include "PlotFile3D.h"
  33.  
  34. // Error handling
  35. void PlotFile3D::error(const char* s) const 
  36. {
  37.   (*lib_error_handler)("PlotFile3D",s);
  38. }
  39.  
  40. //===========================================================================
  41. // Constructor
  42. //===========================================================================
  43.  
  44. static void _find_basis(Vec3D& u, Vec3D& v, Vec3D& w,
  45.                         const double th, const double ph)
  46. {
  47.     double sth = sin(th), cth = cos(th);
  48.     double sph = sin(ph), cph = cos(ph);
  49.     u = Vec3D( sph, cph, 0 );
  50.     v = Vec3D( -sth * cph, sth * sph, cth );
  51.     w = u ^ v;
  52. }
  53.  
  54. void PlotFile3D::initialize(bool st, double th, double ph)
  55. {
  56.     stereo = st;
  57.     ppq = 0;
  58.     valid3D = FALSE;
  59.     if( stereo ) {
  60.         const double stereo_separation = 2.5 * M_PI / 180.0;
  61.         _find_basis( ul,vl,wl, th, ph-stereo_separation );
  62.         _find_basis( ur,vr,wr, th, ph+stereo_separation );
  63.     } else {
  64.         _find_basis( ul,vl,wl, th, ph                   );
  65.         ur = ul; vr = vl; wr = wl;
  66.     }
  67.     undefine3D();
  68. }
  69.  
  70. // destructor
  71.  
  72. PlotFile3D::~PlotFile3D() 
  73. {
  74.   if (ppq != 0) { delete sintab; delete costab; }
  75. }
  76.  
  77. //===========================================================================
  78. // Set plot limits via space()
  79. //===========================================================================
  80.  
  81. PlotFile3D& PlotFile3D::space(const Vec3D& pt0, const Vec3D& pt1)
  82. {
  83.     // Set up 2D plot package
  84.     PlotFile::space(umini,vmini,umaxi,vmaxi);
  85.     
  86.     // Find 2D rectangle bounding the box specified by user
  87.  
  88. #ifdef __GNUC__
  89.     double umin  = (pt0.x() * ul.x()) <? (pt0.x() * ur.x())
  90.                 <? (pt1.x() * ul.x()) <? (pt1.x() * ur.x());
  91.     double umax  = (pt0.x() * ul.x()) >? (pt0.x() * ur.x())
  92.                 >? (pt1.x() * ul.x()) >? (pt1.x() * ur.x());
  93.     double vmin  = (pt0.x() * vl.x()) <? (pt0.x() * vr.x())
  94.                 <? (pt1.x() * vl.x()) <? (pt1.x() * vr.x());
  95.     double vmax  = (pt0.x() * vl.x()) >? (pt0.x() * vr.x())
  96.                 >? (pt1.x() * vl.x()) >? (pt1.x() * vr.x());
  97.  
  98.     umin        += (pt0.y() * ul.y()) <? (pt0.y() * ur.y())
  99.                 <? (pt1.y() * ul.y()) <? (pt1.y() * ur.y());
  100.     umax        += (pt0.y() * ul.y()) >? (pt0.y() * ur.y())
  101.                 >? (pt1.y() * ul.y()) >? (pt1.y() * ur.y());
  102.     vmin        += (pt0.y() * vl.y()) <? (pt0.y() * vr.y())
  103.                 <? (pt1.y() * vl.y()) <? (pt1.y() * vr.y());
  104.     vmax        += (pt0.y() * vl.y()) >? (pt0.y() * vr.y())
  105.                 >? (pt1.y() * vl.y()) >? (pt1.y() * vr.y());
  106.  
  107.     umin        += (pt0.z() * ul.z()) <? (pt0.z() * ur.z())
  108.                 <? (pt1.z() * ul.z()) <? (pt1.z() * ur.z());
  109.     umax        += (pt0.z() * ul.z()) >? (pt0.z() * ur.z())
  110.                 >? (pt1.z() * ul.z()) >? (pt1.z() * ur.z());
  111.     vmin        += (pt0.z() * vl.z()) <? (pt0.z() * vr.z())
  112.                 <? (pt1.z() * vl.z()) <? (pt1.z() * vr.z());
  113.     vmax        += (pt0.z() * vl.z()) >? (pt0.z() * vr.z())
  114.                 >? (pt1.z() * vl.z()) >? (pt1.z() * vr.z());
  115. #else
  116.     double tmin, tmax, tmp;
  117.  
  118.     tmin = tmax = pt0.x() * ul.x();
  119.     tmp = pt0.x() * ur.x();
  120.     if (tmp < tmin) tmin = tmp;
  121.     if (tmp > tmax) tmax = tmp;
  122.     tmp = pt1.x() * ul.x();
  123.     if (tmp < tmin) tmin = tmp;
  124.     if (tmp > tmax) tmax = tmp;
  125.     tmp = pt1.x() * ur.x();
  126.     if (tmp < tmin) tmin = tmp;
  127.     if (tmp > tmax) tmax = tmp;
  128.  
  129.     double umin  = tmin;
  130.     double umax  = tmax;
  131.  
  132.     tmin = tmax = pt0.x() * vl.x();
  133.     tmp = pt0.x() * vr.x();
  134.     if (tmp < tmin) tmin = tmp;
  135.     if (tmp > tmax) tmax = tmp;
  136.     tmp = pt1.x() * vl.x();
  137.     if (tmp < tmin) tmin = tmp;
  138.     if (tmp > tmax) tmax = tmp;
  139.     tmp = pt1.x() * vr.x();
  140.     if (tmp < tmin) tmin = tmp;
  141.     if (tmp > tmax) tmax = tmp;
  142.  
  143.     double vmin  = tmin;
  144.     double vmax  = tmax;
  145.  
  146.     tmin = tmax = pt0.y() * ul.y();
  147.     tmp = pt0.y() * ur.y();
  148.     if (tmp < tmin) tmin = tmp;
  149.     if (tmp > tmax) tmax = tmp;
  150.     tmp = pt1.y() * ul.y();
  151.     if (tmp < tmin) tmin = tmp;
  152.     if (tmp > tmax) tmax = tmp;
  153.     tmp = pt1.y() * ur.y();
  154.     if (tmp < tmin) tmin = tmp;
  155.     if (tmp > tmax) tmax = tmp;
  156.  
  157.     umin        += tmin;
  158.     umax        += tmax;
  159.  
  160.     tmin = tmax = pt0.y() * vl.y();
  161.     tmp = pt0.y() * vr.y();
  162.     if (tmp < tmin) tmin = tmp;
  163.     if (tmp > tmax) tmax = tmp;
  164.     tmp = pt1.y() * vl.y();
  165.     if (tmp < tmin) tmin = tmp;
  166.     if (tmp > tmax) tmax = tmp;
  167.     tmp = pt1.y() * vr.y();
  168.     if (tmp < tmin) tmin = tmp;
  169.     if (tmp > tmax) tmax = tmp;
  170.  
  171.     vmin        += tmin;
  172.     vmax        += tmax;
  173.  
  174.     tmin = tmax = pt0.z() * ul.z();
  175.     tmp = pt0.z() * ur.z();
  176.     if (tmp < tmin) tmin = tmp;
  177.     if (tmp > tmax) tmax = tmp;
  178.     tmp = pt1.z() * ul.z();
  179.     if (tmp < tmin) tmin = tmp;
  180.     if (tmp > tmax) tmax = tmp;
  181.     tmp = pt1.z() * ur.z();
  182.     if (tmp < tmin) tmin = tmp;
  183.     if (tmp > tmax) tmax = tmp;
  184.  
  185.     umin        += tmin;
  186.     umax        += tmax;
  187.  
  188.     tmin = tmax = pt0.z() * vl.z();
  189.     tmp = pt0.z() * vr.z();
  190.     if (tmp < tmin) tmin = tmp;
  191.     if (tmp > tmax) tmax = tmp;
  192.     tmp = pt1.z() * vl.z();
  193.     if (tmp < tmin) tmin = tmp;
  194.     if (tmp > tmax) tmax = tmp;
  195.     tmp = pt1.z() * vr.z();
  196.     if (tmp < tmin) tmin = tmp;
  197.     if (tmp > tmax) tmax = tmp;
  198.  
  199.     vmin        += tmin;
  200.     vmax        += tmax;
  201. #endif
  202.    
  203.     // Map this rectangle on to umini et al. using as much space as
  204.     // possible and centering.  Also, if stereo was requested, split
  205.     // the window right down the middle and try to fit the picture in.
  206.     int pixwidth  = umaxi - umini;
  207.     int pixheight = vmaxi - vmini;
  208.     if (stereo) pixwidth /= 2;
  209.     double uscale = pixwidth  / (umax-umin);
  210.     double vscale = pixheight / (vmax-vmin);
  211.     if (uscale <= vscale) scale = uscale;
  212.     else { scale = vscale; pixwidth = int(pixwidth * vscale/uscale); }
  213.     scale *= margin;
  214.     uorig = 0.5 * (pixwidth  - scale*(umin+umax));
  215.     vorig = 0.5 * (pixheight - scale*(vmin+vmax));
  216.     if (stereo) stereo_offset = pixwidth;
  217.     
  218.     // Set state appropriately
  219.     return home();
  220. }
  221.  
  222.  
  223. //===========================================================================
  224. // move() and cont() are used as primitives by the other routines
  225. //===========================================================================
  226.  
  227. PlotFile3D& PlotFile3D::move(const Vec3D& p)
  228. {
  229.     intp = p;
  230.     define3D();
  231.     return *this;
  232. }
  233.  
  234.  
  235. void PlotFile3D::project(int& ui,int& vi, const Vec3D& u, const Vec3D& v,
  236.                           const Vec3D& rel) const
  237. {
  238.   ui = int(rel * u * scale + uorig);
  239.   vi = int(rel * v * scale + vorig);
  240. }
  241.  
  242. static
  243. #ifdef __GNUC__
  244. inline
  245. #endif
  246. bool _clip(int& a, const int amin, const int amax,
  247.                          int& b, const int afar, const int bfar)
  248. {
  249.     double frac = 0;
  250.     if( a < amin ) {
  251.         if( afar <= amin ) return FALSE;
  252.         frac = ((double)amin-a)/(afar-a);
  253.         a = amin;
  254.     } else if( a > amax ) {
  255.         if( afar >= amax ) return FALSE;
  256.         frac = ((double)amax-a)/(afar-a);
  257.         a = amax;
  258.     }
  259.     if( frac != 0 ) b += int(frac*(bfar-b));
  260.     return TRUE;
  261. }
  262.  
  263. void PlotFile3D::line2D(int u0, int v0, int u1, int v1)
  264. {
  265.     // Clip against each edge in turn
  266.     if( !_clip(u0,umini,umaxi,v0,u1,v1) ) return;
  267.     if( !_clip(v0,vmini,vmaxi,u0,v1,u1) ) return;
  268.     if( !_clip(u1,umini,umaxi,v1,u0,v0) ) return;
  269.     if( !_clip(v1,vmini,vmaxi,u1,v0,u0) ) return;
  270.  
  271.     // Draw clipped segment
  272.     PlotFile::move(u0,v0); PlotFile::cont(u1,v1);
  273. }
  274.  
  275.  
  276. PlotFile3D& PlotFile3D::cont(const Vec3D& p)
  277. {
  278.     must_be_valid3D();
  279.     int u0, v0, u1, v1;
  280.     if (stereo) {
  281.         project(u0,v0, ul,vl, intp);
  282.         project(u1,v1, ul,vl, p);
  283.         line2D(u0,v0,u1,v1);
  284.         project(u0,v0, ur,vr, intp); u0 += stereo_offset;
  285.         project(u1,v1, ur,vr, p);    u1 += stereo_offset;
  286.         line2D(u0,v0,u1,v1);
  287.     } else {
  288.         project(u0,v0, ul,vl, intp);
  289.         project(u1,v1, ul,vl, p);
  290.         line2D(u0,v0,u1,v1);
  291.     }
  292.     intp = p;
  293.     return *this;
  294. }
  295.  
  296.  
  297. //===========================================================================
  298. // These last commands are generated by sequences of move() and cont()
  299. //===========================================================================
  300.  
  301.  
  302. PlotFile3D& PlotFile3D::box(const Vec3D& p0, const Vec3D& p1)
  303. {
  304.   // The corners:
  305.   Vec3D x0y0z1 (p0.x(),p0.y(),p1.z()); 
  306.   Vec3D x0y1z0 (p0.x(),p1.y(),p0.z()); 
  307.   Vec3D x0y1z1 (p0.x(),p1.y(),p1.z()); 
  308.   Vec3D x1y0z0 (p1.x(),p0.y(),p0.z()); 
  309.   Vec3D x1y0z1 (p1.x(),p0.y(),p1.z());
  310.   Vec3D x1y1z0 (p1.x(),p1.y(),p0.z()); 
  311.  
  312.   // Bottom rectangle
  313.   move(p0);
  314.   cont(x1y0z0);
  315.   cont(x1y1z0);
  316.   cont(x0y1z0);
  317.   cont(p0);
  318.  
  319.   // One leg
  320.   cont(x0y0z1);
  321.   
  322.   // Top rectangle
  323.   cont(x1y0z1);
  324.   cont(p1);
  325.   cont(x0y1z1);
  326.   cont(x0y0z1);
  327.  
  328.   // Three other legs
  329.   line(x1y0z0, x1y0z1);
  330.   line(x1y1z0, p1);
  331.   line(x0y1z0, x0y1z1);
  332.  
  333.   // Set state
  334.   undefine3D(); 
  335.   return *this;
  336. }
  337.  
  338.  
  339. PlotFile3D& PlotFile3D::circle(const Vec3D& center, const Vec3D& rad,
  340.                                const int points_per_quadrant )
  341. {
  342.   int i;
  343.   
  344.   // Build trig lookup table if necessary
  345.  
  346.   if( ppq != points_per_quadrant ) 
  347.   {
  348.     if( ppq != 0 ) { delete sintab; delete costab; }
  349.     ppq = points_per_quadrant;
  350.     sintab = new double[ppq]; costab = new double[ppq];
  351.     for( i=0; i<ppq; i++ ) sintab[i] = sin(double(i)*M_PI_2/double(ppq));
  352.     for( i=0; i<ppq; i++ ) costab[i] = cos(double(i)*M_PI_2/double(ppq));
  353.   }
  354.  
  355.   // Find radius vector normal to circle, and two vectors
  356.   // orthogonal to it
  357.  
  358.   Vec3D nz = rad;
  359.  
  360.   double radius = mod(nz); 
  361.   nz /= radius;
  362.  
  363.   double tx, ty, tz;
  364.  
  365.   if( nz.x() <= nz.y() && nz.x() <= nz.z() ) 
  366.   {
  367.     tx = 0; ty = nz.z(); tz = -nz.y();
  368.   }
  369.   else if ( nz.y() <= nz.z() )
  370.   {
  371.     tx = nz.z(); ty = 0; tz = -nz.x();
  372.   }
  373.   else
  374.   {
  375.     tx = nz.y(); ty = -nz.x(); tz = 0;
  376.   }
  377.   
  378.   Vec3D nx (tx, ty, tz);
  379.  
  380.   nx /= mod(nx);      // Unit vector orthogonal to nz
  381.  
  382.   Vec3D ny = nz ^ nx;       // Unit vector orthogonal to ny and nz
  383.  
  384.   nx *= radius; ny *= radius; // Now they're basis vectors for circle
  385.  
  386.  
  387.   // Iterate...
  388.   Vec3D cnx(center+nx);
  389.   move(cnx);
  390.  
  391.   for( i=0; i<ppq; i++ ) cont(center + nx * costab[i] + ny * sintab[i]);
  392.   for( i=0; i<ppq; i++ ) cont(center + ny * costab[i] - nx * sintab[i]);
  393.   for( i=0; i<ppq; i++ ) cont(center - nx * costab[i] - ny * sintab[i]); 
  394.   for( i=0; i<ppq; i++ ) cont(center - ny * costab[i] + nx * sintab[i]); 
  395.  
  396.   cont(cnx);
  397.  
  398.   // Set state
  399.   undefine3D(); 
  400.   return *this;
  401. }
  402.  
  403.  
  404. //===========================================================================
  405. // Methods not found in plot(3)--particular to this package
  406. //===========================================================================
  407.  
  408. PlotFile3D& PlotFile3D::home()  // Cursor to upper left corner
  409. {
  410.     int homex = umini;
  411.     int homey = int(vmaxi + 0.05 * (vmini-vmaxi));
  412.     PlotFile::move(homex,homey);
  413.     undefine3D();
  414.     return *this;
  415. }
  416.  
  417.  
  418. PlotFile3D& PlotFile3D::sphere(const Vec3D& center, const Vec3D& rad,
  419.                                const int points_per_quadrant)
  420. {
  421.   Vec3D eyeline = wl * mod(rad);
  422.   circle( center, eyeline, points_per_quadrant );
  423.   return circle(center, rad, points_per_quadrant );
  424. }
  425.  
  426.  
  427. //===========================================================================
  428. // Raw coordinate versions of things
  429. //===========================================================================
  430.  
  431. PlotFile3D& PlotFile3D::space(const double x0, const double y0, 
  432.                               const double z0,
  433.                               const double x1, const double y1, 
  434.                               const double z1)
  435. {
  436.   Vec3D p0(x0, y0, z0);
  437.   Vec3D p1(x1, y1, z1);
  438.   return space(p0, p1);
  439. };
  440.  
  441. PlotFile3D& PlotFile3D::move(const double xi, const double yi, const double zi)
  442. {
  443.     intp = Vec3D(xi, yi, zi);
  444.     define3D();
  445.     return *this;
  446. }
  447.  
  448. PlotFile3D& PlotFile3D::cont(const double xi, const double yi, const double zi)
  449. {
  450.   Vec3D p(xi, yi, zi);
  451.   return cont(p);
  452. }
  453.  
  454. PlotFile3D&  PlotFile3D::box(const double x0, const double y0, const double z0,
  455.                              const double x1, const double y1, const double z1)
  456. {
  457.   Vec3D x0y0z0(x0, y0, z0);
  458.   Vec3D x1y1z1(x1, y1, z1);
  459.   return box(x0y0z0, x1y1z1);
  460. }
  461.  
  462. PlotFile3D& PlotFile3D::circle(const double cx,const double cy,const double cz,
  463.                    const double rx,const double ry,const double rz,
  464.                    const int points_per_quadrant )
  465.     // cx,cy,cz:  position of circle's center
  466.     // rx,ry,rz:  length is circle's radius; direction is normal to circle
  467. {
  468.   Vec3D center(cx, cy, cz);
  469.   Vec3D rad(rx, ry, rz);
  470.   return circle(center, rad, points_per_quadrant);
  471. }
  472.  
  473. PlotFile3D& PlotFile3D::sphere(const double cx,const double cy,const double cz,
  474.                    const double rx,const double ry,const double rz,
  475.                    const int points_per_quadrant)
  476. {
  477.   Vec3D center(cx, cy, cz);
  478.   Vec3D rad(rx, ry, rz);
  479.   return sphere(center, rad, points_per_quadrant);
  480. }
  481.  
  482. //===========================================================================
  483. // Test program
  484. //===========================================================================
  485.  
  486. #ifdef TEST
  487. #include <stream.h>
  488. main()
  489. {
  490.   PlotFile3D foo(fopen("test.pl", "w"));
  491.   foo.space(-1000,-1000,-1000,1000,1000,1000).erase();
  492.       
  493.   // Draw and label
  494.   foo.linemod("longdashed").box(-1000,-1000,-1500,1000,1000,0);
  495.   foo.linemod("solid").sphere(0,0,500,0,0,500,30);
  496.   foo.linemod("dotted").circle(0,0,500,0,0,800,50);
  497.   char title[] = "Test of PlotFile3D package";
  498.   foo.home().label(title);
  499.   
  500. }
  501. #endif
  502.