"Fossies" - the Fresh Open Source Software Archive

Member "aqsis-1.8.2/prototypes/newcore/invbilin.h" (24 Aug 2012, 7698 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) 2001, Paul C. Gregory and the other authors and contributors
    3 // All rights reserved.
    4 //
    5 // Redistribution and use in source and binary forms, with or without
    6 // modification, are permitted provided that the following conditions are met:
    7 //
    8 // * Redistributions of source code must retain the above copyright notice,
    9 //   this list of conditions and the following disclaimer.
   10 // * Redistributions in binary form must reproduce the above copyright notice,
   11 //   this list of conditions and the following disclaimer in the documentation
   12 //   and/or other materials provided with the distribution.
   13 // * Neither the name of the software's owners nor the names of its
   14 //   contributors may be used to endorse or promote products derived from this
   15 //   software without specific prior written permission.
   16 //
   17 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
   18 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
   19 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
   20 // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
   21 // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
   22 // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
   23 // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
   24 // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
   25 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
   26 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
   27 // POSSIBILITY OF SUCH DAMAGE.
   28 //
   29 // (This is the New BSD license)
   30 
   31 #ifndef AQSIS_INVBILIN_H_INCLUDED
   32 #define AQSIS_INVBILIN_H_INCLUDED
   33 
   34 #include "util.h"
   35 
   36 namespace Aqsis {
   37 
   38 /** \brief Functor for inverse bilinear mapping
   39  *
   40  * Forward bilinear interpolation is extremely staightforward: Given four
   41  * corners of a bilinear patch, A,B,C and D, we define the bilinear surface
   42  * between them at parameter value (u,v) to be
   43  *
   44  * P(u,v) := (1-v)*((1-u)*A + u*B) + v*((1-u)*C + u*D)
   45  *
   46  * That is, linear interpolation in the u-direction along opposite edges,
   47  * followed by linear interpolation in the v-direction between the two
   48  * resulting points.
   49  *
   50  * So, we can get a point P easily from the (u,v) coordinates.  The problem
   51  * solved here is the reverse: get (u,v) from the (x,y) coordinates of P.  This
   52  * class assumes that the user may want to do this lookup operation many times
   53  * per micropolygon, so it holds a cache of coefficients for the lookup
   54  * function.
   55  *
   56  * Here's a schematic showing the bilinear patch ABCD, along with two important
   57  * edge vectors E and F used in the calculation.
   58  *
   59  * \verbatim
   60  *
   61  *                 C-----------------D
   62  *   v-        ^   |                /
   63  *   direction |   |               /
   64  *             |   |              /
   65  *          F=C-A  |   .P        /
   66  *             |   |            /
   67  *             |   |           /
   68  *             |   |          /
   69  *             |   |         /
   70  *                 A--------B
   71  *  
   72  *                  ------->
   73  *                  E = B-A
   74  *  
   75  *                u-direction ->
   76  *
   77  * \endverbatim
   78  *
   79  * The inverse of the equation for P(u,v) may be calculated analytically, but
   80  * there's a bunch of special cases which result in numerical instability, each
   81  * of which needs to be carefully coded in.  Here we use two iterations of
   82  * Newton's method instead to find (u,v) such that the residual function
   83  *
   84  *   f(u,v) := (1-v)*((1-u)*A + u*B) + v*((1-u)*C + u*D) - P
   85  *
   86  * is zero.  This works very reliably when ABCD is pretty much rectangular with
   87  * two iterations giving accuracy to three or four decimal places, and only
   88  * appears to converge slowly when two adacent sides are parallel.
   89  *
   90  * Tests show that the analytical solution to the equations and two iterations
   91  * of Newton's method are very similar in speed, so we take the latter for
   92  * implementation simplicity.  (See the prototypes directory for both
   93  * implementations.)
   94  *
   95  * NOTE: The speed of this code is pretty critical for the micropolygon
   96  * sampling implementation, since it's inside one of the the innermost loops.
   97  * The accuracy however probably isn't so much of a big deal since
   98  * micropolygons typically are very small on screen.
   99  */
  100 class InvBilin
  101 {
  102     public:
  103         /// Construct a lookup functor with zeros for the vertices.
  104         InvBilin();
  105 
  106         /** \brief Reset the micropolygon vertices.
  107          *
  108          * The ordering of vertices is as follows:
  109          *
  110          *   C---D
  111          *   |   |
  112          *   A---B
  113          */
  114         void init(V2f A, V2f B, V2f C, V2f D);
  115 
  116         /** \brief Perform the inverse bilinear mapping
  117          *
  118          * The mapping takes P = (x,y) and returns the original (u,v)
  119          * parameters of the bilinear patch which would map to P under the
  120          * usual bilinear interpolation scheme.
  121          */
  122         V2f operator()(V2f P) const;
  123 
  124     private:
  125         template<bool unsafeInvert>
  126         static V2f solve(V2f M1, V2f M2, V2f b);
  127 
  128         V2f bilinEval(V2f uv) const;
  129 
  130         V2f m_A;
  131         V2f m_E;
  132         V2f m_F;
  133         V2f m_G;
  134         bool m_linear;
  135 };
  136 
  137 
  138 //==============================================================================
  139 // InvBilin implementation.
  140 inline InvBilin::InvBilin()
  141     // Note that we initialize m_A etc. to zero here only to avoid a bogus gcc
  142     // compiler warning.  These should be properly initialized via the init()
  143     // function.
  144     : m_A(0),
  145     m_E(0),
  146     m_F(0),
  147     m_G(0),
  148     m_linear(false)
  149 {}
  150 
  151 inline void InvBilin::init(V2f A, V2f B, V2f C, V2f D)
  152 {
  153     m_A = A,
  154     m_E = B-A;
  155     m_F = C-A;
  156     m_G = -m_E-C+D;
  157     m_linear = false;
  158     // Determine whether the micropolygon is almost-rectangular.  If it is, we
  159     // set the m_linear flag as a hint to the solver that we only need to solve
  160     // linear equations (requiring only one iteration of Newton's method).
  161     float patchSize = std::max(maxNorm(m_F), maxNorm(m_E));
  162     float irregularity = maxNorm(m_G);
  163     if(irregularity < 1e-2*patchSize)
  164         m_linear = true;
  165 }
  166 
  167 inline V2f InvBilin::operator()(V2f P) const
  168 {
  169     // Start at centre of the micropoly & do one or two iterations of Newton's
  170     // method to solve for (u,v).
  171     V2f uv(0.5, 0.5);
  172     uv -= solve<true>(m_E + m_G*uv.y, m_F + m_G*uv.x, bilinEval(uv)-P);
  173     if(!m_linear)
  174     {
  175         // The second iteration is only used if we know that the micropolygon
  176         // is non-rectangular.
  177         uv -= solve<false>(m_E + m_G*uv.y, m_F + m_G*uv.x, bilinEval(uv)-P);
  178     }
  179     // TODO: Investigate exact solution for InvBilin again!  Newton's method
  180     // can fail to produce the right results sometimes, hence the clamping
  181     // below.
  182     //
  183     // Note that the clamping below has a measurable performance effect (a few
  184     // percent) when rendering z-buffers.
  185     if(uv.x < 0) uv.x = 0;
  186     if(uv.x > 1) uv.x = 1;
  187     if(uv.y < 0) uv.y = 0;
  188     if(uv.y > 1) uv.y = 1;
  189     return uv;
  190 }
  191 
  192 /** \brief Solve the linear equation M*x = b  with the matrix M = [M1 M2]
  193  *
  194  * The unsafeInvert parameter indicates whether we should check that an inverse
  195  * exists before trying to find it.  If it doesn't exist, we return the zero
  196  * vector.
  197  */
  198 template<bool unsafeInvert>
  199 inline V2f InvBilin::solve(V2f M1, V2f M2, V2f b)
  200 {
  201     float det = cross(M1, M2);
  202     if(unsafeInvert || det != 0) det = 1/det;
  203     return det * V2f(cross(b, M2), -cross(b, M1));
  204 }
  205 
  206 /// Evaluate the bilinear function at the coordinates (u,v)
  207 inline V2f InvBilin::bilinEval(V2f uv) const
  208 {
  209     return m_A + m_E*uv.x + m_F*uv.y + m_G*uv.x*uv.y;
  210 }
  211 
  212 
  213 } // namespace Aqsis
  214 
  215 #endif // AQSIS_INVBILIN_H_INCLUDED