Celestin Apprentice 7/Apprentice-Release7.iso/Source Code/C ++/Applications/Venus 3.5/view_3d.cc

Text File | 1997-04-23 | 29.9 KB | 726 lines | [TEXT/CWIE]

/*
 ************************************************************************
 *
 *                   2D rendering of a surface in 3D space (a map)
 *
 * Suppose we have a surface in space z=f(x,y) and we want to render it to view
 * on a screen (in 2D). The function z=f(x,y) may be specified as a 2D matrix, or
 * as an image/map. In the latter case we also need to know how to transform a pixel
 * value into z, an elevation of the corresponding point (x,y) on the map.
 *
 * Coordinate systems and notation
 * World coordinate system: axes Ox, Oy, Oz. We choose Oz to go straight up (the higher
 * z, the higher the elevation), and lay Ox, Oy along the map's column/row respectively,
 * in the direction of increasing column/row indices. We assume that z=0 corresponds to
 * the base of a map (base of clouds).
 *
 * We set the view point (eye point, camera point) at point E (xe,ye,ze), with the direction
 * of gaze specified by the vector g (gx,gy,gz). That is, an observer looks from point E
 * in the general direction g, and sees the image of the world projected onto a focal
 * plane ƒ. Thus the local coordinate system (associated with the view plane) thus is:
 * vector g is a "depth" vector, vector v of the plane shows the "up" direction,
 * and vector u of the plane gives a horizontal direction; u is chosen to make
 * a left-handed coordinate system, that is, u = v x g). The easiest way to specify
 * the up direction, v, is to project (Oz) onto the plane:
 *        v = (0,0,1) - g((0,0,1),g)/|g|^2
 * where (a,b) is a scalar (dot) product of two vectors, a and b.
 *
 * We need an equation of the focal plane to project the world onto it. The focal plane
 * is orthogonal to the gaze vector g and is located at a distance f0 off the view
 * point E. This gives us
 *                        (g,(r-E)) - f0*|g| = 0
 * Note that point F, a projection of E onto the plane, can be found from
 *                        (g,(F-E)) - f0*|g| = 0
 * well, since EF || g, F = f0*g/|g| + E
 *
 * To project a point A (x,y,z) of the world onto the focal plane, we draw a line
 * from E to A, whose intersection with the focal plane (point B) is the desired
 * projection.
 * Note that EB || EA, that is, EB = p*EA, and B = p*(A-E) + E, where p is some scalar.
 * To find it, we note that B belongs to ƒ, that is, it has to satisfy the focal plane
 * equation
 *            (g,(B-E)) = f0*|g|,
 * which quickly gives us p:
 *            p = f0*|g|/(g,(A-E))
 * To find "local" coordinates of B on a focal plane, we need to project
 * vector FB onto vectors u and v. But first let's determine a "view depth"
 * of A, which is needed for hidden surface/line removal.
 * One can choose the depth of A to be either a distance from A to the focal plane,
 * or the length of vector BA. The latter seems to be more natural,
 *    |BA| = (1-p)|EA|
 * but we'll run into problems later (see below). So we assume that the depth
 * of A, d(A), is its distance from the focal plane
 *    d(A) = (g,(A-E))/|g| - f0
 * For practical purposes, it's better to define the "depth" simply as
 *    d(A) = (g,(A-E))/|g|
 * which could range from f0 (for points on the focal plane) through some
 * upper limit s0 (sight depth).
 * Coming back to vector FB, it can be estimated as follows
 *    FB = B - F = p*(A-E) + E - f0*g/|g| - E = p*(A-E) - f0*g/|g|
 *
 * Note, that g is orthogonal to both u and v vectors of the local coordinate system,
 * so
 *    u = f0/d(A) * (u,(A-E))
 *    v = f0/d(A) * (v,(A-E))
 * To complete the transformation to screen coordinates u and v, we multiply by a
 * scaling factor and add an appropriate offset
 *  u = beta*f0/d(A) * (u,(A-E)) + Ou
 *  v = beta*f0/d(A) * (v,(A-E)) + Ov
 * (note, scaling factors for u and v could be different, it depends on the screen
 * aspect ratio; but we assume it is 1 for now).
 *
 * The drawing method we're going to use, voxel method, involves scanning of the
 * world from fatherst to closest points, and from left to right. That is, given
 * a fixed depth d and azimuth u, we should be able to find a corresponding point
 * (x,y) on the map. It's easy to see it involves solving a set of equations
 *    (x-xe)*gx + (y-ye)*gy + (z-ze)*gz = d*|g|
 *  (x-xe)*ux + (y-ye)*uy + (z-ze)*uz = const
 *  z = z(x,y)
 * The first equation defines a plane (parallel to ƒ ) that intersects the 2d surface
 * (elevation map) along some line(s), which is to be projected to the view plane.
 * Unfortunately, unless the plane is vertical, the lines are many (consider what a
 * horizontal plane would do slicing through the "mountain peaks"). For the voxel algorithm
 * below, we really need only a single nice intersection (without looping, etc).
 * We have to make a simplification: assume that the gaze g is strictly horizontal,
 * and the up direction v is straight "up", that is, v || (Oz). That means that tilting
 * is out, it has to be done by transforming the map itself.
 * Still, we can _simulate_ tilting by making distant peaks appear taller than they
 * would otherwise be in the perspective projection, see below.
 *
 * If the gaze is strictly horizontal, gz=0. The focal plane ƒ is then vertical, and
 * the up direction (vector v) can be chosen to be (0,0,1). The u vector of the trio
 * is then u = (v x g)/|v x g| = (-gy gx 0)/|g|
 * So, the projection formulas from A(x,y,z) into B(u,v,d) become
 *     u = beta*f0/d(A) * (-gy*(x-xe) + gx*(y-ye))/|g| + Ou
 *     v = beta*f0/d(A) * (z-ze) + Ov
 *     d(A) = (gx*(x-xe) + gy*(y-ye))/|g|
 * Without loss of generality, we can assume f0=1 (that is, it can be absorbed in
 * the scaling factor beta, an arbitrary constant at the moment).
 *
 * We scan the world along the planes d=const, and then for each plane of constant
 * depth we scan our 'view port' from left to right, that is, from u=0 to u=D-1,
 * D being the width of our view scope. For each scan point (u,d) we should be able
 * to figure out the corresponding point (x,y) on the map, find out its elevation
 * z=z(x,y), and compute the projection v.
 * First thing first, given d and u, we need to find x and y from
 *    -gy*(x-xe) + gx*(y-ye) = (u-Ou)*d*|g|/beta
 *   gx*(x-xe) + gy*(y-ye) = d*|g|
 * The solution is obviously
 *    x = xe + (d*gx - gy*(u-Ou)*d/beta)/|g|
 *  y = ye + (d*gy + gx*(u-Ou)*d/beta)/|g|
 * The most efficient way of computing x and y is to evaluate the previous formulas
 * at u=0:
 *  xr = xe + (d*gx + gy*Ou*d/beta)/|g|
 *  yr = ye + (d*gy - gx*Ou*d/beta)/|g|
 * and then evaluate increments xinc, yinc when u changes by one
 *  xinc = - gy*d/beta/|g|
 *  yinc =   gx*d/beta/|g|
 * Keeping in mind that d=const at a u scan.
 * It's easy to see that it's convenient to express xinc, yinc in units of |g|, and
 * select |g| to be something nice, like |g|=1<<14. Thus |g| makes a very suitable units
 * in which xr, xinc, etc quantities would be measured. The neatest thing is that
 * all these quantities would be integral in |g| units (or can be rounded to be integral
 * without losing too much precision). Note that |g| is better be that big. At each step
 * in d in the code below, xinc, yinc are decremented by -gy/beta, gx/beta. We need to make
 * at least 100 of these steps for good representation of depth. It follows then that |g|/beta
 * should be of order 10-100 (so max error of 1 in xinc decrement won't translate into
 * a huge accumulated error 100/|g|*number_of_u_steps). Otherwise the animation would be
 * jerky.
 *
 * So, the algorithm to process scanlines d=const is as follows:
 * 1. fix d=sight_depth and compute
 *      xr = xe*|g| + d*gx + gy*Ou*d/beta                // In units of |g|
 *      yr = ye*|g| + d*gy - gx*Ou*d/beta
 *      xinc = - gy*d/beta
 *      yinc =   gx*d/beta
 *      zmax = ze + (vmax-Ov)*d/beta
 *      zmin = ze + (vmin-Ov)*d/beta
 * 3. set xu=xr, yu=yr
 *    for u=0 to D-1
 * 4.     compute x=xu/|g|, y=yu/|g|
 * 5.     find z=f(x,y) from the map
 *          if z > zmax set v=vmax
 *          else if z < zmin set v=vmin
 *          else v = Ov + beta*(z-ze)/d
 *          Note a scanline point (u,v) and connect to the previous scanline
 *          xu += xinc; yu += yinc
 *    end u-loop
 *      d -= 1                                    // move a step closer to the focal plane
 *      xr -= gx + gy*Ou/beta                        // and recompute the "ref" quantities
 *      yr -= gy - gx*Ou/beta
 *      xinc -= - gy/beta
 *      yinc -=   gx/beta
 *      zmax -= (vmax-Ov)/beta-1                // That'll keep [zmin,zmax] a bit wider (by 2)
 *      zmin -= (vmin-Ov)/beta+1                // than necessary to offset possible "round-off" errors
 *  repeat while d > |g|
 *
 * Again, all the quantities above, xr, yr, xinc, etc. are integral. That means
 * that all the divisions are integer divisions, too (and can be precomputed
 * in all the cases but one, evaluation of v).
 *
 * Since we scan from larger d's to smaller d's (that is, from farther ahead
 * to near, that is, towards the observer), then if a (u,v) point generated by
 * the algorithm should obscure (u,v) point generated at an earlier pass
 * (with larger d), it would just overwrite the old point. However, if we just
 * draw (u,v) dots we've computed we get a dotty picture. Note, that actually we see
 * a segment of a line (x,y,z)-(x1,y-dy,z1). If we connect two points generated
 * in two consecutive scans, (u,v1) and (u,vold), it wouldn't be always
 * right, because the original line could been obscured. Suppose that
 * ze (elevation of the observation point) is high enough (comparable with
 * the hight of the tallest peaks). Then if v(u,d+delta_d) > v(u,d) for some
 * u and delta_d>0, then we have a descending slope, and the line is visible.
 * Otherwise the line is on ascending slope, and it couldn't be visible
 * because it will be obscured by the descending slope (which should be
 * closer to us). Note the line from (u,c1) to (u,vold) is a vertical one,
 * and therefore, is very easy to draw (and clip, if necessary).
 *
 * When we create (u,v) map, we assign to the point (u,v) the intensity
 * (color) of the f(x,y) point of the original map. We can use Gouraud
 * interpolation in drawing the line.
 *
 * Colors and elevations: we assume that the pixel value of the map image,
 * map(j,i) is an "elevation code" for the corresponding point x,y. That means,
 * that the true elevation of the point can be computed via some simple formula
 * z = map(j,i) < z_cut_off ? 0 : (map(i,j)*elevation_factor)>>9;
 * This formula can scale the pixel values both up and down.
 * Moreover, each "elevation code" is associated with a corresponding color via a
 * private colormap, see class. We assume that this colormap is arranged in some order
 * to allow "interpolation". That is, the color for the entry k should be in a sense
 * "in-between" the colors for entries k-1 and k+1.
 * Note if we decrease the ze as we project closer and closer point,
 * we can simulate a "camera tilt". it's not going to a be a genuine tilt (say,
 * we can't look into a distant valley this way); still, it seems to create a
 * good impression.
 *
 * The color-interpolation algorithm warrants some explanation. Suppose we are to
 * draw a vertical line from v1 to v2 (v2>v1). Point v1 has color index c1,
 * point v2 has color index c2. Suppose c2 > c1 and c_span = c2-c1 < v_span = v2-v1
 * We will use an algorithm that is in spirit of Bresenham's line-drawing algorithm
 * (and does not use any multiplications/divisions)
 * The algorithm is simple: we start at point v1 going up to v2, assigning each
 * pixel, as we go along, color curr_color, or ++curr_color. At each point we have
 * to choose whether to increment the curr_color or not. Notice that curr_color
 * has to be incremented every round(v_span/c_span) steps. If we make an accumulator
 * and increment it by c_span at every v step, curr_color has to be incremented when the
 * accumulator reaches v_span. Thus the algorithm is simple:
 *  acc = 1            // equivalent to "rounding up" at round(v_span/c_span)
 *  curr_color = c1
 *  for v=v1 to v2
 *    pixel = curr_color
 *    acc += c_span
 *    if( acc >= v_span )
 *      acc -= v_span, ++curr_color
 *  end
 *
 * Note that in accessing pixels of an IMAGE, map(j,i), the first argument is the
 * row index (corresponding to y), and the second argument is the col index
 * (corresponding to x).
 *
 * Inspiration:
 * Tim Clarke, [e-mail]
 * But this version of Venus is *very* far away from MARS.
 *
 * Generalization: draw boxes (quadraluzation), than we can increase the
 * scanning step in u (or in y) in the algorithm above.
 * step on u in two instead of one, and interpolate v in between
 *
 ************************************************************************
 */

#include "ImageViews.h"

#define SIMULATE_TILT 0

                                // Kind of a formal constructor, merely an initialization
ThreeDViewBase::ThreeDViewBase(const ElevationMap& image_map, const ViewerPosition& _viewer_pos,
                       const ProjectionParameters& _projection_parms)
  : OffScreenBuffer(_projection_parms.q_view_rect(),image_map.q_clut()),
    map(image_map),
    viewer_pos(_viewer_pos),
    projection_parms(_projection_parms),
    scanline1(_projection_parms.q_width()),
    scanline2(_projection_parms.q_width())
{
  assert( ViewerPosition::g_unit == (1<<14) );
}

template <class LineProjector>
ThreeDView<LineProjector>::ThreeDView(const ElevationMap& image_map, const ViewerPosition& _viewer_pos,
                       const ProjectionParameters& _projection_parms)
  : ThreeDViewBase(image_map,_viewer_pos,_projection_parms)
{
}

ProjectionParameters::ProjectionParameters
    (const ScreenRect& bounds, const int _ze, const int _beta)
    : ze(_ze), beta(_beta),
      u_sup(bounds.q_width()),
      v_sup(bounds.q_height()),
      Ou(u_sup/2), Ov(v_sup/2)                // Pick up somewhere in the middle
{}

                            // Dump the contents of a ViewerPosition
void ViewerPosition::dump(void) const
{
  message("\nObserver is at (%d,%d) gazing in the direction (%d,%d)\n",
            xe,ye,gx,gy);
}

                            // Dump the contents of ProjectionParameters
void ProjectionParameters::dump(void) const
{
  message("\n-->Projection parameters:");
  message("\n  observers' elevation %d",ze);
  message("\n  zoom factor (local frame scaling factor) %d",beta);
  message("\n  local frame %dx%d",u_sup,v_sup);
  message("\n  local frame origin (%d,%d)\n",Ou,Ov);
}

/*
 *----------------------------------------------------------------------
 *                Dealing with scanlines of the view window buffer
 */

                                // Handling a scanline: allocate one scanline
ThreeDViewBase::OneViewScanline::OneViewScanline(const int _width)
  : width(_width), scan_points(new ScanPoint[_width])
{
  assert( scan_points != nil );
}

                                // Destroy the scanline
ThreeDViewBase::OneViewScanline::~OneViewScanline(void)
{
  assert( scan_points != nil );
  delete [] scan_points;
}

                                // print a scanline point
void ThreeDViewBase::OneViewScanline::print(const card i) const
{
  assert( i < width );
  message("scanline point %d is (%d,%d)\n",i,scan_points[i].v,scan_points[i].color);
}


                                // The following is just an ugly kludge in lieu of member
                                // function templates
class ScanLineAccess
{
protected:
  ThreeDViewBase::OneViewScanline& the_scanline;
public:
  ScanLineAccess(ThreeDViewBase::OneViewScanline& a_scanline) : the_scanline(a_scanline) {}
  int q_width(void) const                { return the_scanline.q_width(); }
                                    // This is the real reason for this class
  ThreeDViewBase::OneViewScanline::ScanPoint * q_scan_begin(void) const
                                          { return the_scanline.scan_points; }
};

                                // Note that the entire iteration would be done "in-line"
                                // Class T (an iteratee) is supposed to have a member
                                //    operator () (ScanPoint& curr_point)
                                // to do something meaningful with the current point
                                // (say, fill it in)
template <class T> struct ScanLineIterator : ScanLineAccess
{
  ScanLineIterator(ThreeDViewBase::OneViewScanline& a_scanline) : ScanLineAccess(a_scanline) {}
  ThreeDViewBase::OneViewScanline& get_scanline(void) const { return the_scanline; }
  void for_each(T& iteratee)            // Get the iteratee to do something for each point
   {
     ThreeDViewBase::OneViewScanline::ScanPoint * p = q_scan_begin();
     const ThreeDViewBase::OneViewScanline::ScanPoint * p_end = p + q_width();
     while( p < p_end )
       iteratee(*p++);
   }
};

                                // This is an iterator fow two scanlines
                                // class T has to have a member function 'operator()'
                                // that takes two scanline points (from the previous and
                                // the current scanline) and 'handles' them
template <class T> class TwoScanLinesIterator
{
  const ScanLineAccess from_scanline;
  const ScanLineAccess to_scanline;
public:
  TwoScanLinesIterator(ThreeDViewBase::OneViewScanline& _from_scanline,
                         ThreeDViewBase::OneViewScanline& _to_scanline)
    : from_scanline(_from_scanline), to_scanline(_to_scanline) {}
  void for_each(T& iteratee)
   {
     const ThreeDViewBase::OneViewScanline::ScanPoint * from_p = from_scanline.q_scan_begin();
     const ThreeDViewBase::OneViewScanline::ScanPoint * from_p_end = from_p + from_scanline.q_width();
     const ThreeDViewBase::OneViewScanline::ScanPoint * to_p = to_scanline.q_scan_begin();
     while( from_p < from_p_end )
       iteratee(*from_p++,*to_p++);
   }
};


/*
 *----------------------------------------------------------------------
 *                        Projection algorithm itself
 */

                                    // This is the class that fills out a scanline
                                    // of the view plane for a current depth d
class ThreeDProjector
{
  enum {                                // A range d, the depth, can vary in
         Sight_depth = 140,                        // The farthest the observer can see
         Focus_distance = 40,                    // The closest he can see
                                         // classification of depths into far, medium, near
         far_tier  = Focus_distance+90,            // d > far_tier: far away objects
         near_tier = Focus_distance+70 };

protected:

  typedef ThreeDViewBase::OneViewScanline::ScanPoint Point;        // Just an abbreviation
  int d;                                        // The current projection depth
  const IMAGE& map;
  const int map_width_u, map_height_u;            // Unnormalized (in "units" |g|, that is,
                                                // shifted by 14) xmax and ymax

  const int xe, ye, ze;                            // World coordinates of a view point
  const int gx, gy;                                // Directions of a gaze vector (gz=0)
  const int beta;                                // Local coordinates scaling factor
  const int Ou, Ov;                                // Origin of the local coordinate ref Pt
  const int vmin, vmax;                            // Allowable span of the local v coordinate
  const int z_cutoff, elevation_factor;            // Coefficients to determine elevation from
                                                  // the map's "elevation code"


  int xr, yr;                                    // In units of |g|
  int xinc, yinc;                                // Increment in xu, yu, when u changes by 1
  int zmax, zmin;                                // When z lies within [zmin, zmax], v would
                                                  // _almost_ surely fall within [vmin,vmax]
  const int beta_times_elev;                    // (beta*elevation_factor)
  int beta_over_d;                                // (beta/d)*elevation_factor
  int ze_u;                                        // ze/elevation_factor

                                                // The following are changes (derivatives)
                                                // in the corresponding quantities when
                                                // d changes by 1, precomputed for easy
                                                // reference
  const int xr_d, yr_d, xinc_d, yinc_d, zmax_d, zmin_d;

  int xu, yu;                                    // Current x, y unnormalized


  int clip_v(const int v) const            { return v < vmin ? vmin : v > vmax ? vmax : v; }

                                              // Project the current point on the map into 'point'
                                              // Return FALSE if the projection appears to be
                                              // off-limits
  inline bool project_to_point(Point& point)
  {
    if( xu < 0 || xu >= map_width_u || yu < 0 || yu >= map_height_u )
      return false;                        // The corresponding world point was outside the map
    const int z = point.color = (unsigned)map(yu>>14,xu>>14); // first coordinate is row number
    if( z < z_cutoff )
      return false;                        // The elevation was below the horizon

                                              // Note, even if z was within limits, v still can
                                              // fall outside [vmin,vmax] because of rounding
                                              // errors in computing zmax, zmin
      //    point.v = vmin;
    point.v = z >= zmax ? vmax : z <= zmin ? vmin : clip_v( Ov + ((beta_over_d * (z-ze_u))>>9) );
//    if( Ov > 100 & point.v > 0 ) _error("point.v %d xu %d, yu %d, zmin %d, zmax %d,ze_u %d, z %d d %d", point.v,xu,yu,zmin,zmax,ze_u,z,d);
    return true;
  }

public:
    ThreeDProjector(const ElevationMap& _map, const ViewerPosition& viewer_pos,
                    const ProjectionParameters& projp)
      : map(_map), d(Sight_depth),
        map_width_u(_map.q_ncols()<<14), map_height_u(_map.q_nrows()<<14),
        xe(viewer_pos.xe), ye(viewer_pos.ye), ze(projp.ze),
        gx(viewer_pos.gx), gy(viewer_pos.gy), beta(projp.beta), Ou(projp.Ou), Ov(projp.Ov),
        vmin(0), vmax(projp.q_height()-1),
        z_cutoff(_map.q_z_cutoff()),
        elevation_factor(_map.q_elev_scale()),
        beta_times_elev(beta*elevation_factor),
        ze_u((ze<<9)/elevation_factor),
        xr_d(gx + (gy*Ou)/beta),
        yr_d(gy - (gx*Ou)/beta),
        xinc_d(-gy/beta), yinc_d(gx/beta),
        zmax_d(((vmax-Ov)<<9)/beta_times_elev-1),
        zmin_d(((vmin-Ov)<<9)/beta_times_elev+1)            // to account for roundoff errors

      { xu = xr = (xe<<14) + d*xr_d;                            // In units of |g|
        yu = yr = (ye<<14) + d*yr_d;
        xinc = - (gy*d)/beta;
        yinc = (gx*d)/beta;
        beta_over_d = beta_times_elev/d;
        zmax = ze_u + ((vmax-Ov)<<9)/beta_over_d + 1;
        zmin = ze_u + ((vmin-Ov)<<9)/beta_over_d - 1;
      }

    inline bool can_go_closer(void) const         { return d >= Focus_distance; }
    inline bool q_looking_far(void) const        { return d > far_tier; }
    inline bool q_looking_near(void) const        { return d <= near_tier; }

    void move_closer(void)                    // Prepare to handle a new scanline one step closer
      { d -= 1; xu = (xr -= xr_d); yu = (yr -= yr_d);
        xinc -= xinc_d; yinc -= yinc_d;
        zmax -= zmax_d; zmin -= zmin_d;
        #if SIMULATE_TILT
        ze_u += 1;
        #endif
        beta_over_d = beta_times_elev/d;
      }

    inline void operator () (Point& point)
    {
      if( !project_to_point(point) )
        point.v = vmin, point.color = 0;
      xu += xinc, yu += yinc;
    }
  };

class ThreeDProjectorCached : public ThreeDProjector
{
  typedef ThreeDViewBase::OneViewScanline::ScanPoint Point;        // Just an abbreviation

  Point cached_point;
  int use_cached;                    // when it is 0 or negative, recompute the point
  int cache_inc;

public:
  ThreeDProjectorCached(const ElevationMap& map, const ViewerPosition& viewer_pos,
                        const ProjectionParameters& projp)
      : ThreeDProjector(map,viewer_pos,projp),
        use_cached(0), cache_inc(1) {}

  void move_closer(void)                    // Prepare to handle a new scanline one step closer
  {
    ThreeDProjector::move_closer();
    use_cached = 0;
    if( q_looking_near() )
      cache_inc = 0;                        // disables skipping of the point on close range
  }

  inline void operator () (Point& point)
  {
    if( --use_cached < 0 )                    // if use_cached > 0, return the previous (cached)
    {                                        // point rather than really computing it...
      if( !project_to_point(cached_point) )
        cached_point.v = vmin, cached_point.color = 0;
      use_cached = cache_inc;
    }
    point = cached_point;
    xu += xinc, yu += yinc;
  }
};

                                        // Clean the last scan line
struct LastScanLineCleaner
{
  inline void operator () (ThreeDViewBase::OneViewScanline::ScanPoint& point)
  {
    point.v = 0, point.color = 0;
  }
};


                                // Connect two scanlines, sl0 to sl1, drawing
                                // vertical lines between the corresponding points
                                // (if the lines are visible)
                                // We really need some trust here, that is, we want to
                                // believe the projector made v lie within necessary limits
                                // See the explanation for color interpolation algorithm
                                // above
class ScanLines_Connector
{
  const PixMapHandle pixmap;
  const int maxrow;                            // max row number of the pixmap
  const int bytes_per_row;                    // bytes per row in the pixmap
  unsigned char * orig_pixels_ptr;            // Ptr to the beginning of the PixMap
  unsigned char * pixels_ptr;                // current pixmap pointer

                                        // This is a real scanpoint connector
                                        // connects from_point to to_point
                                        // (provided the line is visible, that is,
                                        // on the descending slope)
                                        // Note that we assumed that larger v mean
                                        // points close to the top. On a Mac though
                                        // larger row number corresponds to the
                                        // bottom points. That is, we need a flip
                                        // !! Further optimization: build a scanlines
                                        // row index (to get rid of multiplications by
                                        // bytes_per_row below)
  inline void draw_vert_line(const ThreeDViewBase::OneViewScanline::ScanPoint& from_point,
                                 const ThreeDViewBase::OneViewScanline::ScanPoint& to_point)
  {
    if( to_point.v > from_point.v )
      return;                        // The line is on the ascending slope and *will* be obscured later
                                    // Keep in mind that to_v <= from_v now
    if( from_point.v <= 0 )
      return;                        // means both v and vold are out-of-picture

    if( from_point.color == 0 && to_point.color == 0 )    // The line is in "background" color
      return;                                            // won't waste our time

    assert( to_point.v >= 0 && from_point.v <= maxrow );
                                    // Now, 0 <= to_v <= from_v
    int v_span = from_point.v - to_point.v;
    if( v_span == 0 )                        // one-pixel-long line
    {
      pixels_ptr[(maxrow-to_point.v) * bytes_per_row] = to_point.color;
      return;
    }
    if( v_span == 1 )                // The vertical line is only two-pixel tall
    {                                // No color-interpolation necessary
      unsigned char * pixp = pixels_ptr + (maxrow-to_point.v) * bytes_per_row;
      *pixp = to_point.color;
      *(pixp-bytes_per_row) = from_point.color;
      return;
    }
    int color_span = from_point.color - to_point.color;
                    // The line is of the same color, no color-interpolation
                    // or the line goes to the background: won't interpolate either
    if( color_span == 0 || to_point.color == 0)
    {
      const int color = from_point.color;
      unsigned char * pixp_beg = pixels_ptr + (maxrow-to_point.v) * bytes_per_row;
      for(register int v=v_span; v >= 0; v--, pixp_beg -= bytes_per_row)
          *pixp_beg = color;
      return;
    }

    int acc_dec = v_span;                // by how much to decrement acc
    int curr_color_inc = 1;                // and increment curr_color for the next point
    if( color_span > v_span )
    {
      acc_dec = 0; curr_color_inc = 0;
      for(register int acc=1+color_span; acc >= v_span; acc -= v_span)
         acc_dec += v_span, curr_color_inc++;
    }
    if( color_span < 0 )
    {
      color_span = -color_span;
      if( color_span > v_span )
      {
        acc_dec = 0; curr_color_inc = 0;
        for(register int acc=1+color_span; acc >= v_span; acc -= v_span)
           acc_dec += v_span, curr_color_inc--;
      }
      else
        curr_color_inc = -1;
    }
    int curr_color = to_point.color;
    int acc = 1;
    unsigned char * pixp_beg = pixels_ptr + (maxrow-to_point.v) * bytes_per_row;
    for(register int v=v_span; v >= 0; v--, pixp_beg -= bytes_per_row)
    {
      *pixp_beg = curr_color;
      acc += color_span;
      while( acc >= v_span )                // Note this loop runs at most twice
        acc -= acc_dec,
        curr_color += curr_color_inc;
    }
    #if 0
    if( abs(*(pixp_beg+bytes_per_row) - from_point.color) >= abs(curr_color_inc) )
    _error("from (%d,%d) to (%d,%d), color_span %d, acc_dec %d, inc %d last %d\n",
            from_point.v,from_point.color,to_point.v,to_point.color,
            color_span, acc_dec, curr_color_inc,
            *(pixp_beg+bytes_per_row));
    #endif
  }

public:
  ScanLines_Connector(const OffScreenBuffer& canvas)
    : pixmap(canvas.get_pixmap()), maxrow(canvas.height()-1),
      bytes_per_row(canvas.bytes_per_row())
  {
    assert( LockPixels(pixmap) );
    pixels_ptr = orig_pixels_ptr = (unsigned char *)GetPixBaseAddr(pixmap);
  }
  ~ScanLines_Connector(void)    {  UnlockPixels(pixmap); }

  void reset(void)            { pixels_ptr = orig_pixels_ptr; }    // Reset for the new run

                                        // A connector-iteratee
                                        // It is given to scan points (for every column in turn)
                                        // and its' job to connect the point via a vertical
                                        // line
  void operator() (const ThreeDViewBase::OneViewScanline::ScanPoint& from_point,
                     const ThreeDViewBase::OneViewScanline::ScanPoint& to_point)
  {
    draw_vert_line(from_point,to_point);    // The real job is done by somebody else
    pixels_ptr++;                            // we merely advance the (cached) curr scan pointer
  }

};

/*
 *----------------------------------------------------------------------
 *                Root module to draw a projection of the world
 *                        in an offscreen graf world
 */

                                    // Prepare the buffer (eg, clean it)
void ThreeDViewBase::prepaint(void)
{
  clear();
}

template <class LineProjector>
void ThreeDView<LineProjector>::project(void)
{
  prepaint();

  LineProjector projector(map, viewer_pos, projection_parms);

  ScanLineIterator<LineProjector> fill_scanline_1(scanline1);
  ScanLineIterator<LineProjector> fill_scanline_2(scanline2);
  ScanLineIterator<LineProjector> * prev_scanline = &fill_scanline_1,
                                        * curr_scanline = &fill_scanline_2;

  ScanLines_Connector scanlines_connectee(*this);
                                                // note curr_scanline ^ pointer_diff
                                                // gives prev_scanline, and
                                                // vice versa
  const int pointer_diff = (long int)curr_scanline ^ (long int)prev_scanline;

  (*prev_scanline).for_each(projector);

  while( projector.can_go_closer() )
  {
    (*curr_scanline).for_each(projector);
    TwoScanLinesIterator<ScanLines_Connector>
        scanlines_connector(prev_scanline->get_scanline(),curr_scanline->get_scanline());
    scanlines_connector.for_each(scanlines_connectee);

                                            // exchange pointers
    (long int&)prev_scanline ^= pointer_diff;
    (long int&)curr_scanline ^= pointer_diff;

    projector.move_closer();
    if( !projector.q_looking_near() )        // Unless we're projecting near-view objects,
      projector.move_closer();                // we decrement the depth at a faster rate...
    if( projector.q_looking_far() )
      projector.move_closer();
    scanlines_connectee.reset();
  }



  {                                            // Handle the last scanline
   ScanLineIterator<LastScanLineCleaner> last_scanline_filler(curr_scanline->get_scanline());
   last_scanline_filler.for_each(LastScanLineCleaner());
    TwoScanLinesIterator<ScanLines_Connector>
        scanlines_connector(prev_scanline->get_scanline(),last_scanline_filler.get_scanline());
    scanlines_connector.for_each(scanlines_connectee);
  }
}

                        // Tell me about this view
void ThreeDViewBase::dump(const char title []) const
{
  map.dump();
  viewer_pos.dump();
  projection_parms.dump();
}