"Fossies" - the Fresh Open Source Software Archive

Member "aqsis-1.8.2/libs/core/bilinear.h" (24 Aug 2012, 9103 Bytes) of package /linux/privat/aqsis-1.8.2.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file.

    1 // Aqsis
    2 // Copyright (C) 1997 - 2001, Paul C. Gregory
    3 //
    4 // Contact: pgregory@aqsis.org
    5 //
    6 // This library is free software; you can redistribute it and/or
    7 // modify it under the terms of the GNU General Public
    8 // License as published by the Free Software Foundation; either
    9 // version 2 of the License, or (at your option) any later version.
   10 //
   11 // This library is distributed in the hope that it will be useful,
   12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
   13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   14 // General Public License for more details.
   15 //
   16 // You should have received a copy of the GNU General Public
   17 // License along with this library; if not, write to the Free Software
   18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
   19 
   20 
   21 /** \file
   22  * \brief Forward and reverse bilinear mappings.
   23  * \author Paul C. Gregory (pgregory@aqsis.org)
   24  * \author Chris Foster
   25  */
   26 
   27 //? Is .h included already?
   28 #ifndef BILINEAR_H_INCLUDED
   29 #define BILINEAR_H_INCLUDED 1
   30 
   31 #include    <aqsis/aqsis.h>
   32 #include    <aqsis/math/vector2d.h>
   33 
   34 namespace Aqsis {
   35 
   36 
   37 /** \brief Functor for inverse bilinear mapping
   38  *
   39  * Forward bilinear interpolation is extremely staightforward: Given four
   40  * corners of a bilinear patch, A,B,C and D, we define the bilinear surface
   41  * between them at parameter value (u,v) to be
   42  *
   43  * P(u,v) := (1-v)*((1-u)*A + u*B) + v*((1-u)*C + u*D)
   44  *
   45  * That is, linear interpolation in the u-direction along opposite edges,
   46  * followed by linear interpolation in the v-direction between the two
   47  * resulting points.
   48  *
   49  * So, we can get a point P easily from the (u,v) coordinates.  The problem
   50  * solved here is the reverse: get (u,v) from the (x,y) coordinates of P.  This
   51  * class assumes that the user may want to do this lookup operation many times
   52  * per micropolygon, so it holds a cache of coefficients for the lookup
   53  * function.
   54  *
   55  * Here's a schematic showing the bilinear patch ABCD, along with two important
   56  * edge vectors E and F used in the calculation.
   57  *
   58  * \verbatim
   59  *
   60  *                 C-----------------D
   61  *   v-        ^   |                /
   62  *   direction |   |               /
   63  *             |   |              /
   64  *          F=C-A  |   .P        /
   65  *             |   |            /
   66  *             |   |           /
   67  *             |   |          /
   68  *             |   |         /
   69  *                 A--------B
   70  *  
   71  *                  ------->
   72  *                  E = B-A
   73  *  
   74  *                u-direction ->
   75  *
   76  * \endverbatim
   77  *
   78  * The inverse of the equation for P(u,v) may be calculated analytically, but
   79  * there's a bunch of special cases which result in numerical instability, each
   80  * of which needs to be carefully coded in.  Here we use two iterations of
   81  * Newton's method instead to find (u,v) such that the residual function
   82  *
   83  *   f(u,v) := (1-v)*((1-u)*A + u*B) + v*((1-u)*C + u*D) - P
   84  *
   85  * is zero.  This works very reliably when ABCD is pretty much rectangular with
   86  * two iterations giving accuracy to three or four decimal places, and only
   87  * appears to converge slowly when two adacent sides are parallel.
   88  *
   89  * Tests show that the analytical solution to the equations and two iterations
   90  * of Newton's method are very similar in speed, so we take the latter for
   91  * implementation simplicity.  (See the prototypes directory for both
   92  * implementations.)
   93  *
   94  * NOTE: The speed of this code is pretty critical for the micropolygon
   95  * sampling implementation, since it's inside one of the the innermost loops.
   96  * The accuracy however probably isn't so much of a big deal since
   97  * micropolygons typically are very small on screen.
   98  */
   99 class CqInvBilinear
  100 {
  101     public:
  102         /// Construct a lookup functor with zeros for the vertices.
  103         CqInvBilinear();
  104 
  105         /** Use the given vertices for inverse bilinear lookups.
  106          *
  107          * \see setVertices for the vertex ordering.
  108          */
  109         CqInvBilinear(CqVector2D A, CqVector2D B, CqVector2D C, CqVector2D D);
  110 
  111         /** \brief Reset the micropolygon vertices.
  112          *
  113          * The ordering of vertices is as follows:
  114          *
  115          *   C---D
  116          *   |   |
  117          *   A---B
  118          */
  119         void setVertices(CqVector2D A, CqVector2D B, CqVector2D C, CqVector2D D);
  120 
  121         /** \brief Perform the inverse bilinear mapping
  122          *
  123          * The mapping takes P = (x,y) and returns the original (u,v)
  124          * parameters of the bilinear patch which would map to P under the
  125          * usual bilinear interpolation scheme.
  126          */
  127         CqVector2D operator()(CqVector2D P) const;
  128 
  129     private:
  130         template<bool unsafeInvert>
  131         static CqVector2D solve(CqVector2D M1, CqVector2D M2, CqVector2D b);
  132 
  133         CqVector2D bilinEval(CqVector2D uv) const;
  134 
  135         CqVector2D m_A;
  136         CqVector2D m_E;
  137         CqVector2D m_F;
  138         CqVector2D m_G;
  139         bool m_linear;
  140 };
  141 
  142 
  143 //---------------------------------------------------------------------
  144 /** Bilinearly evalute the four specified values at the specified intervals.
  145  * \attention The type to be interpolated must support operator+, operator- and operator*.
  146  * \param A min u min v.
  147  * \param B max u min v.
  148  * \param C min u max v.
  149  * \param D max u max v.
  150  * \param s the fraction in the u direction.
  151  * \param t the fraction in the v direction.
  152  * \return the interpolated value.
  153  */
  154 
  155 template <class T>
  156 T BilinearEvaluate( const T& A, const T& B, const T& C, const T& D, TqFloat s, TqFloat t );
  157 
  158 
  159 /** \brief Bilinear interpolation, fast version.
  160  *
  161  * Bilinear interpolation over a "quadrilateral" of values arranged in the
  162  * standard order as follows:
  163  *
  164  * \verbatim
  165  *
  166  *   C---D    ^
  167  *   |   |    | v-direction [ uv.y() ]
  168  *   A---B    |
  169  *
  170  *   u-direction -->
  171  *   [ uv.x() ]
  172  *
  173  * \endverbatim
  174  *
  175  * This version performs no checks for whether u and v are between 0 and 1 and
  176  * may give odd results if they're not.
  177  *
  178  * \param A,B,C,D - quadrilateral corners.
  179  * \param uv - Parametric (u,v) coordinates, should be in the interval [0,1]
  180  */
  181 template<typename T>
  182 inline T bilerp(T A, T B, T C, T D, CqVector2D uv);
  183 
  184 
  185 
  186 //==============================================================================
  187 // Implementation details.
  188 //==============================================================================
  189 template <class T>
  190 inline T BilinearEvaluate( const T& A, const T& B, const T& C, const T& D, TqFloat s, TqFloat t )
  191 {
  192     T AB, CD;
  193     // Work out where the u points are first, then linear interpolate the v value.
  194     if ( s <= 0.0 )
  195     {
  196         AB = A;
  197         CD = C;
  198     }
  199     else
  200     {
  201         if ( s >= 1.0 )
  202         {
  203             AB = B;
  204             CD = D;
  205         }
  206         else
  207         {
  208             AB = static_cast<T>( ( B - A ) * s + A );
  209             CD = static_cast<T>( ( D - C ) * s + C );
  210         }
  211     }
  212 
  213     T R;
  214     if ( t <= 0.0 )
  215         R = AB;
  216     else
  217     {
  218         if ( t >= 1.0 )
  219             R = CD;
  220         else
  221             R = static_cast<T>( ( CD - AB ) * t + AB );
  222     }
  223 
  224     return ( R );
  225 }
  226 
  227 
  228 template<typename T>
  229 inline T bilerp(T A, T B, T C, T D, CqVector2D uv)
  230 {
  231     TqFloat w0 = (1-uv.y())*(1-uv.x());
  232     TqFloat w1 = (1-uv.y())*uv.x();
  233     TqFloat w2 = uv.y()*(1-uv.x());
  234     TqFloat w3 = uv.y()*uv.x();
  235     return w0*A + w1*B + w2*C + w3*D;
  236 }
  237 
  238 
  239 //------------------------------------------------------------------------------
  240 // CqInvBilinear implementation
  241 inline CqInvBilinear::CqInvBilinear()
  242     : m_A(),
  243     m_E(),
  244     m_F(),
  245     m_G(),
  246     m_linear(false)
  247 {}
  248 
  249 inline CqInvBilinear::CqInvBilinear(CqVector2D A, CqVector2D B,
  250                                     CqVector2D C, CqVector2D D)
  251     : m_A(),
  252     m_E(),
  253     m_F(),
  254     m_G(),
  255     m_linear(false)
  256 {
  257     setVertices(A,B,C,D);
  258 }
  259 
  260 inline void CqInvBilinear::setVertices(CqVector2D A, CqVector2D B,
  261                                        CqVector2D C, CqVector2D D)
  262 {
  263     m_A = A,
  264     m_E = B-A;
  265     m_F = C-A;
  266     m_G = -m_E-C+D;
  267     m_linear = false;
  268     // Determine whether the micropolygon is almost-rectangular.  If it is, we
  269     // set the m_linear flag as a hint to the solver that we only need to solve
  270     // linear equations (requiring only one iteration of Newton's method).
  271     TqFloat patchSize = max(maxNorm(m_F), maxNorm(m_E));
  272     TqFloat irregularity = maxNorm(m_G);
  273     if(irregularity < 1e-2*patchSize)
  274         m_linear = true;
  275 }
  276 
  277 inline CqVector2D CqInvBilinear::operator()(CqVector2D P) const
  278 {
  279     // Start at centre of the micropoly & do one or two iterations of Newton's
  280     // method to solve for (u,v).
  281     CqVector2D uv(0.5, 0.5);
  282     uv -= solve<true>(m_E + m_G*uv.y(), m_F + m_G*uv.x(), bilinEval(uv)-P);
  283     if(!m_linear)
  284     {
  285         // The second iteration is only used if we know that the micropolygon
  286         // is non-rectangular.
  287         uv -= solve<false>(m_E + m_G*uv.y(), m_F + m_G*uv.x(), bilinEval(uv)-P);
  288     }
  289     return uv;
  290 }
  291 
  292 /** \brief Solve the linear equation M*x = b  with the matrix M = [M1 M2]
  293  *
  294  * The unsafeInvert parameter indicates whether we should check that an inverse
  295  * exists before trying to find it.  If it doesn't exist, we return the zero
  296  * vector.
  297  */
  298 template<bool unsafeInvert>
  299 inline CqVector2D CqInvBilinear::solve(CqVector2D M1, CqVector2D M2,
  300                                               CqVector2D b)
  301 {
  302     TqFloat det = cross(M1, M2);
  303     if(unsafeInvert || det != 0) det = 1/det;
  304     return det * CqVector2D(cross(b, M2), -cross(b, M1));
  305 }
  306 /// Evaluate the bilinear function at the coordinates (u,v)
  307 inline CqVector2D CqInvBilinear::bilinEval(CqVector2D uv) const
  308 {
  309     return m_A + m_E*uv.x() + m_F*uv.y() + m_G*uv.x()*uv.y();
  310 }
  311 
  312 
  313 } // namespace Aqsis
  314 
  315 #endif  // !BILINEAR_H_INCLUDED