TFRasterize.h

Go to the documentation of this file.

//------------------------------------------------------------------------
//
//   Joe Kniss
//     3-20-03
//                   ________    ____   ___ 
//                  |        \  /    | /  /
//                  +---+     \/     |/  /
//                  +--+|  |\    /|     < 
//                  |  ||  | \  / |  |\  \ 
//                  |      |  \/  |  | \  \ 
//                   \_____|      |__|  \__\
//                       Copyright  2003 
//                      Joe Michael Kniss
//                   <<< jmk@cs.utah.edu >>>
//               "All Your Base are Belong to Us"
//-------------------------------------------------------------------------

/// TFRasterize.h

#ifndef __SIMIAN_TF_RASTERIZE_DOT_H
#define __SIMIAN_TF_RASTERIZE_DOT_H

#include <mathGutz.h>
#include <iostream>

using namespace std; /// bad

/// a triangle class if you don't already have somthing
/// to represent a triangle

template<class T>
class Triangle {
public:
   typedef gutz::vec3<T> Vec3;
   enum N_TEX_COORDS {
      N_TEX = 8
   };

   Triangle()
   {
      for(int i=0; i<3; ++i)
         pts[i] = 0;
      for(int i=0; i<N_TEX; ++i)
         for(int j=0; j<3; ++j)
            tex[i][j] = 0;
   }

   Triangle(Vec3 *p1, Vec3 *p2, Vec3 *p3)
   {
      if((!pi) || (!p2) || (!p3))
      {
         std::cerr << "Triangle(), got null point!" << std::endl;
      }
      pts[0] = p1; pts[1] = p2; pts[2] = p3;
      for(int i=0; i<N_TEX; ++i)
         for(int j=0; j<3; ++j)
            tex[i][j] = 0;
   }
   Triangle(const Triangle &t)
   {
      for(int i=0; i<3; ++i) pts[i] = t.pts[i];
      for(int i=0; i<N_TEX; ++i)
         for(int j=0; j<3; ++j)
            tex[i][j] = t.tex[i][j];
   }
   ///~Triangle();  /// just using default for now
   Triangle &operator=(const Triangle &t)
   {
      for(int i=0; i<3; ++i) pts[i] = t.pts[i];
      for(int i=0; i<N_TEX; ++i)
         for(int j=0; j<3; ++j)
            tex[i][j] = t.tex[i][j];
      return *this;
   }
   void setPoints(Vec3 *p1, Vec3 *p2, Vec3 *p3)
   {
      if((!p1) || (!p2) || (!p3))
      {
         std::cerr << "Triangle(), got null point!" << std::endl;
      }
      pts[0] = p1; pts[1] = p2; pts[2] = p3;
   }
   void setTexCoord(int tc, Vec3 *t1, Vec3 *t2, Vec3 *t3)
   {
      if((!t1) || (!t2) || (!t3))
      {
         std::cerr << "Triangle(), got null tc!" << std::endl;
      }
      tex[tc][0] = t1; tex[tc][1] = t2; tex[tc][2] = t3;
   }

   Vec3 &getTex(int tcNum, int ptNum)             { return *(tex[tcNum][ptNum]); }
   const Vec3 &getTex(int tcNum, int ptNum) const { return *(tex[tcNum][ptNum]); }

   /// a convenient accessor!
   Vec3       &operator[](int i)       { return *(pts[i]); }
   const Vec3 &operator[](int i) const { return *(pts[i]); }

   /// three points, verticies of triangle
   Vec3  *pts[3];
   /// N_TEX coordinate sets, one Vec3 coordinate for each point
   Vec3  *tex[N_TEX][3];
};


/// TriRasterInfo:
///   a class to help you setup rasterization.
///  Works kind of like an iterator, here is an example:
/// \code
///  TriRasterInfo ri(imageSizeX, imageSizeY);
///  for(...)/// for each triangle
///  {  
///     ri.setTri(someTriangle);
///     while(!ri.end())
///     {
///         /// position of fragment is ri.x() and ri.y()
///         /// you can interpolate values using the interp*() functions
///         ++ri;  /// increment to the next pixel
///     }
///  }
/// \endcode
///
/// any structure that returns a point using "[]" will 
/// work as the triangle type, you could just use an 
/// array of vec3<>'s.
///   TODO: make this work for c arrays of c arrays??
template< class T, class TT = Triangle<T> >
class TriRasterInfo {
public:
   typedef gutz::vec4<T> Vec4;
   typedef gutz::vec3<T> Vec3;
   typedef gutz::vec2<T> Vec2;
   typedef gutz::ray3<T> Ray3;
   typedef TT   Tri3;

   TriRasterInfo(int sx, int sy) : yPos(0), maxY(0), _sx(sx), _sy(sy) {}

   inline int  getSizeX() const { return _sx; }
   inline int  getSizeY() const { return _sy; }

   /// uncomment this if there are issues with non virtual 
   /// destructors
   ///~TriRasterInfo() {}

   /// use it like an iterator!... kind of
   inline bool end() const { return yPos >= maxY; }

   inline TriRasterInfo &operator++() 
   {
      incX();
      if(xPos >= maxX)
      {
         //cerr << " scan line y = " << yPos << " : xmin, xmax => " << minX << " "  << maxX << endl;
         incY();
      }
      return *this;
   }

   /// returns true if the triangle has NO renderable pixels
   /// returns false if the has any renderable pixels
   inline bool        setTri(Tri3 const *t) { return setUp(t); }
   inline Tri3 const *getTri() const { return tri; }

   inline T x() const { return xPos; }
   inline T y() const { return yPos; }
   
   inline void reset() { setUp(this->tri); }

   ///////////////////////////////////////////////////////////
   /// get the minimum x for the CURRENT scan line (at y())
   int   getMinX() const { return minX; }
   /// get the starting scan line (min y pos)
   int   getMinY() const { return minY; }

   ///////////////////////////////////////////////////////////
   /// get the maximum x for the CURRENT scan line (at y())
   int   getMaxX() const { return maxX; }
   /// get the ending scan line (max y pos)
   int   getMaxY() const { return maxY; }

   ///////////////////////////////////////////////////////////
   /// get the max x position in image for some scan line yp
   inline int getMaxX(T yp) const
   {
      if(yPos > yswapEnd) /// lower segment
      {
         const T t = intersectY(yPos, r[2]);
         return clamp(0, int(t)-1, _sx-1);
      }
      else  /// upper segment
      {
         const T t = intersectY(yPos, r[ee]);
         return clamp(0, int(t)-1, _sx-1);
      }
   }

   ///////////////////////////////////////////////////////////
   /// get the min x position in image for some scanline yp
   inline int getMinX(T yp) const 
   {
      if(yPos > yswapStart)
      {
         const T t = intersectY(yPos, r[2]);
         return clamp(0, int(t), _sx-1);
      }
      else
      {
         T t = intersectY(yPos, r[se]);
         return clamp(0, int(t), _sx-1);
      }
   }


   /// increment y
   inline void incY() 
   {
      ++yPos;
      minX = getMinX(yPos);
      xPos = T(minX);
      maxX = getMaxX(yPos);
      setBary(T(xPos), T(yPos));
   }

   /// increment x
   inline void incX()
   {
      alpha += xa;
      beta  += xb;
      gamma += xg;
      ++xPos;
   }

   inline T getAlpha(const T x, const T y) const
   { return x * xa + y * ya + ac; }

   inline T getBeta(const T x, const T y) const
   { return x * xb + y * yb + bc; }

   inline T getGamma(const T x, const T y) const
   { return x * xg + y * yg + gc; }

   ////////////////////////////
   /// update barycentric coordinates for the current x(),y() positions
   /// alpha, beta, gamma will be set based on
   ///   x,y pos you give
   ////////////////////////////
   inline void setBary(const T x, const T y)
   {
      alpha = getAlpha(x,y);
      beta  = getBeta(x,y);
      gamma = getGamma(x,y);
   }
   /// quick check if current pos (x(),y()) is inside triangle
   inline bool validBary() const
   {
      return (alpha >= 0) && (beta >= 0) && (gamma >= 0);
   }

   /// interpolation on a triangle using barycentric coordinates
   
   /// interpolate a 4 vector using the current barycentric coordinates
   template< class VT > //< any vector type even a c array
   inline Vec4 interp4(const VT &v0, const VT &v1, const VT &v3) const
   {
      return Vec4(v0[0] * alpha + v1[0] * beta + v2[0] * gamma,
                  v0[1] * alpha + v2[1] * beta + v2[1] * gamma,
                  v0[2] * alpha + v2[2] * beta + v2[3] * gamma,
                  v0[3] * alpha + v2[3] * beta + v2[3] * gamma);
   }
   /// interpolate a 3 vector using the current barycentric coordinates
   template< class VT > //< any vector type even a c array
   inline Vec3 interp3(const VT &v0, const VT &v1, const VT &v2) const
   {
      return Vec3(v0[0] * alpha + v1[0] * beta + v2[0] * gamma,
                  v0[1] * alpha + v2[1] * beta + v2[1] * gamma,
                  v0[2] * alpha + v2[2] * beta + v2[3] * gamma);
   }
   /// interpolate a 2 vector using the current barycentric coordinates
   template< class VT > //< any vector type even a c array
   inline Vec2 interp2(const VT &v0, const VT &v1, const VT &v2) const
   {
      return Vec3(v0[0] * alpha + v1[0] * beta + v2[0] * gamma,
                  v0[1] * alpha + v2[1] * beta + v2[1] * gamma);
   }
   /// interpolate a scalar using the current barycentric coordinates
   template< class TP > ///< any type
   inline TP interp(const TP &v0, const TP &v1, const TP &v2) const
   {
      return v0 * alpha + v1 * beta + v2 * gamma;
   }



public: /// but NOT FOR GENERAL USE

   ////////////////////////////
   /// public data
   ////////////////////////////
   T    alpha,beta,gamma;
   T    yPos;  //< current Y position
   T    xPos;  //< current X position
   /// min/max values for setting up loops
   int  minY;
   int  maxY;
   int  minX;
   int  maxX;
   Tri3 const *tri;

   /////////////////////////////////////////////////////////////
   /// setup for rasterizaton, set the triangle being rasterized
   ///   NOT FOR GENERAL USE, use setTri()
   /////////////////////////////////////////////////////////////
   /// \code
   ///  2 basic cases for "tight rendering" 
   ///      case 1    case 2            Loose/"de-implemented"
   ///       p0=sp    p0=sp             x---x
   ///     r0  o        o  r0           |o  |
   ///     se /| r1  r1 |\ ee           ||\ |
   ///  p1-> o | ee  se | o <- p1=mp    || o|
   ///     r2 \|        |/ r2           ||/ |
   ///     me  o        o  me           |o  |
   ///       p2=ep    p2=ep             x---x
   /// \endcode 
   ///   There is one additional check made to insure that the 
   ///     start and end edges are not horizontal, if they are we
   ///     swap immediately. If two are horizontal, we don't have
   ///     any pixels to be rendered.
   /// Users need not worry about this, it's just an explanation of
   ///  how this function sets things up.
   /////////////////////////////////////////////////////////////
   
   inline bool setUp(Tri3 const *t3)
   {
      /// set internal triangle
      this->tri = t3;
      /// alias Tri* to a Tri, for easier access
      const Tri3 &t = (*t3);
      //cerr << " rasterizing triangle 0 = " << t[0] << "  1= " << t[1] << "  2= " << t[2] << endl;
      /// which of the triangle points is smallest in Y (by index)
      int sp = mini3(t[0].y, t[1].y, t[2].y);
      int ep = maxi3(t[0].y, t[1].y, t[2].y); //< largest y (by index)
      int mp = 0; //< middle y (by index)
      if(((0==sp) && (2==ep)) || ((2==sp) && (0==ep))) mp = 1;
      if(((1==sp) && (0==ep)) || ((0==sp) && (1==ep))) mp = 2;
      //cerr << " order " << sp << " " << mp << " " << ep << endl;
      /// Init bounding box ( and starting point for loose ) (by integer value)
      minY = clamp(0,int(min3(t[0].y, t[1].y, t[2].y)), _sy-1);
      maxY = clamp(0,int(max3(t[0].y, t[1].y, t[2].y))-1, _sy-1);
      /// need to capture this to see if we have any usefull pixels
      int tminX = clamp(0,int(min3(t[0].x, t[1].x, t[2].x)), _sx-1);
      int tmaxX = clamp(0,int(max3(t[0].x, t[1].x, t[2].x))-1, _sx-1);
      /// starting Y position, first scan line
      yPos = T(minY);
      /// build rays for tight rendering.
      ///  rays always point in +y direction
      ///  these correspond to case 1 above
      r[0] = Ray3(t[sp], t[mp] - t[sp]);
      r[1] = Ray3(t[sp], t[ep] - t[sp]);
      r[2] = Ray3(t[mp], t[ep] - t[mp]);
      /// which case are we in for "tight" rendering
      //cerr << " swaps " << t[sp].y << " " << t[mp].y << " " << t[ep].y << endl;
      yswapStart = T( int(t[mp].y) );   /// case 1
      yswapEnd   = T( maxY );  
      se = 0;
      ee = 1;
      if(r[se].d.y == 0)  /// start edge is horizontal
      {
         se = 2;
      }
      if(r[ee].d.y == 0)  /// end edge is horizontal
      {
         ee = 2;
      }
      /// init the tmp ray, only the .p.y value will be modified...
      /// should ride the x=0 axis, point toward (1,0,0)
      ///tmpr = Ray3(Vec3(T(0.0),r[se].p.y + r[se].d.y/T(2.0),T(0.0)), Vec3(T(1.0), T(0.0), T(0.0)));
      //cerr << " edge start " << r[se].p << "     " << r[se].d << endl;
      //cerr << " edge next  " << r[ee].p << "     " << r[ee].d << endl;
      //cerr << " edge swap  " << r[2].p << "     " << r[2].d << endl;
      //cerr << endl;
      //cerr << " tmp ray " << tmpr.p << "  =  " << tmpr.d << endl;
      //cerr << " intersect " << intersect2D(tmpr,r[se]) << "   " << intersect2D(tmpr,r[ee]) << endl;
      T typ = r[se].p.y + r[se].d.y/T(2.0);
      if(intersectY(typ,r[se]) > intersectY(typ,r[ee])) /// case 2
      {
         //cerr << " swaping edge order " << endl;
         int tmp = se;
         se = ee;
         ee = tmp;
         yswapStart = T( maxY );
         yswapEnd   = T( int(t[mp].y) );
      }
      /// now that we have edge orders figured out... 
      /// set the current x position
      minX = int(clamp(0,intersectY(yPos,r[se]),_sx-1));
      xPos = T(minX);
      maxX = int(clamp(0,intersectY(yPos,r[ee])-1,_sx-1));
      //cerr << " edge start " << r[se].p << "     " << r[se].d << endl;
      //cerr << " edge next  " << r[ee].p << "     " << r[ee].d << endl;
      //cerr << " edge swap  " << r[2].p << "     " << r[2].d << endl;
      //cerr << " yswapS = " << yswapStart << endl;
      //cerr << " yswapE = " << yswapEnd << endl;

      //cerr << " xdims " << minX << "  " << maxX << endl;
      //cerr << " ydims " << minY << "  " << maxY << endl;
      /// precompute barycentric constants
      setupBary(t[0],t[1],t[2]);
      /// initialize them
      setBary(T(xPos),T(yPos));
      /// report if there are any pixels that need to be rendered
      return !(minY < maxY) || !(tminX < tmaxX);
   }

   /// shouldn't create without a size!!!
   TriRasterInfo() : yPos(0), maxY(0), _sx(0), _sy(0) {}

   inline T intersectY(const T yp, const Ray3 &r) const
   {
      /// since gutz::intersect2D() is more general than 
      /// we need for this algorithm, I optimized it by 
      /// removing terms that will always be zero. And assume that
      /// we don't need to check for divide by zero
      /// r1: p = (0,yp), d=(1,0)
      /// [d1x, d1y] dot Perp([d2x,d2y])
      ///const T div = -r.d.y; /// removed check for x/0
      ///if(!div) return std::numeric_limits<T>::infinity();
      return r.p.x + (r.p.y * r.d.x - yp*r.d.x)/(-r.d.y);
   }

   inline 
   T intersect2D(const Ray3 &ri, const Ray3 &ro) const
   {
      return gutz::intersect2D(ri.p.x, ri.p.y, ri.d.x, ri.d.y,
                               ro.p.x, ro.p.y, ro.d.x, ro.d.y);
   }

   /// image size:
   int _sx, _sy;

   /// pre-computed alpha, beta, gamma constants:
   /// ex: alpha = x * xa + y * ya + ca
   T xa, ya, ac, xb, yb, bc, xg, yg, gc;
   void setupBary(const Vec3 &p0, const Vec3 &p1, const Vec3 &p2)
   {
      /// alpha
      xa = p1.y - p2.y;
      ya = p2.x - p1.x;
      ac  = p1.x * p2.y - p2.x * p1.y;
      /// beta
      xb = p2.y - p0.y;
      yb = p0.x - p2.x;
      bc  = p2.x * p0.y - p0.x * p2.y;
      /// gamma
      xg = p0.y - p1.y;
      yg = p1.x - p0.x;
      gc  = p0.x * p1.y - p1.x * p0.y;
      /// take care of divide
      /// alpha
      T tmp;
      tmp = T(1.0)/(p0.x * xa + p0.y * ya + ac);
      xa *= tmp; ya *= tmp; ac *= tmp;
      tmp = T(1.0)/(p1.x * xb + p1.y * yb + bc);
      xb *= tmp; yb *= tmp; bc *= tmp;
      tmp = T(1.0)/(p2.x * xg + p2.y * yg + gc);
      xg *= tmp; yg *= tmp; gc *= tmp;
   }

   int  se, ee;
   T  yswapStart, yswapEnd;

   Ray3 r[3];

};

#if 0
class TFRasterInfo : public TriRasterInfo<STF::tfSType, Triangle<STF::tfSType> >
{
public:
   typedef STF::tfSType      SType;
   typedef Triangle<SType>   TriType;
   typedef gutz::vec2<SType> Vec2;
   typedef gutz::mat3<SType> Mat3;

   TFRasterInfo() : TriRasterInfo<SType,TriType>() {}
   TFRasterInfo(const TriType &t) : TriRasterInfo<SType,TriType>(t) {}
   TFRasterInfo(const TriType t)  : TriRasterInfo<SType,TriType>(t) {}

   ~TFRasterInfo();

   /// optimize rasterization by leverating these to
   /// speed up transformations from map space to 
   ///  tf-element space
   SType xIncS;
   SType yIncS;
   Vec2  xIncV;
   Vec2  yIncV;
   Mat3  xform;


};
#endif


#endif




#if 0


#include <limits>
#include "TFParams.h"
#include <iostream>

namespace STFR { /// namespaced to prevent pollution (give a hoot?)

typedef STF::tfSType      ST;
typedef STF::tfRangeType  VT;
typedef STF::tfDomainType PT;
typedef STF::tfVec2       V2;

/// a generic rasterization algorithm.
///
/// This is super nasty! But it is so for several reasons... 
/// A.
///   We need to be able to play with blend ops since:
///     1) they are different for each optical property
///     2) there are multiple ways to blend TFElements and 
///        I don't know which one is right
/// B.
///   We also need to be able to generate a 2D TF image from any 
///     combination of 4 optical properties. I choose 4 since it 
///     is the universial type of images & textures
///
/// C.
///   I don't, and you shouldn't, want to write thousands of 
///     rasterization algorithms; one for each Element type, 
///     and one for each combination of values...
///
/// D.
///   We need this to go very very fast.
///
/// This is at least my third iteration on a generic method for rasterizing elements.
/// The first iterations used virtual functions for eval and blend, they went very slow.
///  Unaceptably slow!  We need to avoid virtual function calls in the main loop if we
///  can help it!  I had to sacrifice generality in blending, we must select the blendops
///  at compile time, but they can now be inlined.   This is the primary reason why this
///  algorithm isn't general with respect to the number of optical properties rasterized.
/// There is a helper dispatch #def below, it should be used by a subclass of 
///  TFElement.  It is preferiable that the eval() function be inlined, but it is 
///  absolutely nessesary that it NOT be virtual, that will make rasterization very slow.

template<
   class B1,  ///blenders
   class B2,
   class B3,
   class B4,
   class NT,   ///nrro type 
   class EL    ///element type
>   ///blend type
inline
void rasterize2D4(Nrro::NrroIter<NT> ni, 
                  const gutz::vec2i &posIdx, 
                  const STF::tfRangeIdx &valIdx,
                  const EL *elt)
{
   std::cerr << " hi type = " << typeid(NT).name() << std::endl;

   vec2<NT> mmm; /// min & max for the map

   if(numeric_limits<NT>::epsilon()) /// float/double type
   {
      mmm = vec2<NT>(NT(0.0),NT(1.0));
   }
   else /// integral type
   {
      mmm = vec2<NT>(numeric_limits<NT>::min(), numeric_limits<NT>::max());
   }

   /// idicies 
   const int i1 = valIdx[0];
   const int i2 = valIdx[1];
   const int i3 = valIdx[2];
   const int i4 = valIdx[3];

   /// complete bounding box
   const PT bmin = elt->getMinBox();
   const PT bmax = elt->getMaxBox();

   /// bounding box for the current axies
   const int minX = clamp(0, int(bmin[posIdx.x] * ni->dim(1)), ni->dim(1));
   const int minY = clamp(0, int(bmin[posIdx.y] * ni->dim(2)), ni->dim(2));
   const int maxX = clamp(0, int(bmax[posIdx.x] * ni->dim(1)), ni->dim(1));
   const int maxY = clamp(0, int(bmax[posIdx.y] * ni->dim(2)), ni->dim(2));

   /// tf domain positions and deltas for rasterization
   const ST px = minX/float(ni->dim(1));
   const ST py = minY/float(ni->dim(2));
   const ST dx = 1.0f/float(ni->dim(1));
   const ST dy = 1.0f/float(ni->dim(2));

   /// temps used in inner loop, so we are not constructing tons of stuff
   ST u;
   VT o;
   V2 pos(px,py);

   /// Rasterize!!!
   for(int y = minY; y<= maxY; ++y, pos.y+=dy)  /// y axis
   {
      for(int x = minX; x<= maxX; ++x, pos.x+=dx) /// x axis
      {
         /// over values (values from this element)
         o = elt->eval(pos, posIdx);
#if 1
         /// update the map

         u = g_affine(NT(mmm.x), ni(0,x,y), NT(mmm.y), ST(0.0), ST(1.0));
         //B1::blend(o,u,i1);
         ni(0,x,y) = g_affine(ST(0.0), o[i1], ST(1.0), NT(mmm.x), NT(mmm.y));

         u = g_affine(NT(mmm.x), ni(1,x,y), NT(mmm.y), ST(0.0), ST(1.0));
         //B2::blend(o,u,i2);
         ni(1,x,y) = g_affine(ST(0.0), o[i2], ST(1.0), NT(mmm.x), NT(mmm.y));

         u = g_affine(NT(mmm.x), ni(2,x,y), NT(mmm.y), ST(0.0), ST(1.0));
         //B3::blend(o,u,i3);
         ni(2,x,y) = g_affine(ST(0.0), o[i3], ST(1.0), NT(mmm.x), NT(mmm.y));

         u = g_affine(NT(mmm.x), ni(3,x,y), NT(mmm.y), ST(0.0), ST(1.0));
         //B4::blend(o,u,i4);
         ni(3,x,y) = g_affine(ST(0.0), o[i4], ST(1.0), NT(mmm.x), NT(mmm.y));

#endif
      }
      pos.x = px;
   }

}

/////////////////// Pure nasty! /////////////////////
/// copied and modified from nrroDispatch.h

#define dispatchRast2D4(B1,B2,B3,B4,NSP,PI,VI,EL) \
   \
   switch( NSP->getType() ) { \
   case Nrro::CHAR: \
      STFR::rasterize2D4< B1 , B2 , B3 , B4 >( NSP ->begin<char>(), PI , VI , EL ); \
      break; \
   case Nrro::UCHAR:  \
      STFR::rasterize2D4< B1 , B2 , B3 , B4 >( NSP ->begin<unsigned char>(), PI , VI , EL ); \
      break; \
   case Nrro::SHORT: \
      STFR::rasterize2D4< B1 , B2 , B3 , B4 >( NSP ->begin<short>(), PI , VI , EL ); \
      break; \
   case Nrro::USHORT: \
      STFR::rasterize2D4< B1 , B2 , B3 , B4 >( NSP ->begin<unsigned short>(), PI , VI , EL ); \
      break; \
   case Nrro::INT: \
      STFR::rasterize2D4< B1 , B2 , B3 , B4 >( NSP ->begin<int>(), PI , VI , EL ); \
      break; \
   case Nrro::UINT: \
      STFR::rasterize2D4< B1 , B2 , B3 , B4 >( NSP ->begin<unsigned int>(), PI , VI , EL ); \
      break; \
   case Nrro::FLOAT: \
      STFR::rasterize2D4< B1 , B2 , B3 , B4 >( NSP ->begin<float>(), PI , VI , EL ); \
      break; \
   case Nrro::DOUBLE: \
      STFR::rasterize2D4< B1 , B2 , B3 , B4 >( NSP ->begin<double>(), PI , VI , EL ); \
      break; \
   }

} /// end namespace STFR

#endif

---

Send questions, comments, and bug reports to: