---

gunu_mba_approximate.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 "gul_matrix.h"
#include "guge_normalize.h"
#include "guma_minimize.h"
#include "gunu_basics.h"
#include "gunu_refine.h"
#include "gunu_make_compatible.h"
#include "gunu_mba_approximate.h"

namespace gunu {

using gul::rtr;
using gul::Max;
using gul::Ptr;
using gul::point;
using gul::point1;
using gul::point2;
using gul::vec4;
using gul::mat4x4;
using guge::CalcBoundingBoxE;
using guma::GoldenSectionSearch;
using gul::set;

/*-------------------------------------------------------------------------*//**
  BSpline approximation,
  calculates the points of a control lattice 'delta', which approximates
  data points in 'P'. (P[].z contains the function value, P[].x and P[].y
  the location in the domain)                                                 */
/*----------------------------------------------------------------------------*/
template< class T >
void BAapproximate( 
             int nP, Ptr< point<T> >& P, 
             int nu, int pu, Ptr<T>& U,  
             int nv, int pv, Ptr<T>& V,
             Ptr< Ptr< point1<T> > >& delta )
{
  Ptr<T> Nu, Nv;
  Ptr< Ptr<T> > omega, w, w2;
  int i,j,k,l,C,uspan,vspan,uind,vind; 
  T sum_w2;

/* ---- reserve local arrays ------------------- */ 

  omega.reserve_pool(nv+1);
  for( i = 0; i < nv+1; i++ )
    omega[i].reserve_pool(nu+1);

  w.reserve_pool(pu+1);
  w2.reserve_pool(pu+1);
  for( i = 0; i < pu+1; i++ )
  {
    w[i].reserve_pool(pv+1);
    w2[i].reserve_pool(pv+1);
  }
  Nu.reserve_pool(pu+1);
  Nv.reserve_pool(pv+1);
  
  for( j = 0; j <= nv; j++ )
  {
    for( i = 0; i <= nu; i++ )
    {
      delta[j][i].x = (T)0;
      omega[j][i] = (T)0;
    }  
  }          
  for( C = 0; C < nP; C++ )
  {
    uspan = FindSpan( P[C].x, nu, pu, U );
    CalcBasisFunctions( P[C].x, uspan, pu, U, Nu );

    vspan = FindSpan( P[C].y, nv, pv, V );
    CalcBasisFunctions( P[C].y, vspan, pv, V, Nv );

    sum_w2 = (T)0;
    for( k = 0; k <= pu; k++ ) 
    {
      for( l = 0; l <= pv; l++ )
      {
        w[k][l] = Nu[k] * Nv[l];
        w2[k][l] = w[k][l] * w[k][l];
        sum_w2 += w2[k][l];      
      }
    }
    uind = uspan - pu;
    vind = vspan - pv;
    for( k = 0; k <= pu; k++ )
    {
      for( l = 0; l <= pv; l++ )
      {
        delta[vind+l][uind+k].x += w2[k][l] * w[k][l] * P[C].z / sum_w2;
        omega[vind+l][uind+k] += w2[k][l];
      }
    }
  }
  for( j = 0; j <= nv; j++ )
  {
    for( i = 0; i <= nu; i++ )
    {
      if( omega[j][i] != (T)0 )
        delta[j][i].x /= omega[j][i];
      else
        delta[j][i].x = (T)0;
    }
  }                 
}
template void BAapproximate( 
             int nP, Ptr< point<float> >& P, 
             int nu, int pu, Ptr<float>& U,  
             int nv, int pv, Ptr<float>& V,
             Ptr< Ptr< point1<float> > >& delta );
template void BAapproximate( 
             int nP, Ptr< point<double> >& P, 
             int nu, int pu, Ptr<double>& U,  
             int nv, int pv, Ptr<double>& V,
             Ptr< Ptr< point1<double> > >& delta );

/*-------------------------------------------------------------------------*//**
  BSpline approximation,
  calculates the points of a control lattice 'delta', which approximates
  data points in 'P'. (P[].z contains the function value, P[].x and P[].y
  the location in the domain). An individual standard deviation for each data
  point must be given in 'Sigma'                                              */
/*----------------------------------------------------------------------------*/
template< class T >
void BAapproximate( 
             int nP, Ptr< point<T> >& P, Ptr<T>& Sigma,
             int nu, int pu, Ptr<T>& U,  
             int nv, int pv, Ptr<T>& V,
             Ptr< Ptr< point1<T> > >& delta )
{
  Ptr<T> Nu, Nv;
  Ptr< Ptr<T> > omega, w, w2;
  int i,j,k,l,C,uspan,vspan,uind,vind; 
  T factor, sum_w2;

/* ---- reserve local arrays ------------------- */ 

  omega.reserve_pool(nv+1);
  for( i = 0; i < nv+1; i++ )
    omega[i].reserve_pool(nu+1);

  w.reserve_pool(pu+1);
  w2.reserve_pool(pu+1);
  for( i = 0; i < pu+1; i++ )
  {
    w[i].reserve_pool(pv+1);
    w2[i].reserve_pool(pv+1);
  }
  Nu.reserve_pool(pu+1);
  Nv.reserve_pool(pv+1);
  
  for( j = 0; j <= nv; j++ )
  {
    for( i = 0; i <= nu; i++ )
    {
      delta[j][i].x = (T)0;
      omega[j][i] = (T)0;
    }  
  }          
  for( C = 0; C < nP; C++ )
  {
    uspan = FindSpan( P[C].x, nu, pu, U );
    CalcBasisFunctions( P[C].x, uspan, pu, U, Nu );

    vspan = FindSpan( P[C].y, nv, pv, V );
    CalcBasisFunctions( P[C].y, vspan, pv, V, Nv );

    sum_w2 = (T)0;
    for( k = 0; k <= pu; k++ ) 
    {
      for( l = 0; l <= pv; l++ )
      {
        w[k][l] = Nu[k] * Nv[l];
        w2[k][l] = w[k][l] * w[k][l];
        sum_w2 += w2[k][l];      
      }
    }
    factor = (T)1 / (Sigma[C] * Sigma[C]); 
      
    uind = uspan - pu;
    vind = vspan - pv;
    for( k = 0; k <= pu; k++ )
    {
      for( l = 0; l <= pv; l++ )
      {
        w2[k][l] *= factor;             /* (weight/standard_deviation)^2 */
        
        delta[vind+l][uind+k].x += w2[k][l] * w[k][l] * P[C].z / sum_w2;
        omega[vind+l][uind+k] += w2[k][l];
      }
    }
  }
  for( j = 0; j <= nv; j++ )
  {
    for( i = 0; i <= nu; i++ )
    {
      if( omega[j][i] != (T)0 )
        delta[j][i].x /= omega[j][i];
      else
        delta[j][i].x = (T)0;
    }
  }                 
}
template void BAapproximate( 
             int nP, Ptr< point<float> >& P, Ptr<float>& Sigma,
             int nu, int pu, Ptr<float>& U,  
             int nv, int pv, Ptr<float>& V,
             Ptr< Ptr< point1<float> > >& delta );
template void BAapproximate( 
             int nP, Ptr< point<double> >& P, Ptr<double>& Sigma,
             int nu, int pu, Ptr<double>& U,  
             int nv, int pv, Ptr<double>& V,
             Ptr< Ptr< point1<double> > >& delta );


/*-------------------------------------------------------------------------*//**
  Multilevel BSpline approximation,
  Calculates a NURBS surface (one-dimensional), which approximates
  data points in 'P', (P[].z contains the function value, P[].x and P[] the
  location in the domain). An individual standard deviation for each data
  point can be given in 'Sigma'.
  (Remarks: the data point array P is changed - after the algorithm the P[].z
  contain the difference between the function values of the calculated
  surface and and the original value!!!,
  Output arrays U,V,Psi must be reserved by the caller:
  U: 2*pu + 2^(nIter-1) + 1,
  V: 2*pv + 2^(nIter-1) + 1,
  Psi: (pv + 2^(nIter-1))*(pu + 2^(nIter-1))                                  */   
/*----------------------------------------------------------------------------*/
template< class T >
void MBAapproximate( 
       int nP, Ptr< point<T> >& P, bool useSigma, Ptr<T>& Sigma, int nIter,
       int pu, int pv, Ptr<T>& U, Ptr<T>& V, Ptr< Ptr< point1<T> > >& Psi )                      
{
  Ptr<T> U1,V1,X;
  Ptr< Ptr< point1<T> > > Phi;
  int xp,i,j,nu,nv,s,maxSpans;
  point1<T> S;
  T d;

  if( nIter < 1 ) return;
  maxSpans = 1 << (nIter-1);

/* --- reserve memory for local arrays ----------------------------- */
  
  U1.reserve_pool(pu + maxSpans + pu + 1);
  V1.reserve_pool(pv + maxSpans + pv + 1);
  i = Max(pu,pv);
  X.reserve_pool(i + maxSpans + i + 1);

  Phi.reserve_pool(pv + maxSpans);
  for( i = 0; i < maxSpans + pv; i++ )
    Phi[i].reserve_pool(pu + maxSpans);

/* ----------------------------------------------------------------- */

  for( i = 0; i <= pu; i++ )
  {
    U[i] = (T)0;
    U[pu+1+i] = (T)1;
  }
  for( i = 0; i <= pv; i++ )
  {
    V[i] = (T)0;
    V[pv+1+i] = (T)1;
  }  
  
  for( j = 0; j <= pv; j++ )
    for( i = 0; i <= pu; i++ )
      set( Psi[j][i], (T)0 );

  s = 1;
  nu = pu;
  nv = pv;

  do
  {
    if( useSigma )
       BAapproximate( nP, P, Sigma, nu, pu, U, nv, pv, V, Phi );
    else
       BAapproximate( nP, P, nu, pu, U, nv, pv, V, Phi );

    for( j = 0; j <= nv; j++ )
      for( i = 0; i <= nu; i++ )
        Psi[j][i].x += Phi[j][i].x;
       
    for( i = 0; i < nP; i++ )
    {
      SurfacePoint( P[i].x, P[i].y, nu, pu, U, nv, pv, V, Phi, &S );               
      P[i].z -= S.x;
    }

    if( s >= maxSpans )          /* --- ready --- */
      break;         

    s <<= 1;
    
    d = (T)1/(T)s;
    for( i = pu+1, xp = 0; i <= nu+1; i++, xp++ ) 
      X[xp] = U[i] - d;
    RefineSurface( nu, pu, U, nv, pv, V, Psi, 
                   X, xp-1, gul::u_direction, U1, V1, Phi );      
    nu += xp;
    
    for( i = pv+1, xp = 0; i <= nv+1; i++, xp++ ) 
      X[xp] = V[i] - d;
    RefineSurface( nu, pu, U1, nv, pv, V1, Phi,
                   X, xp-1, gul::v_direction, U, V, Psi );
    nv += xp;
      
  }while(1);
}
template void MBAapproximate( 
       int nP, Ptr< point<float> >& P, bool useSigma, Ptr<float>& Sigma,
       int nIter, int pu, int pv, Ptr<float>& U, Ptr<float>& V,
       Ptr< Ptr< point1<float> > >& Psi );
template void MBAapproximate( 
       int nP, Ptr< point<double> >& P, bool useSigma, Ptr<double>& Sigma,
       int nIter, int pu, int pv, Ptr<double>& U, Ptr<double>& V,
       Ptr< Ptr< point1<double> > >& Psi );


/*------------------------------------------------------------------------*//**
  utility class for finding a rotation, so that the volume of the bounding 
  box of a set of points gets minimal volume. this can be done much more
  optimal (by using RAPID's method for calculating oriented bounding boxes)! */
/*---------------------------------------------------------------------------*/
template< class T >
class minbbox_rec
{ 
public:
  Ptr< point<T> > P;
  int             nP;
  T theta0_x;
  T theta0_y;
  T theta0_z;

  minbbox_rec( int anP, Ptr< point<T> >& aP, T theta0x, T theta0y, T theta0z )
  {
    nP = anP; P = aP; 
    theta0_x = theta0x; theta0_y = theta0y; theta0_z = theta0z;
  }
  T volume( T theta_x, T theta_y, T theta_z )
  {
    mat4x4<T> m1 = mat4x4<T>::rotate_x(theta0_x + theta_x);
    mat4x4<T> m2 = mat4x4<T>::rotate_y(theta0_y + theta_y);
    mat4x4<T> m3 = mat4x4<T>::rotate_z(theta0_z + theta_z);
    mat4x4<T> m;
    m = m1 * m2 * m3;
    vec4<T> v,vt;
    T minx,maxx,miny,maxy,minz,maxz;

    v[0] = P[0].x;  v[1] = P[0].y;  v[2] = P[0].z;  v[3] = 1.0;
    vt = v * m;
    minx = maxx = vt[0];
    miny = maxy = vt[1];
    minz = maxz = vt[2];

    for( int i = 1; i < nP; i++ )
    {
      v[0] = P[i].x;  v[1] = P[i].y;  v[2] = P[i].z;  v[3] = 1.0;
      vt = v * m;
      if( vt[0] < minx ) minx = vt[0];
      else if( vt[0] > maxx ) maxx = vt[0];
      if( vt[1] < miny ) miny = vt[1];
      else if( vt[1] > maxy ) maxy = vt[1];
      if( vt[2] < minz ) minz = vt[2];
      else if( vt[2] > maxz ) maxz = vt[2];
    }
    return (maxx-minx)*(maxy-miny)*(maxz-minz);  /* volume */  
  }
  static T minimize_about_z_cb( T theta, void *data )
  {
    minbbox_rec *m = (minbbox_rec *)data;
    return m->volume( (T)0, (T)0, theta );
  }
};
template class minbbox_rec<float>;
template class minbbox_rec<double>;


/*-----------------------------------------------------------------------*//**
  calculate the rotation about the Z-axis, so that the bounding box
  of the given Points 'P' gets minimal volume. 'mintol' is the minimum change
  in the domaain between two iterations (in degree), and 'maxits' the maximum
  number of iterations                                                      */
/*------------------------------------------------------------------------- */
template< class T >
T MinimizeBBoxAboutZ( int nP, Ptr< point<T> >& P, T mintol, int maxits )
{
  T minx,maxx,miny,maxy,minz,maxz,theta,V,V1;
  minbbox_rec<T> mr( nP, P, (T)0, (T)0, (T)0 );
  
  CalcBoundingBoxE( nP, P, minx, maxx, miny, maxy, minz, maxz );
  V = (maxx-minx)*(maxy-miny)*(maxz-minz);  /* volume at 0 and 90 degree */
                                            /* between these: */
  V1 = mr.volume( (T)0, (T)0, rtr<T>::golden_c() * rtr<T>::pi()/(T)2 );  
  if( V1 < V )  
  {
    GoldenSectionSearch( 
       (T)0, rtr<T>::golden_c() * rtr<T>::pi()/(T)2, rtr<T>::pi()/(T)2,
       minbbox_rec<T>::minimize_about_z_cb, (void *)&mr, 
       mintol*rtr<T>::pi()/(T)180, maxits, &theta,
       &V );
  }
  else
    theta = (T)0;

  return theta;
}  
template float MinimizeBBoxAboutZ( int nP, Ptr< point<float> >& P, 
                                   float mintol, int maxits );
template double MinimizeBBoxAboutZ( int nP, Ptr< point<double> >& P,
                                    double mintol, int maxits );


/*-------------------------------------------------------------------*//**
  executes the MBA algorithm and constructs a 3-dimensional surface
  from the results (x(),y() are identity mappings).
  when 'minimize' is set, the base rectangle in the xy-plane is
  chosen so that its area is minimal                                    */
/* ---------------------------------------------------------------------*/
template< class T, class HP >
GULAPI void SurfaceOverXYPlane( 
          int nDatPoints, Ptr< point<T> >& datPoints,
          bool useStdDevs, Ptr<T>& StdDevs,
          bool minimize, int nIter,
          int pu, int pv,
          int *ret_nu, Ptr<T> *retU,
          int *ret_nv, Ptr<T> *retV,
          Ptr< Ptr<HP> > *retPw )
{

  T minx,maxx,miny,maxy,minz,maxz;
  Ptr< point<T> > mbaPoints;
  int ref_nu, ref_nv, mba_nu, mba_nv, nu, nv;
  int ref_pu, ref_pv, mba_pu, mba_pv, i, j; 
  Ptr<T> refU, refV, mbaU, mbaV, U, V;
  Ptr< Ptr< point<T> > > refP, tmp1P, P;
  Ptr< Ptr< point1<T> > > mbaPz, tmp2P;
  Ptr<T> tmp1U, tmp1V, tmp2U, tmp2V;
  int tmp1_nu, tmp1_nv, tmp2_nu, tmp2_nv;
  int tmp1_pu, tmp1_pv, tmp2_pu, tmp2_pv;
  T x,y,z,ox,oy,sint,cost,theta;
  Ptr< Ptr<HP> > Pw;
  point<T> CP;

  mbaPoints.reserve_pool(nDatPoints);

  if( minimize != 0 )
  {
    /* --- minimize area of base rectangle in xy-plane --------- */ 

    theta = MinimizeBBoxAboutZ( nDatPoints, datPoints, (T)1, 20 );

    sint = rtr<T>::sin(theta);
    cost = rtr<T>::cos(theta);

    // 1. point:
    ox = datPoints[0].x;
    oy = datPoints[0].y;
    z = datPoints[0].z;
    minz = maxz = z;
    minx = maxx = x = cost * ox + sint * oy;
    miny = maxy = y = cost * oy - sint * ox;
    mbaPoints[0].x = x;
    mbaPoints[0].y = y;
    mbaPoints[0].z = z;         
      
    // remaining points:      
    for( i = 1; i < nDatPoints; i++ )
    {
      ox = datPoints[i].x;
      oy = datPoints[i].y;
      z = datPoints[i].z;
      x = cost * ox + sint * oy;
      y = cost * oy - sint * ox;

      if( x < minx ) minx = x;
      else if( x > maxx ) maxx = x;
      if( y < miny ) miny = y;
      else if( y > maxy ) maxy = y;
      if( z < minz ) minz = z;
      else if( z > maxz ) maxz = z;

      mbaPoints[i].x = x;
      mbaPoints[i].y = y;
      mbaPoints[i].z = z;         
    }
    // normalise coordinates:
    for( i = 0; i < nDatPoints; i++ )
    {
      mbaPoints[i].x = (mbaPoints[i].x - minx) / (maxx-minx);
      if( mbaPoints[i].x < (T)0 ) mbaPoints[i].x = (T)0;
      else if( mbaPoints[i].x > (T)1 ) mbaPoints[i].x = (T)1;
      mbaPoints[i].y = (mbaPoints[i].y - miny) / (maxy-miny);
      if( mbaPoints[i].y < (T)0 ) mbaPoints[i].y = (T)0;
      else if( mbaPoints[i].y > (T)1 ) mbaPoints[i].y = (T)1;
      mbaPoints[i].z = (mbaPoints[i].z - minz) / (maxz-minz);
    } 
  }
  else
  {
    /* --- do not minimize area of base rectangle in xy-plane --------- */ 

    CalcBoundingBoxE( nDatPoints, datPoints,
                      minx, maxx, miny, maxy, minz, maxz );
    // normalise coordinates:
    for( i = 0; i < nDatPoints; i++ )
    {
      mbaPoints[i].x = (datPoints[i].x - minx) / (maxx-minx);
      if( mbaPoints[i].x < (T)0 ) mbaPoints[i].x = (T)0;
      else if( mbaPoints[i].x > (T)1 ) mbaPoints[i].x = (T)1;
      mbaPoints[i].y = (datPoints[i].y - miny) / (maxy-miny);
      if( mbaPoints[i].y < (T)0 ) mbaPoints[i].y = (T)0;
      else if( mbaPoints[i].y > (T)1 ) mbaPoints[i].y = (T)1;
      mbaPoints[i].z = (datPoints[i].z - minz) / (maxz-minz);
    }
  }  
  /* ------ do the MBA-Algorithm: ------------------------------------ */

  // reserve return arrays
  mba_nu = pu + (1<<(nIter-1)) - 1;
  mbaU.reserve_pool( 2*pu + (1<<(nIter-1)) + 1 );
  mba_nv = pv + (1<<(nIter-1)) - 1;
  mbaV.reserve_pool( 2*pv + (1<<(nIter-1)) + 1 );
  mbaPz.reserve_pool(mba_nv+1);
  Pw.reserve_pool(mba_nv+1);
  for( i = 0; i < mba_nv+1; i++ )
  {
    mbaPz[i].reserve_pool(mba_nu+1);
    Pw[i].reserve_pool(mba_nu+1);
  }             
  MBAapproximate<T>( nDatPoints, mbaPoints, useStdDevs, StdDevs,
                     nIter, pu, pv, mbaU, mbaV, mbaPz ); 

  /* ---- construct a 3d return surface: --------------------------------- */

  // build a bilinear basis surface
  refU.reserve_pool(4);
  for( i = 0; i < 2; i++ )
  {
    refU[i] = (T)0;
    refU[i+2] = (T)1;
  }     
  refV = refU;
  ref_nu = 1;
  ref_pu = 1;
  ref_nv = 1;
  ref_pv = 1;

  refP.reserve_pool(2);
  refP[0].reserve_pool(2);
  refP[1].reserve_pool(2);

  refP[0][0].x = (T)0; refP[0][0].y = (T)0; refP[0][0].z = (T)0;
  refP[0][1].x = (T)1; refP[0][1].y = (T)0; refP[0][1].z = (T)0;
  refP[1][0].x = (T)0; refP[1][0].y = (T)1; refP[1][0].z = (T)0;
  refP[1][1].x = (T)1; refP[1][1].y = (T)1; refP[1][1].z = (T)0;

  /* --- make base surface and surface created by MBA compatible ------ */

  MakeSurfacesCompatible<T,point1<T>,point<T> >( 
              mba_nu, pu, mbaU,  mba_nv, pv, mbaV, mbaPz, 
              ref_nu, ref_pu, refU, ref_nv, ref_pv, refV, refP,
              gul::v_direction,
              &tmp2_nu, &tmp2_pu, &tmp2U, &tmp2_nv, &tmp2_pv, &tmp2V, &tmp2P,
              &tmp1_nu, &tmp1_pu, &tmp1U, &tmp1_nv, &tmp1_pv, &tmp1V, &tmp1P );

  MakeSurfacesCompatible<T,point1<T>,point<T> >( 
              tmp2_nu, tmp2_pu, tmp2U, tmp2_nv, tmp2_pv, tmp2V, tmp2P,
              tmp1_nu, tmp1_pu, tmp1U, tmp1_nv, tmp1_pv, tmp1V, tmp1P,
              gul::u_direction,
              &mba_nu, &mba_pu, &mbaU, &mba_nv, &mba_pv, &mbaV, &mbaPz,
              &nu, &pu, &U, &nv, &pv, &V, &P );

  /* set z-coordinate of result surface control points: */

  for( j = 0; j <= nv; j++ )
    for( i = 0; i <= nu; i++ ) 
      P[j][i].z = mbaPz[j][i].x;

  /* denormalise control points: */

  for( j = 0; j <= nv; j++ )
    for( i = 0; i <= nu; i++ ) 
    {
      CP.x = (P[j][i].x * (maxx-minx)) + minx;
      CP.y = (P[j][i].y * (maxy-miny)) + miny;
      CP.z = (P[j][i].z * (maxz-minz)) + minz;
      Pw[j][i] = CP;
   }  

  if( minimize != 0 )
  {
    /* --- rotate surface to the original orientation ------------------ */

    sint = rtr<T>::sin(-theta);
    cost = rtr<T>::cos(-theta);

    for( j = 0; j <= nv; j++ )
      for( i = 0; i <= nu; i++ ) 
      {
        ox = Pw[j][i].x;
        oy = Pw[j][i].y; 

        x = cost * ox + sint * oy;
        y = cost * oy - sint * ox;

        Pw[j][i].x = x;
        Pw[j][i].y = y;        
      }  
  }        
    
  /* ready: */

  *ret_nu = nu;
  *retU = U;
  *ret_nv = nv;
  *retV = V;
  *retPw = Pw;
}
template GULAPI void SurfaceOverXYPlane( 
          int nDatPoints, Ptr< point<float> >& datPoints,
          bool useStdDevs, Ptr<float>& StdDevs,
          bool minimize, int nIterations,
          int pu, int pv,
          int *ret_nu, Ptr<float> *retU,
          int *ret_nv, Ptr<float> *retV,
          Ptr< Ptr< point<float> > > *retPw );
template GULAPI void SurfaceOverXYPlane( 
          int nDatPoints, Ptr< point<double> >& datPoints,
          bool useStdDevs, Ptr<double>& StdDevs,
                bool minimize, int nIterations,
          int pu, int pv,
          int *ret_nu, Ptr<double> *retU,
          int *ret_nv, Ptr<double> *retV,
          Ptr< Ptr< point<double> > > *retPw );


/*-------------------------------------------------------------------*//**
  creates a surface with the MBA algorithm. 'Dat' contains the 3d
  values and 'Dom' the locations in the parametric domain 
  (output arrays are reserved automatically)                            */
/* ---------------------------------------------------------------------*/

template< class T >
void GULAPI MBASurface( int nDat, Ptr< point<T> > Dat, Ptr< point2<T> > Dom, 
                 int nIter, int pu, int pv,          
                 int *ret_nu, Ptr<T> *retU,
                 int *ret_nv, Ptr<T> *retV,
                 Ptr< Ptr< point<T> > > *retP )
{
  T minu,maxu,minv,maxv;
  int nu,nv,i,j;
  Ptr<T> U,V,null;
  Ptr< Ptr< point1<T> > > F;
  Ptr< point<T> > coord;
  Ptr< Ptr< point<T> > > Pw;

  // reserve work arrays

  coord.reserve_pool(nDat);

  nu = pu + (1<<(nIter-1)) - 1;
  U.reserve_pool( 2*pu + (1<<(nIter-1)) + 1 );
  nv = pv + (1<<(nIter-1)) - 1;
  V.reserve_pool( 2*pv + (1<<(nIter-1)) + 1 );
  F.reserve_pool(nv+1);
  Pw.reserve_pool(nv+1);
  for( i = 0; i < nv+1; i++ )
  {
    F[i].reserve_pool(nu+1);
    Pw[i].reserve_pool(nu+1);
  }  

  // normalize domain points

  CalcBoundingBoxE( nDat, Dom, minu, maxu, minv, maxv );

  for( i = 0; i < nDat; i++ )
  {
    coord[i].x = (Dom[i].x - minu) / (maxu-minu);
    if( coord[i].x < (T)0 ) coord[i].x = (T)0;
    else if( coord[i].x > (T)1 ) coord[i].x = (T)1;

    coord[i].y = (Dom[i].y - minv) / (maxv-minv);
    if( coord[i].y < (T)0 ) coord[i].y = (T)0;
    else if( coord[i].y > (T)1 ) coord[i].y = (T)1;
  }

  /* -------- calculate x-coordinate function ----------------------------- */

  for( i = 0; i < nDat; i++ )
    coord[i].z = Dat[i].x;

  MBAapproximate<T>( nDat, coord, false, null,
                     nIter, pu, pv, U, V, F ); 

  for( i = 0; i <= nv; i++ )
    for( j = 0; j <= nu; j++ )
      Pw[i][j].x = F[i][j].x;  

  /* -------- calculate y-coordinate function ----------------------------- */

  for( i = 0; i < nDat; i++ )
    coord[i].z = Dat[i].y;

  MBAapproximate<T>( nDat, coord, false, null,
                     nIter, pu, pv, U, V, F ); 

  for( i = 0; i <= nv; i++ )
    for( j = 0; j <= nu; j++ )
      Pw[i][j].y = F[i][j].x;  

  /* -------- calculate z-coordinate function ----------------------------- */

  for( i = 0; i < nDat; i++ )
    coord[i].z = Dat[i].z;

  MBAapproximate<T>( nDat, coord, false, null,
                     nIter, pu, pv, U, V, F ); 

  for( i = 0; i <= nv; i++ )
    for( j = 0; j <= nu; j++ )
      Pw[i][j].z = F[i][j].x;  

  /* ----------- ready ---------------------------------------------------- */

  *ret_nu = nu;
  *retU = U;
  *ret_nv = nv;
  *retV = V;
  *retP = Pw;
}

// template instantiation
template GULAPI void MBASurface( 
                     int nDat, const Ptr< point<float> > Dat, const Ptr< point2<float> > Dom, 
                     int nIter, int pu, int pv,          
                     int *ret_nu, Ptr<float> *retU,
                     int *ret_nv, Ptr<float> *retV,
                     Ptr< Ptr< point<float> > > *retP );
// template instantiation
template GULAPI void MBASurface( 
                     int nDat, const Ptr< point<double> > Dat, const Ptr< point2<double> > Dom, 
                     int nIter, int pu, int pv,          
                     int *ret_nu, Ptr<double> *retU,
                     int *ret_nv, Ptr<double> *retV,
                     Ptr< Ptr< point<double> > > *retP );



}

---