OSGQuadTreeCreator.cpp

00001 /*---------------------------------------------------------------------------*\
00002  *                           OpenSG NURBS Library                            *
00003  *                                                                           *
00004  *                                                                           *
00005  * Copyright (C) 2001-2006 by the University of Bonn, Computer Graphics Group*
00006  *                                                                           *
00007  *                         http://cg.cs.uni-bonn.de/                         *
00008  *                                                                           *
00009  * contact: edhellon@cs.uni-bonn.de, guthe@cs.uni-bonn.de, rk@cs.uni-bonn.de *
00010  *                                                                           *
00011 \*---------------------------------------------------------------------------*/
00012 /*---------------------------------------------------------------------------*\
00013  *                                License                                    *
00014  *                                                                           *
00015  * This library is free software; you can redistribute it and/or modify it   *
00016  * under the terms of the GNU Library General Public License as published    *
00017  * by the Free Software Foundation, version 2.                               *
00018  *                                                                           *
00019  * This library is distributed in the hope that it will be useful, but       *
00020  * WITHOUT ANY WARRANTY; without even the implied warranty of                *
00021  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU         *
00022  * Library General Public License for more details.                          *
00023  *                                                                           *
00024  * You should have received a copy of the GNU Library General Public         *
00025  * License along with this library; if not, write to the Free Software       *
00026  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.                 *
00027  *                                                                           *
00028 \*---------------------------------------------------------------------------*/
00029 /*---------------------------------------------------------------------------*\
00030  *                                Changes                                    *
00031  *                                                                           *
00032  *                                                                           *
00033  *                                                                           *
00034  *                                                                           *
00035  *                                                                           *
00036  *                                                                           *
00037 \*---------------------------------------------------------------------------*/
00038 #include "OSGQuadTreeCreator.h"
00039 #include <stdio.h>
00040
00041 OSG_USING_NAMESPACE
00042
00043 #ifdef WIN32
00044 #pragma warning (disable : 985)
00045 #endif
00046 
00047 #ifdef _DEBUG
00048  #ifdef WIN32
00049   #undef THIS_FILE
00050 static char THIS_FILE[] = __FILE__;
00051  #endif
00052 #endif
00053 
00054 double QuadTreeCreator::
00055 computeApproximationError(DCTPFace *f)
00056 {
00057 //the approximate error comes from to error,
00058 //I: interpolation with a bilinear surface
00059 //II: interpolating this bilinear surface with two triangles
00060
00061 //first gotta get needed controlpoints from the index-ed quadtreeleaf
00062     BezierTensorSurface *btsp =
00063         static_cast<BezierTensorSurface*>(f->faceinfo);
00064
00065     const DCTPVec4dmatrix &          cps = btsp->getControlPointMatrix();
00066     const DCTPVec4dmatrix::size_type m   = cps.size() - 1;
00067     //size of control point matrix in x direction
00068     const DCTPVec4dvector::size_type n = cps[0].size() - 1;
00069     //size of cpm. in y direction
00070     Vec3d b00;
00071     Vec3d b0n;
00072     Vec3d bm0;
00073     Vec3d bmn;
00074     b00[0] = cps[0][0][0] / cps[0][0][3];
00075     b00[1] = cps[0][0][1] / cps[0][0][3];
00076     b00[2] = cps[0][0][2] / cps[0][0][3];
00077     b0n[0] = (cps[0][n][0] / cps[0][n][3]);
00078     b0n[1] = (cps[0][n][1] / cps[0][n][3]);
00079     b0n[2] = (cps[0][n][2] / cps[0][n][3]);
00080     bm0[0] = (cps[m][0][0] / cps[m][0][3]);
00081     bm0[1] = (cps[m][0][1] / cps[m][0][3]);
00082     bm0[2] = (cps[m][0][2] / cps[m][0][3]);
00083     bmn[0] = (cps[m][n][0] / cps[m][n][3]);
00084     bmn[1] = (cps[m][n][1] / cps[m][n][3]);
00085     bmn[2] = (cps[m][n][2] / cps[m][n][3]);
00086     //now computin' I. error:
00087     //we hafta decide a maxmimum value:
00088     double                     max_quad = -1.0;
00089     Vec3d                      cij, bij;
00090     double                     quad_size;
00091     DCTPVec3dmatrix::size_type i;
00092     DCTPVec3dvector::size_type j;
00093     const double               rn = 1.0 / n;
00094     const double               rm = 1.0 / m;
00095
00096     Vec3d       ci0 = b00;
00097     Vec3d       cin = b0n;
00098     const Vec3d dc0 = (bm0 - b00) * rm;
00099     const Vec3d dcn = (bmn - b0n) * rm;
00100 //  Vec2d cl_uv;
00101 //  int i_err = 0;
00102
00103 //  if( m_bForTrimming )
00104     {
00105 //      cl_uv[0] = 0.0;
00106         for(i = 0; i <= m; ++i)
00107         {
00108 //          cl_uv[1] = 0.0;
00109             cij = ci0;
00110             const Vec3d dci = (cin - ci0) * rn;
00111
00112             for(j = 0; j <= n; ++j)
00113             {
00114                 if(cps[i][j][3] < DCTP_EPS)
00115                 {
00116                     // the weight is zero so we return a practically
00117                     // infinite error enforcing a subdivision of the patch
00118                     quad_size = 1.0e300;
00119                 }
00120                 else
00121                 {
00122                     bij[0] = (cps[i][j][0] / cps[i][j][3]) - cij[0];
00123                     bij[1] = (cps[i][j][1] / cps[i][j][3]) - cij[1];
00124                     bij[2] = (cps[i][j][2] / cps[i][j][3]) - cij[2];
00125 //              bij = btsp->computewdeCasteljau( cl_uv, i_err ) - cij;
00126                     quad_size = bij.squareLength();
00127                 }
00128                 if(quad_size > max_quad)
00129                     max_quad = quad_size;
00130                 cij += dci;
00131 //              cl_uv[1] += rn;
00132             } //end for
00133
00134             ci0 += dc0;
00135             cin += dcn;
00136 //          cl_uv[0] += rm;
00137         }
00138
00139         const double error_I = sqrt(max_quad);
00140
00141         //now get error2
00142         const Vec3d  linerr_vect = (b00 - b0n + bmn - bm0);
00143         const double error_II    = sqrt(linerr_vect.squareLength() ) * 0.25;
00144         //now the error iz the sum o' I and II:
00145         return error_I + error_II;
00146     }
00147 /*  else
00148     {
00149         Vec3d       cl_norm;
00150         const Vec3d cdc0 = ( b0n - b00 ) * rn;
00151         const Vec3d cdcm = ( bmn - bm0 ) * rn;
00152         Vec3d       c0j;
00153         Vec3d       cmj;
00154 
00155         for( i = 0; i <= m; ++i )
00156         {
00157             cij = ci0;
00158             c0j = b00;
00159             cmj = bm0;
00160 
00161             const Vec3d dci = ( cin - ci0 ) * rn;
00162 
00163 //          std::cerr << dci << std::endl;
00164 
00165             for( j = 0; j <= n; ++j )
00166             {
00167                 const Vec3d dcj = cmj - c0j;
00168 
00169                 bij = cps[ i ][ j ] - cij;
00170                 cl_norm = dci.cross( dcj );
00171 //              std::cerr << dci << " x " << dcj << std::endl;
00172 //              std::cerr << cl_norm << "," << bij << std::endl;
00173                 quad_size = cl_norm.squareLength( );
00174                 if( quad_size > 0.0 )
00175                 {
00176                     cl_norm *= 1.0 / sqrt( quad_size );
00177                     quad_size = osgabs( bij.dot( cl_norm ) );
00178 //                  quad_size = bij.squareLength();
00179                     if ( quad_size > max_quad ) max_quad = quad_size;
00180                 }
00181                 cij += dci;
00182                 c0j += cdc0;
00183                 cmj += cdcm;
00184             } //end for
00185             ci0 += dc0;
00186             cin += dcn;
00187         }
00188 //      std::cerr << max_quad << std::endl;
00189         return max_quad;
00190     }*/
00191 }
00192
00193
00194 int QuadTreeCreator::setInitialLeaves(
00195     const beziersurfacematrix& patches,
00196     const DCTPdvector&         intervals_u,
00197     const DCTPdvector&         intervals_v)
00198 {
00199 //sets up the basic leaves array of the quadtree
00200 //FIXME: no size & sanity checks made to input, they may be wrong
00201     beziersurfacematrix::size_type u_size = patches.size();
00202     beziersurfacevector::size_type v_size = patches[0].size();
00203     Vec3d                          a,b,c,d;
00204
00205     for(beziersurfacematrix::size_type u = 0; u < u_size; ++u)
00206     {
00207         for(beziersurfacevector::size_type v = 0; v < v_size; ++v)
00208         {
00209             a[0] = intervals_u[u]; a[1] = intervals_v[v + 1]; a[2] = 0.0;
00210             b[0] = intervals_u[u + 1]; b[1] = intervals_v[v + 1]; b[2] = 0.0;
00211             c[0] = intervals_u[u + 1]; c[1] = intervals_v[v]; c[2] = 0.0;
00212             d[0] = intervals_u[u]; d[1] = intervals_v[v]; d[2] = 0.0;
00213
00214 //          std::cerr << a << b << c << d << std::endl;
00215
00216             DCTPFace *face = qtm->AddQuad(a, b, c, d, 0.0);
00217             face->faceinfo = new BezierTensorSurface;
00218 #ifdef OSG_ADAPTIVE_QUAD_TREE
00219             face->norm = -1.0;
00220 #endif
00221             *(static_cast<BezierTensorSurface*>(face->faceinfo)) = patches[u][v];
00222             BezierTensorSurface temp_surface = patches[u][v];
00223         } //end inner for
00224
00225     }
00226
00227     return 0;
00228 }
00229
00230 int QuadTreeCreator::finishSubdivisions(DCTPFace *f)
00231 {
00232 /*        BezierTensorSurface *btsp = (BezierTensorSurface*)f->faceinfo;
00233         beziersurfacematrix created_surfaces( 0 );
00234         if ( btsp->midPointSubDivision( created_surfaces ) ) {
00235                 std::cerr << "QuadTreeCreator::finisSubdivisions:" <<
00236                         "Cannot midpointsubdivide BezierSurf'!" << std::endl;
00237                 return -1;
00238         }
00239         dctpfacevector::size_type idx = qtm->faces.size() - 3;
00240         (qtm->faces[ idx ])->faceinfo = new BezierTensorSurface;
00241         *((BezierTensorSurface*)(qtm->faces[ idx++ ])->faceinfo) =
00242                 created_surfaces[ 1 ][ 1 ];
00243         (qtm->faces[ idx ])->faceinfo = new BezierTensorSurface;
00244         *((BezierTensorSurface*)(qtm->faces[ idx++ ])->faceinfo) =
00245                 created_surfaces[ 1 ][ 0 ];
00246         (qtm->faces[ idx ])->faceinfo = new BezierTensorSurface;
00247         *((BezierTensorSurface*)(qtm->faces[ idx++ ])->faceinfo) =
00248                 created_surfaces[ 0 ][ 0 ];
00249 
00250         delete btsp;
00251         f->faceinfo = new BezierTensorSurface;
00252         *((BezierTensorSurface*)f->faceinfo) = created_surfaces[ 0 ][ 1 ];
00253         return 0;*/
00254
00255     BezierTensorSurface *btsp = static_cast<BezierTensorSurface*>(f->faceinfo);
00256     beziersurfacevector  created_surfaces;
00257     if(btsp->midPointSubDivision(created_surfaces) )
00258     {
00259         std::cerr << "QuadTreeCreator::finisSubdivisions:" <<
00260         "Cannot midpointsubdivide BezierSurf'!" << std::endl;
00261         return -1;
00262     }
00263     dctpfacevector::size_type idx = qtm->faces.size() - 3;
00264
00265     (qtm->faces[idx])->faceinfo                                         = new BezierTensorSurface;
00266     *(static_cast<BezierTensorSurface*>((qtm->faces[idx++])->faceinfo)) =
00267         created_surfaces[2];
00268
00269 /*      std::cerr << std::endl;
00270         std::cerr << ( qtm->faces[ idx ] )->orig_face[ 0 ]->coords << std::endl;
00271         std::cerr << ( qtm->faces[ idx ] )->orig_face[ 1 ]->coords << std::endl;
00272         std::cerr << ( qtm->faces[ idx ] )->orig_face[ 2 ]->coords << std::endl;
00273         std::cerr << ( qtm->faces[ idx ] )->orig_face[ 3 ]->coords << std::endl;
00274         std::cerr << created_surfaces[ 1 ].getControlPointMatrix( )[ 0 ][ 0 ] << std::endl;*/
00275
00276     (qtm->faces[idx])->faceinfo                                         = new BezierTensorSurface;
00277     *(static_cast<BezierTensorSurface*>((qtm->faces[idx++])->faceinfo)) =
00278         created_surfaces[1];
00279
00280     (qtm->faces[idx])->faceinfo = btsp;
00281
00282     f->faceinfo                                       = new BezierTensorSurface;
00283     *(static_cast<BezierTensorSurface*>(f->faceinfo)) = created_surfaces[0];
00284     return 0;
00285 }
00286
00287
00288 int QuadTreeCreator::createQuadTree(void)
00289 {
00290 //now this is the essentials of this whole class
00291 //all other are 'eye-candys'
00292
00293     for(dctpfacevector::size_type f = 0;
00294         f < qtm->faces.size(); ++f)
00295     {
00296         DCTPFace *face = qtm->faces[f];
00297         while(computeApproximationError(face) > error_epsilon)
00298         {
00299             qtm->SubdivideQuad(face);
00300             if(finishSubdivisions(face) )
00301                 return -1;
00302         }
00303         face->norm = error_epsilon / computeBilinearNorm(face);
00304     }
00305
00306     return 0;     //all OK, whew...
00307 }
00308
00309 double QuadTreeCreator::computeBilinearNorm(DCTPFace *face)
00310 {
00311     BezierTensorSurface *btsp =
00312         static_cast<BezierTensorSurface*>(face->faceinfo);
00313     DCTPVec4dmatrix&           cps = btsp->getControlPointMatrix();
00314     DCTPVec4dmatrix::size_type m   = cps.size() - 1;
00315     DCTPVec4dvector::size_type n   = cps[0].size() - 1;
00316     Vec3d                      b00;
00317     Vec3d                      b0n;
00318     Vec3d                      bm0;
00319     Vec3d                      bmn;
00320     b00[0] = cps[0][0][0] / cps[0][0][3];
00321     b00[1] = cps[0][0][1] / cps[0][0][3];
00322     b00[2] = cps[0][0][2] / cps[0][0][3];
00323     b0n[0] = (cps[0][n][0] / cps[0][n][3]) - b00[0];
00324     b0n[1] = (cps[0][n][1] / cps[0][n][3]) - b00[1];
00325     b0n[2] = (cps[0][n][2] / cps[0][n][3]) - b00[2];
00326     bm0[0] = (cps[m][0][0] / cps[m][0][3]) - b00[0];
00327     bm0[1] = (cps[m][0][1] / cps[m][0][3]) - b00[1];
00328     bm0[2] = (cps[m][0][2] / cps[m][0][3]) - b00[2];
00329     bmn[0] = (cps[m][n][0] / cps[m][n][3]) - b00[0];
00330     bmn[1] = (cps[m][n][1] / cps[m][n][3]) - b00[1];
00331     bmn[2] = (cps[m][n][2] / cps[m][n][3]) - b00[2];
00332     double tmp1,tmp2;
00333     tmp1 = sqrt(b0n.squareLength() );
00334     tmp2 = sqrt(bm0.squareLength() );
00335     if(tmp1 < tmp2)
00336         tmp1 = tmp2;
00337     tmp2 = sqrt(bmn.squareLength() ) / sqrt(2.f);        //1.414213; //divided by length
00338     if(tmp1 < tmp2)
00339         return tmp2;
00340     else
00341         return tmp1;
00342 }
00343
00344 void QuadTreeCreator::resetMesh(void)
00345 {
00346     dctpfacevector::iterator fe = qtm->faces.end();
00347
00348     for(dctpfacevector::iterator f = qtm->faces.begin(); f != fe; ++f)
00349         delete static_cast<BezierTensorSurface*>((*f)->faceinfo);
00350
00351     qtm->reinit();
00352 }