---

gunu_interpolate.cpp

Go to the documentation of this file.

/* LIBGUL - Geometry Utility Library
 * Copyright (C) 1998-1999 Norbert Irmer
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 */

 
#include "stdafx.h"

#include <iostream>

#include "gul_types.h"
#include "gul_vector.h"
#include "gunu_basics.h"
#include "gunu_interpolate.h"

namespace gunu {

using gul::rtr;
using gul::Ptr;
using gul::point;
using gul::hpoint;
using gul::Max;
using gul::length;
using gul::cross_product;
using gul::set;

/*----------------------------------------------------------------------
  Calculates the control points of a cubic curve, interpolating
  the data points Q0,..,Qn. (see "The NURBS Book")
  The knot vector U must be computed by the caller.
  The first and last two control points P0,P1,Pn+1,Pn+2, depend
  on the directions & magnitudes of the curve endpoint derivatives,
  and must be filled in by the caller too.
 ----------------------------------------------------------------------*/
template< class T >
GULAPI void CubicCurveInterpolation( int n, Ptr<point<T> > Q, Ptr<T> U,
                                  Ptr<point<T> > P )
{
  Ptr<point<T> > R;
  Ptr<T>         dd,abc;
  T   deni;
  int i;

  R.reserve_pool(n+1);
  dd.reserve_pool(n+1);
  abc.reserve_pool(4);
 
  for( i = 3; i < n; i++ )
    R[i] = Q[i-1];

  CalcBasisFunctions( U[4], 4, 3, U, abc );

  deni = (T)1/abc[1];
  P[2] = deni*(Q[1]-abc[0]*P[1]);
 
  for( i = 3; i < n; i++ )
  {
    dd[i] = deni*abc[2];
    
    CalcBasisFunctions( U[i+2], i+2, 3, U, abc );

    deni = (T)1/(abc[1]-abc[0]*dd[i]);
    P[i] = deni*(R[i]-abc[0]*P[i-1]);
  }

  dd[n] = deni*abc[2];

  CalcBasisFunctions( U[n+2], n+2, 3, U, abc );

  deni = (T)1/(abc[1]-abc[0]*dd[n]);
  P[n] = deni*(Q[n-1]-abc[2]*P[n+1]-abc[0]*P[n-1]);

  for( i = n-1; i >= 2; i-- )
  {
    P[i] = P[i]-dd[i+1]*P[i+1];
  }
}
template
GULAPI void CubicCurveInterpolation( int n, Ptr<point<float> > Q, Ptr<float> U,
                                  Ptr<point<float> > P );
template
GULAPI void CubicCurveInterpolation( int n, Ptr<point<double> > Q, Ptr<double> U,
                                  Ptr<point<double> > P );

/*-----------------------------------------------------------------------
  Creates a non-homogeneous NURBS curve, which interpolates a set of 'n+1'
  data points 'Q' with "Local Bicubic Interpolation" (see "The NURBS Book").
  For each data point two control points are generated. Storage for
  the output arrays 'U' and 'P' must be reserved by the caller:

  For P: (n+1)*2
  For U: ((n+1)*2 + 4)
-------------------------------------------------------------------- */  
template< class K >
GULAPI void CubicLocalCurveInterpolation( 
                     int n, Ptr< point<K> > Q, bool cornerflag,
                     Ptr<K> U, Ptr< point<K> > P )
{
  Ptr< point<K> > q,T;
  point<K> v1;
  K a,b,c,alpha,l1,l2,u,ui;
  int upos, cpts_pos,k;
  bool no_alpha;
  
  q.reserve_pool(n+4);
  T.reserve_pool(n+1);
  
 /* ---------- calculate tangent vectors T{k} -----------------------------*/

// it's a bit confusing, but to be able to have a q[-1], i am using:
// q'[k] = Q{k} - Q{k-1}
// q[k+1] -> q'[k]

  for( k = 1; k <= n; k++ )
    q[k+1] = Q[k] - Q[k-1];

  q[1] = (K)2 * q[2] - q[3];        /* q'{0}  =  2 * q'{1} - q'{2} */
  q[0] = (K)2 * q[1] - q[2];        /* q'[-1} =  2 * q'{0} - q'{1} */

  q[n+2] = (K)2 * q[n+1] - q[n];    /* q'{n+1} = 2 * q'{n} - q'{n-1}   */
  q[n+3] = (K)2 * q[n+2] - q[n+1];  /* q'{n+2} = 2 * q'{n+1} - q'{n}   */
  
  for( k = 0; k <= n; k++ )
  {    
    /* alpha{k} = |q'{k-1} x q'{k}| / (|q'{k-1} x q'{k}| + |q'{k+1} x q'{k+2}|)
    */
    l1 = length(cross_product(q[k], q[k+1]));
    l2 = length(cross_product(q[k+2], q[k+3]));

    no_alpha = false;
    if( l1 + l2 != 0.0 )
      alpha = l1 / (l1 + l2);   
    else
      if( !cornerflag )
        alpha = 0.5;              /* make corners round */
      else
        no_alpha = true;            /* take 0 as tangent vector */
        
    if( no_alpha == false )
    {
                                    /* V{k} = (1-alpha)*q{k} + alpha*q{k+1} */
      T[k] = ((K)1-alpha)*q[k+1] + alpha*q[k+2];
      l1 = length(T[k]);
      if( l1 != 0.0 )
        T[k] = ((K)1/l1) * T[k];    /*  T{k} = V{k} / |V{k}| */
      else
        set( T[k], (K)0 );
    }
    else
     set( T[k], (K)0 );    
  }

/* ----------- calculate Bezier segments ------------------------------- */

  P[0] = Q[0];
  cpts_pos = 1;
  U[0] = 0;
  U[1] = 0;
  U[2] = 0;
  U[3] = 0;
  upos = 4;
  
  for( k = 0; k < n; k++ )
  {
    v1 = T[k] + T[k+1];       /* a = 16 - |T{k}{0} + T{k}{3}|^2 */
    a = ((K)16) - (v1 * v1);  /* (with T{k}{0} = T{k}, T{k}{3} = T{k+1} */
    
    l1 = q[k+2] * v1;         /* b = 12 * (P{3}-P{0})*(T{k}{0}+T{k}{3}) */
    b = ((K)12) * l1;         /* (with P{0} = Q{k}, P{3} = Q{k+1}) */

    c = ((K)-36) * (q[k+2] * q[k+2]);    /* c = -36 * |P{3} - P{0}|^2 */
    
    alpha = 
      ( -(b/(((K)2)*a)) + rtr<K>::sqrt(((b*b)/(((K)4)*a*a)) - (c/a)) ) / ((K)3);

                                        
    P[cpts_pos++] = Q[k] + (alpha * T[k]);     /* P1 = P0 + 1/3 * alpha' * T0 */    
    P[cpts_pos++] = Q[k+1] - (alpha * T[k+1]); /* P2 = P3 - 1/3 * alpha' * T3 */

                                /* u{k+1} = u{k} + 3 * |(P{k}{1} - P{k}{0})| */
    v1 = P[cpts_pos-2] - Q[k];
    l1 = length(v1);
    u = U[upos-1] + ((K)3) * l1;   
    U[upos++] = u;
    U[upos++] = u;
  }
 
  P[cpts_pos] = Q[n];  
  upos -= 2;
  ui = ((K)1) / U[upos];
  for( k = 4; k < upos; k++ )
    U[k] *= ui;

  U[upos++] = (K)1;
  U[upos++] = (K)1;
  U[upos++] = (K)1;
  U[upos] = (K)1;
}
template
GULAPI void CubicLocalCurveInterpolation( 
       int n, Ptr< point<float> > Q, bool cornerflag,
       Ptr<float> U, Ptr< point<float> > P );
template
GULAPI void CubicLocalCurveInterpolation( 
       int n, Ptr< point<double> > Q, bool cornerflag,
       Ptr<double> U, Ptr< point<double> > P );

/*-----------------------------------------------------------------------
  Creates a non-homogeneous NURBS surface, which interpolates a net of 
  'm+1' * 'n+1' data points 'Q' with "Local Bicubic Interpolation" 
  (see "The NURBS Book"). For each data point two control points are generated. Storage for
  the output arrays 'U','V' and 'controlPoints' must be reserved by the caller:

  For controlPoints: (m+1)*(n+1)*2*sizeof(NUPoint)  
  For U: ((n+1)*2 + 4)
  For V: ((m+1)*2 + 4)
-------------------------------------------------------------------- */  
template< class T, class HP >
GULAPI void CubicLocalSurfaceInterpolation( 
           int n, int m, Ptr< Ptr< point<T> > > Q, bool cornerflag,
           Ptr<T> U, Ptr<T> V, Ptr< Ptr< HP > > controlPoints )
{
  T gamma, alpha, a, l1, l2, total, d;
  int rpos, cpos, k, l, i;
  bool no_alpha;
  Ptr< Ptr< Ptr< Ptr< point<T> > > > > P;
  Ptr<T> r, s, alpha_k, beta_l, delta_uk, delta_vl, ub, vb;
  Ptr< Ptr< point<T> > > D_uv_, d_vu_, d_uv_, T_u_, T_v_;
  Ptr< point<T> > q, D_v_, d_v_, D_u_, d_u_;
  point<T> v1;
  
/* ---- reserve memory for local arrays ----------------------------- */

  P.reserve_pool(m+1);
  T_u_.reserve_pool(m+1);
  T_v_.reserve_pool(m+1);
  d_vu_.reserve_pool(m+1);
  d_uv_.reserve_pool(m+1);
  D_uv_.reserve_pool(m+1);

  for( l = 0; l <= m; l++ )
  {
    P[l].reserve_pool(n+1);
    T_u_[l].reserve_pool(n+1);
    T_v_[l].reserve_pool(n+1);
    d_vu_[l].reserve_pool(n+1);
    d_uv_[l].reserve_pool(n+1);
    D_uv_[l].reserve_pool(n+1);

    for( k = 0; k <= n; k++ )
    {
      P[l][k].reserve_pool(3);

      for( i = 0; i < 3; i++ )
        P[l][k][i].reserve_pool(3);
    }
  }   
  ub.reserve_pool(n+1);
  vb.reserve_pool(m+1);
  k = Max(m,n);
  q.reserve_pool(k+4);
  r.reserve_pool(m+1);
  s.reserve_pool(n+1);
  delta_vl.reserve_pool(m+1);
  delta_uk.reserve_pool(n+1);
  beta_l.reserve_pool(m+1);
  alpha_k.reserve_pool(n+1);
  D_v_.reserve_pool(n+1);
  d_v_.reserve_pool(n+1);
  D_u_.reserve_pool(m+1);
  d_u_.reserve_pool(m+1);

/* -------------- calculate T_u_{k,l} ------------------------------------- */  

  total = 0.0;
    
  for( k = 0; k <= n; k++ )
    ub[k] = 0.0;

  for( l = 0; l <= m; l++ )
  {
    for( k = 1; k <= n; k++ )
      q[k+1] = Q[l][k] - Q[l][k-1];

    q[1] = ((T)2) * q[2] - q[3];
    q[0] = ((T)2) * q[1] - q[2];
    q[n+2] = ((T)2) * q[n+1] - q[n];
    q[n+3] = ((T)2) * q[n+2] - q[n+1];
  
    for( k = 0; k <= n; k++ )
    {    
      l1 = length(cross_product(q[k], q[k+1]));
      l2 = length(cross_product(q[k+2], q[k+3]));

      no_alpha = false;
      if( l1 + l2 != 0.0 )
        alpha = l1 / (l1 + l2);   
      else
        if( !cornerflag )
          alpha = (T)0.5;
        else
          no_alpha = true;

      if( no_alpha == false )                
      {
        T_u_[l][k] = ((T)1-alpha)*q[k+1] + alpha*q[k+2];

        l1 = length(T_u_[l][k]);           
        if( l1 != 0.0 )
          T_u_[l][k] = ((T)1/l1) * T_u_[l][k];
        else
          set( T_u_[l][k], (T)0 );
      }
      else
        set( T_u_[l][k], (T)0 );
    }
/* ------------------------------------------------------------------------- */

    r[l] = (T)0;

    for( k = 1; k <= n; k++ )
    {
      d = length(q[k+1]);     /* d = |Q{l,k} - Q{l,k-1}| */
      ub[k] += d;
      r[l] += d;
    }  
    total += r[l];    /* add total chord length of this row */
  }  
  for( k = 1; k < n; k++ )
    ub[k] = ub[k-1] + (ub[k]/total);    

  ub[n] = 1.0;

/* -------------- calculate T_v_{k,l} ------------------------------------- */  

  total = 0.0;
    
  for( l = 0; l <= m; l++ )
    vb[l] = 0.0;

  for( k = 0; k <= n; k++ )
  {
    for( l = 1; l <= m; l++ )
      q[l+1] = Q[l][k] - Q[l-1][k];
 
    q[1] = ((T)2) * q[2] - q[3];
    q[0] = ((T)2) * q[1] - q[2];
    q[m+2] = ((T)2) * q[m+1] - q[m];
    q[m+3] = ((T)2) * q[m+2] - q[m+1];

    for( l = 0; l <= m; l++ )
    {    
      l1 = length(cross_product(q[l], q[l+1]));
      l2 = length(cross_product(q[l+2], q[l+3]));

      no_alpha = false;
      if( l1 + l2 != (T)0 )
        alpha = l1 / (l1 + l2);   
      else
        if( !cornerflag )
          alpha = (T)0.5;              /* make sharp corners round */
        else
          no_alpha = true;            /* else set tangent vector to 0 */
        
      if( no_alpha == false )
      {
        T_v_[l][k] = ((T)1-alpha) * q[l+1] + alpha*q[l+2];
        l1 = length( T_v_[l][k] );            
        if( l1 != (T)0 )
          T_v_[l][k] = ((T)1/l1) * T_v_[l][k];
        else
          set( T_v_[l][k], (T)0 );
      }
      else
        set( T_v_[l][k], (T)0 );
    }
/* ------------------------------------------------------------------------- */

    s[k] = 0.0;

    for( l = 1; l <= m; l++ )
    {
      d = length( q[l+1] );     /* d = |Q{l,k} - Q{l-1,k}| */
      vb[l] += d;
      s[k] += d;
    }  
    total += s[k]; /* add total chord length of this column */ 
  }  
  for( l = 1; l < m; l++ )
    vb[l] = vb[l-1] + (vb[l]/total);    

  vb[m] = (T)1;

/* -------------------------------------------------------------------  
   Calculate the control points of the Bezier patches which lie in the rows
   and columns containing the data points

  Formula for rows:
   P{l,k}{0,0] = Q{l,k}
   P{l,k}{0,1} = Q{l,k} + a * T_u_{l,k} 
   P{l,k}{0,2} = Q{l,k+1} - a * T_u_{l,k+1} 
  with
   a = r{l} * (ub{k+1} - ub{k}) / 3

  Formula for columns:
   P{l,k}{1,0} = Q{l,k} + a * T_v_{l,k} 
   P{l,k}{2,0} = Q{l+1,k} - a * T_v_{l+1,k} 
  with
   a = s{k} * (vb{l+1} - vb{l}) / 3
------------------------------------------------------------------------ */
  for( l = 0; l <= m; l++ )
  {
    for( k = 0; k < n; k++ )
    { 
      P[l][k][0][0] = Q[l][k];
      a = (r[l] * (ub[k+1] - ub[k])) / (T)3;
      P[l][k][0][1] = Q[l][k] + a * T_u_[l][k];
      P[l][k][0][2] = Q[l][k+1] - a * T_u_[l][k+1];
    }
    P[l][n][0][0] = Q[l][n];
  }
  
  for( k = 0; k <= n; k++ )
  {
    for( l = 0; l < m; l++ )
    {
      a = (s[k] * (vb[l+1] - vb[l])) / (T)3;
      P[l][k][1][0] = Q[l][k] + a * T_v_[l][k];
      P[l][k][2][0] = Q[l+1][k] - a * T_v_[l+1][k];
    }
  } 
/* ----------------------------------------------------------------------
 Calculate the 4 inner control points of the Bezier patches.
 For this the mixed partial derivatives D_uv_{k,l} at the 4 corners 
 of the patch are needed
----------------------------------------------------------------------- */
  for( k = 1; k <= n; k++ )
    delta_uk[k] = ub[k] - ub[k-1];
  for( k = 1; k < n; k++ )
    alpha_k[k] = delta_uk[k] / (delta_uk[k] + delta_uk[k+1]);
  alpha_k[0] = 0.5;
  alpha_k[n] = 0.5;  

  for( l = 1; l <= m; l++ )
    delta_vl[l] = vb[l] - vb[l-1];
  for( l = 1; l < m; l++ )
    beta_l[l] = delta_vl[l] / (delta_vl[l] + delta_vl[l+1]);  
  beta_l[0] = 0.5;
  beta_l[m] = 0.5;
   
/* --------- calculate d_vu_{k,l} --------------------------------- */
    
  for( l = 0; l <= m; l++ )
  {
    for( k = 0; k <= n; k++ )
      D_v_[k] = s[k] * T_v_[l][k];    /* D_v_{l,k} = s{k] * T_v_{l,k} */

    for( k = 1; k <= n; k++ )
    {
                      /* d_v_{l,k} := (D_v_{l,k} - D_v_{l,k-1}) / delta_u{k} */
      d_v_[k] = D_v_[k] - D_v_[k-1];

      if( delta_uk[k] != (T)0 )
        d_v_[k] = ((T)1 / delta_uk[k]) * d_v_[k];
      else
        set( d_v_[k], (T)0 );
    }       
    for( k = 1; k < n; k++ )
    {
      d_vu_[l][k] = ((T)1 - alpha_k[k]) * d_v_[k] + alpha * d_v_[k+1];
    }
/* special cases at the borders (k=0,k=n): */

    d_vu_[l][0] = (T)2 * d_v_[1] - d_vu_[l][1];
    d_vu_[l][n] = (T)2 * d_v_[n] - d_vu_[l][n-1];
  }    
/* --------- calculate d_uv_{k,l} --------------------------------- */
    
  for( k = 0; k <= n; k++ )
  {
    for( l = 0; l <= m; l++ )
      D_u_[l] = r[l] * T_u_[l][k];    /* D_u_{l,k} = r{l] * T_u_{l,k} */

    for( l = 1; l <= m; l++ )
    {
                      /* d_u_{l,k} := (D_u_{l,k} - D_u_{l-1,k}) / delta_v{l} */
      d_u_[l] = D_u_[l] - D_u_[l-1];
      if( delta_vl[l] != (T)0 )
        d_u_[l] = ((T)1 / delta_vl[l]) * d_u_[l];
      else
        set( d_u_[l], (T)0 );
    }  
    for( l = 1; l < m; l++ )
    {
      d_uv_[l][k] = ((T)1 - beta_l[l]) * d_u_[l] + beta_l[l] * d_u_[l+1];
    }
/* special cases at the borders (l=0,l=n): */

    d_uv_[0][k] = (T)2 * d_u_[1] - d_uv_[1][k];
    d_uv_[m][k] = (T)2 * d_u_[m] - d_uv_[m-1][k];
  }    
/* --- calculate the mixed partial derivatives D_uv_{k,l} ------ */  

  for( l = 0; l <= m; l++ )
    for( k = 0; k <= n; k++ )
    {
      v1 = alpha_k[k] * d_uv_[l][k] + beta_l[l] * d_vu_[l][k];
      if( alpha_k[k] + beta_l[l] != (T)0 )
        D_uv_[l][k] = ((T)1/(alpha_k[k] + beta_l[l])) * v1; 
      else
        set( D_uv_[l][k], (T)0 );
    }
/*--- calculate the 4 inner control points of the Bezier patches --------- */
/*  used formulas:
  P{l,k}{1,1} = (gamma * D_uv_{l,k}) + P{l,k}{1,0} + P{l,k}{0,1} - P{l,k}{0,0}      
  ... (see NURBS Book, page 403)
*/
  for( l = 0; l < m; l++ )
    for( k = 0; k < n; k++ )
    {
      gamma = (delta_uk[k+1] * delta_vl[l+1]) / (T)9;

      P[l][k][1][1] = gamma * D_uv_[l][k] + 
                      P[l][k][1][0] + P[l][k][0][1] - P[l][k][0][0];
      P[l][k][1][2] = -gamma * D_uv_[l][k+1] + 
                      P[l][k+1][1][0] - P[l][k+1][0][0] + P[l][k][0][2];
      P[l][k][2][1] = -gamma * D_uv_[l+1][k] + 
                      P[l+1][k][0][1] - P[l+1][k][0][0] + P[l][k][2][0];
      P[l][k][2][2] = gamma * D_uv_[l+1][k+1] + 
                      P[l+1][k][0][2] + P[l][k+1][2][0] - P[l+1][k+1][0][0];
    } 
/* ---- construct the knot vectors U,V ---------------------------- */

  U[0] = 0.0;
  U[1] = 0.0;
  U[2] = 0.0;
  U[3] = 0.0;
  cpos = 4;
  for( k = 1; k < n; k++ )
  {
    U[cpos++] = ub[k];
    U[cpos++] = ub[k];
  }   
  U[cpos++] = 1.0;
  U[cpos++] = 1.0;
  U[cpos++] = 1.0;
  U[cpos] = 1.0;

  V[0] = 0.0;
  V[1] = 0.0;
  V[2] = 0.0;
  V[3] = 0.0;
  cpos = 4;
  for( l = 1; l < m; l++ )
  {
    V[cpos++] = vb[l];
    V[cpos++] = vb[l];
  }   
  V[cpos++] = 1.0;
  V[cpos++] = 1.0;
  V[cpos++] = 1.0;
  V[cpos] = 1.0;

/* ---------- fill control point array ---------------- */ 

  rpos = 0;
  cpos = 0;

/* ------------- first row ------------------------------------------- */

  controlPoints[rpos][cpos++] = P[0][0][0][0];
  for( k = 0; k < n; k++ )
  {
    controlPoints[rpos][cpos++] = P[0][k][0][1];
    controlPoints[rpos][cpos++] = P[0][k][0][2];
  }
  controlPoints[rpos][cpos++] = P[0][n][0][0];

/* -------------- inner rows ------------------------------------ */

  for( l = 0; l < m; l++ )
  {
    rpos++;
    cpos = 0;    
    controlPoints[rpos][cpos++] = P[l][0][1][0];
    for( k = 0; k < n; k++ )
    {
      controlPoints[rpos][cpos++] = P[l][k][1][1];
      controlPoints[rpos][cpos++] = P[l][k][1][2];
    }
    controlPoints[rpos][cpos++] = P[l][n][1][0];
    
    rpos++;
    cpos = 0;
    controlPoints[rpos][cpos++] = P[l][0][2][0];
    for( k = 0; k < n; k++ )
    {
      controlPoints[rpos][cpos++] = P[l][k][2][1];
      controlPoints[rpos][cpos++] = P[l][k][2][2];
    }
    controlPoints[rpos][cpos++] = P[l][n][2][0];
  }
/* --------------- last row --------------------------------------- */

  rpos++;
  cpos = 0;
  controlPoints[rpos][cpos++] = P[m][0][0][0];
  for( k = 0; k < n; k++ )
  {
    controlPoints[rpos][cpos++] = P[m][k][0][1];
    controlPoints[rpos][cpos++] = P[m][k][0][2];
  }
  controlPoints[rpos][cpos++] = P[m][n][0][0];

/*------------- ready ------------------------------------------------- */
}
template GULAPI void CubicLocalSurfaceInterpolation( 
    int n, int m, Ptr< Ptr< point<float> > > Q, bool cornerflag,
    Ptr<float> U, Ptr<float> V, Ptr< Ptr< point<float> > > controlPoints );
template GULAPI void CubicLocalSurfaceInterpolation( 
    int n, int m, Ptr< Ptr< point<double> > > Q, bool cornerflag,
    Ptr<double> U, Ptr<double> V, Ptr< Ptr< point<double> > > controlPoints );
template GULAPI void CubicLocalSurfaceInterpolation( 
    int n, int m, Ptr< Ptr< point<float> > > Q, bool cornerflag,
    Ptr<float> U, Ptr<float> V, Ptr< Ptr< hpoint<float> > > controlPoints );
template GULAPI void CubicLocalSurfaceInterpolation( 
    int n, int m, Ptr< Ptr< point<double> > > Q, bool cornerflag,
    Ptr<double> U, Ptr<double> V, Ptr< Ptr< hpoint<double> > > controlPoints );
}

---