---

gunu_project.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 "gul_types.h"
#include "gul_float.h"
#include "gul_vector.h"
#include "guma_newton.h"
#include "gunu_basics.h"
#include "gunu_refine.h"
#include "gunu_sample.h"
#include "gunu_bezier_derivatives.h"
#include "guma_random.h"
#include "guma_rkdtree.h"
#include "gunu_project.h"

#include <iostream>
#include "gul_io.h"

namespace gunu {

using gul::rtr;
using gul::Ptr;
using gul::point;
using gul::hpoint;
using gul::kdpoint;
using gul::List;
using gul::ListNode;
using gul::rel_equal;
using guma::rkdtree;
using guma::kdnnrec;
using guma::kdrec;
using guma::newton_info;
using guma::random_des;
using guma::Newton;
using gul::distance;

/*---------------------------------------------------------------------
  help classes, used by Newton iteration code
---------------------------------------------------------------------*/

template< class T, class HP >
struct patch_ninfo : public gul::pool_object
{
  int pu;
  int pv;
  Ptr<T> U;
  Ptr<T> V;
  Ptr< Ptr<HP> > Pw;

  int row;  /* position of bezier segment in nurbs surface */
  int col;    
};

template< class T, class HP, class EP >
class point_ninfo : public kdpoint<T,EP>
{
public:
  Ptr< patch_ninfo<T,HP> > S;

  T u;
  T v;

  point_ninfo *next;
  point_ninfo *prev;

  point_ninfo( EP& inP, Ptr< patch_ninfo<T,HP> >& inS, T in_u, T in_v )
    : kdpoint<T,EP>( inP, 0 )
  {
    S = inS;
    u = in_u;
    v = in_v;
  }

};

template< class T, class HP, class EP >
class bezier_ninfo : public newton_info<T>
{
public:
  T                           last_u;
  T                           last_v;
  bool                        last_result;
  const EP                   *P;
  const point_ninfo<T,HP,EP> *SP;
  Ptr< Ptr<EP> >              SKL;
  T                           tol1;
  T                           tol2;
  bool                        init;

  bezier_ninfo( T in_tol1, T in_tol2 )
          : newton_info<T>( 2, (T)0.0001, in_tol2 )
  {
    int i;
    
    tol1 = in_tol1;
    tol2 = in_tol2;
  
    init = false;
    P = 0;
    SP = 0;
   
    SKL.reserve_pool(3);
    for( i = 0; i < 3; i++ )
      SKL[i].reserve_pool(3);
  }
  
  virtual bool evaluate( const gul::Ptr<T>& dom );

  void setCurrentPoint( const EP *inP, const point_ninfo<T,HP,EP> *inSP )
  {
    init = false;

    P = inP;
    SP = inSP; 
  }
};


template< class T, class HP, class EP >
bool bezier_ninfo<T,HP,EP>::evaluate( const Ptr<T>& dom )
{
  T u = dom[0], v = dom[1], rootPS;
  EP S, Su, Sv, Suu, Svv, Suv, PS;
  T SuSu, SvSv, b1, b2;

  if( init && (u == last_u) && (v == last_v) )
    return last_result;
  else
  {
    init = true;
    last_u = u;
    last_v = v;
  }   

  if( u < 0.0 ) u = 0.0;
  else if( u > 1.0 ) u = 1.0;  

  if( v < 0.0 ) v = 0.0;
  else if( v > 1.0 ) v = 1.0;  

  BezierSurfaceDerivatives( 
              u, v,  SP->S[0].pu, SP->S[0].U, SP->S[0].pv, SP->S[0].V, 
              SP->S[0].Pw, 2, SKL );

  S = SKL[0][0];
  Su = SKL[1][0];
  Sv = SKL[0][1];
  Suu = SKL[2][0];
  Suv = SKL[1][1];
  Svv = SKL[0][2];

  PS = S - *P;
  SuSu = Su*Su;
  SvSv = Sv*Sv;

  fjac[0][0] = PS*Suu + SuSu;
  fjac[0][1] = fjac[1][0] = PS*Suv + Su*Sv;
  fjac[1][1] = PS*Svv + SvSv;
  
  fvec[0] = Su*PS;
  fvec[1] = Sv*PS;

  if( rel_equal(S,*P,tol1) )
  {
    last_result = true;
    return true;
  }

  rootPS = rtr<T>::sqrt(PS*PS);
  b1 = rtr<T>::sqrt(SuSu)*rootPS;
  b2 = rtr<T>::sqrt(SvSv)*rootPS;

  if( rel_equal(fvec[0]+b1,b1,tol2) &&
      rel_equal(fvec[1]+b2,b2,tol2) )
  {
    last_result = true;
    return true;
  }

  last_result = false;
  return false;
}

/*---------------------------------------------------------------------
  Projects a set of points onto a nurbs surface
---------------------------------------------------------------------*/

template< class T, class HP, class EP >
GULAPI int ProjectToSurface(
           int nDatPoints, const Ptr<EP>& datPoints,
           int ref_nu, int ref_pu, const Ptr<T>& refU,
           int ref_nv, int ref_pv, const Ptr<T>& refV,
           const Ptr< Ptr<HP> >& refPw, 
           Ptr< pts_point<T,EP> >& retPts )
{
  Ptr<T> BezU,BezV;
  Ptr< Ptr< Ptr<HP> > > strips;
  Ptr< Ptr< Ptr<HP> > > segs;
  Ptr<T> sampU, sampV;
  Ptr< Ptr<EP> > sampP;
  int nStrips, nSegs, i, j, k, l;
  int nSampU, nSampV;
  int nKu, nKv;
  Ptr<int> Ku,Mu,Kv,Mv;
  bezier_ninfo<T,HP,EP> I((T)0.000001,(T)0.000001);
  Ptr< patch_ninfo<T,HP> > S;
  point_ninfo<T,HP,EP> *P;
  List< point_ninfo<T,HP,EP> > LP;
  rkdtree< kdpoint<T,EP>,T,random_des > KDT(3);
  Ptr< kdnnrec< kdpoint<T,EP>,T > > N; 
  ListNode< kdrec< kdpoint<T,EP>,T > > *r;
  Ptr<T> dom,x1,x2;  
  point_ninfo<T,HP,EP> *pi;
  int n;
  bool found;
  int index,nFound;
  T u,v,nurbs_u,nurbs_v,min_u,min_v,dist,min_dist;
  EP min_S;
  int nRetPts;

  nRetPts = 0;

  dom.reserve_pool(2);
  x1.reserve_pool(2);
  x2.reserve_pool(2);

  /* form Bezier knot vectors */

  BezU.reserve_pool(ref_pu+ref_pu+2);
  BezV.reserve_pool(ref_pv+ref_pv+2);
  
  BezierKnotVector(ref_pu,BezU);
  BezierKnotVector(ref_pv,BezV);

  /* ----- divide surface into Bezier segments ----------------------- */

  nSampU = ref_pu*3;
  nSampV = ref_pv*3;

  sampP.reserve_pool(nSampV);
  for( i = 0; i < nSampV; i++ )
    sampP[i].reserve_pool(nSampU);
  
  sampU.reserve_pool(nSampU);
  sampV.reserve_pool(nSampV);
  EqualSpacedParameters( nSampU, BezU[0], BezU[ref_pu+1], sampU );
  EqualSpacedParameters( nSampV, BezV[0], BezV[ref_pv+1], sampV );

  nKu = BezierPositions( ref_nu, ref_pu, refU, Ku, Mu );
  nKv = BezierPositions( ref_nv, ref_pv, refV, Kv, Mv );
  nStrips = nKv;
  nSegs = nKu;

  // form horizontal Bezier strips:

  BezierDecomposeSurfaceV( nKv, Kv, Mv, ref_nu, ref_pu, refU, 
                           ref_nv, ref_pv, refV, refPw, &strips );
                           
  // divide horizontal strips into Bezier segments:

  segs.reserve_pool( nStrips );

  for( i = 0; i < nStrips; i++ )
  {
    BezierDecomposeSurfaceU( nKu, Ku, Mu, ref_nu, ref_pu, refU, 
                             ref_pv, ref_pv, BezV, strips[i], &segs );

    strips[i].reset(); // delete memory of strips, not needed anymore

    for( j = 0; j < nSegs; j++ )
    {
      CalcSurfaceMesh( ref_pu, ref_pu, BezU, ref_pv, ref_pv, BezV, segs[j],
                       nSampU, sampU, nSampV, sampV, sampP );

      S.reserve_pool(1);

      S[0].pu = ref_pu;
      S[0].pv = ref_pv;
      S[0].U = BezU;
      S[0].V = BezV;
      S[0].Pw = segs[j];
      S[0].row = i;
      S[0].col = j;
      
      for( k = 0; k < nSampV; k++ )
      {
        for( l = 0; l < nSampU; l++ )
        {
          P = new point_ninfo<T,HP,EP>( sampP[k][l], S, sampU[l], sampV[k] );
          LP.Append( P ); 

          KDT.insert( P );
        }
      }
    }
    // control points are referenced now by the patch infos, and the 'strips'
    // array is still needed for the next strip, so it is neither necessary 
    // nor possible to free something ! 
  }

  // search for each data point a nearest neighbor, and then do the newton
  // iteration

  x1[0] = x1[1] = (T)0;
  x2[0] = x2[1] = (T)1;
  
  for( i = 0; i < nDatPoints; i++ )
  {
    N.reserve_pool(1); // this ensures that memory from previous runs is freed
                     // (though this isn't necessary since N only has 1 element)

    n = KDT.nearest_neighbors( kdpoint<T,EP>(datPoints[i],0), 1, N );
    if( !n ) continue;

    /*
    cout << "data point #" << i << ": " << datPoints[i] << "\n";
    cout << "nearest neighbors with rkdtree\n";
    for( j = 0; j < n; j++ )
      N[j].dump(3);
    */

    nFound = 0;
    for( r = N[0].m_recs.First(); r; r = r->next )
    {
      pi = (point_ninfo<T,HP,EP> *)r->el.m_base;
      I.setCurrentPoint( &datPoints[i], pi ); 

      dom[0] = pi->u;
      dom[1] = pi->v;
      found = Newton( dom, x1, x2, I, 200 );
      if( !found ) continue;

      // calculate position in nurbs surface domain

      u = I.last_u;
      v = I.last_v;

      index = pi->S[0].row-1;
      if( index < 0 )
        nurbs_v = ((T)1-v)*refV[ref_pu] + v*refV[Kv[index+1]];
      else
        nurbs_v = ((T)1-v)*refV[Kv[index]] + v*refV[Kv[index+1]];    

      index = pi->S[0].col-1;
      if( index < 0 )
        nurbs_u = ((T)1-u)*refU[ref_pu] + u*refU[Ku[index+1]];
      else
        nurbs_u = ((T)1-u)*refU[Ku[index]] + u*refU[Ku[index+1]];    

      // if multiple points were found, take the nearest

      dist = distance( I.SKL[0][0], datPoints[i] );

      if( !nFound || (dist < min_dist) )
      {
        min_dist = dist;
        min_u = nurbs_u;
        min_v = nurbs_v;
        min_S = I.SKL[0][0];
      }
      nFound++;
      
      /*
      cout << "point calculated by newton iteration: " << I.SKL[0][0] << "\n";
      */
    }

    // append a record to the result list
    
    if( nFound )
    {
      retPts[nRetPts].i = i;
      retPts[nRetPts].u = min_u;
      retPts[nRetPts].v = min_v;
      retPts[nRetPts].S = min_S;
      nRetPts++;
    }
  }

  return nRetPts;
}
// template instantiation
template GULAPI int ProjectToSurface(
           int nDatPoints, const Ptr< point<float> >& datPoints,
           int ref_nu, int ref_pu, const Ptr<float>& refU,
           int ref_nv, int ref_pv, const Ptr<float>& refV,
           const Ptr< Ptr< hpoint<float> > >& Pw,
           Ptr< pts_point< float,point<float> > >& retPts );
template GULAPI int ProjectToSurface(
           int nDatPoints, const Ptr< point<double> >& datPoints,
           int ref_nu, int ref_pu, const Ptr<double>& refU,
           int ref_nv, int ref_pv, const Ptr<double>& refV,
           const Ptr< Ptr< hpoint<double> > >& Pw,
           Ptr< pts_point< double,point<double> > >& retPts );
template GULAPI int ProjectToSurface(
           int nDatPoints, const Ptr< point<float> >& datPoints,
           int ref_nu, int ref_pu, const Ptr<float>& refU,
           int ref_nv, int ref_pv, const Ptr<float>& refV,
           const Ptr< Ptr< point<float> > >& Pw,
           Ptr< pts_point< float,point<float> > >& retPts );
template GULAPI int ProjectToSurface(
           int nDatPoints, const Ptr< point<double> >& datPoints,
           int ref_nu, int ref_pu, const Ptr<double>& refU,
           int ref_nv, int ref_pv, const Ptr<double>& refV,
           const Ptr< Ptr< point<double> > >& Pw,
           Ptr< pts_point< double,point<double> > >& retPts );

}

---