src/mesh/cfMesh/utilities/helperClasses/geometry/quadricFitting/quadricFittingI.H

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | cfMesh: A library for mesh generation
   4    \\    /   O peration     |
   5     \\  /    A nd           | Author: Franjo Juretic (franjo.juretic@c-fields.com)
   6      \\/     M anipulation  | Copyright (C) Creative Fields, Ltd.
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of cfMesh.
  10
  11     cfMesh is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     cfMesh is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with cfMesh.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25
  26 \*---------------------------------------------------------------------------*/
  27
  28 #include "quadricFitting.H"
  29
  30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  31
  32 namespace Foam
  33 {
  34
  35 // * * * * * * * * * Private Member Functions * * * * * * * * * * * * * * * * //
  36
  37 void quadricFitting::calculateNormalVector()
  38 {
  39     symmTensor mat(symmTensor::zero);
  40
  41     forAll(otherPoints_, i)
  42     {
  43         const vector d = (otherPoints_[i] - origin_);
  44
  45         mat += symm(d * d);
  46     }
  47
  48     mat /= otherPoints_.size();
  49
  50     //- find the eigenvalues of the tensor
  51     const vector ev = eigenValues(mat);
  52
  53     //- estimate the normal as the eigenvector associated
  54     //- to the smallest eigenvalue
  55     normal_ = eigenVector(mat, ev[0]);
  56 }
  57
  58 void quadricFitting::calculateCoordinateSystem()
  59 {
  60     if( mag(normal_) < VSMALL )
  61     {
  62         vecX_ = vecY_ = vector::zero;
  63         forAll(transformedPoints_, i)
  64             transformedPoints_[i] = vector::zero;
  65
  66         return;
  67     }
  68
  69     const plane pl(origin_, normal_);
  70
  71     bool valid(false);
  72
  73     label pI(0);
  74     while( !valid )
  75     {
  76         const point vp = pl.nearestPoint(otherPoints_[pI]);
  77         vecX_ = vp - origin_;
  78
  79         if( magSqr(vecX_) < VSMALL )
  80         {
  81             ++pI;
  82             continue;
  83         }
  84
  85         vecX_ /= mag(vecX_);
  86         vecY_ = normal_ ^ vecX_;
  87         valid = true;
  88     }
  89
  90     if( !valid )
  91     {
  92         FatalErrorIn
  93         (
  94             "void quadricFitting::calculateCoordinateSystem()"
  95         ) << "Cannot find a valid coordinate system!" << endl;
  96
  97         normal_ = vecX_ = vecY_ = vector::zero;
  98     }
  99
 100     //- transform other points to this new coordinate system
 101     transformedPoints_.setSize(otherPoints_.size());
 102
 103     forAll(otherPoints_, i)
 104     {
 105         const vector delta = (otherPoints_[i] - origin_);
 106         point& p = transformedPoints_[i];
 107
 108         p.x() = delta & vecX_;
 109         p.y() = delta & vecY_;
 110         p.z() = delta & normal_;
 111     }
 112 }
 113
 114 void quadricFitting::calculateQuadricCoeffs()
 115 {
 116     if( mag(normal_) < VSMALL )
 117     {
 118         coefficients_.setSize(5);
 119         coefficients_ = 0.0;
 120         return;
 121     }
 122
 123     simpleMatrix<scalar> mat(5, 0.0, 0.0);
 124     DynList<scalar> helper;
 125     helper.setSize(5);
 126
 127     forAll(transformedPoints_, i)
 128     {
 129         const point& p = transformedPoints_[i];
 130
 131         helper[0] = sqr(p.x());
 132         helper[1] = sqr(p.y());
 133         helper[2] = p.x() * p.y();
 134         helper[3] = p.x();
 135         helper[4] = p.y();
 136
 137         for(label rowI=0;rowI<5;++rowI)
 138         {
 139             for(label colI=rowI;colI<5;++colI)
 140                 mat[rowI][colI] += helper[rowI] * helper[colI];
 141
 142             mat.source()[rowI] += helper[rowI] * p.z();
 143         }
 144     }
 145
 146     for(label rowI=1;rowI<5;++rowI)
 147     {
 148         for(label colI=0;colI<rowI;++colI)
 149             mat[rowI][colI] = mat[colI][rowI];
 150
 151         if( mag(mat[rowI][rowI]) < SMALL )
 152             mat[rowI][rowI] = SMALL;
 153     }
 154
 155     coefficients_ = mat.LUsolve();
 156 }
 157
 158 void quadricFitting::calculateBestFit()
 159 {
 160     label iteration(0);
 161
 162     while( iteration++ < 10 )
 163     {
 164         calculateCoordinateSystem();
 165
 166         calculateQuadricCoeffs();
 167
 168         if( (mag(coefficients_[3]) > SMALL) || (mag(coefficients_[4]) > SMALL) )
 169         {
 170             //- correct the normal
 171             const scalar d =
 172                 Foam::sqrt(1 + sqr(coefficients_[3]) + sqr(coefficients_[4]));
 173
 174             const vector newNormal
 175             (
 176                 normal_ / d -
 177                 coefficients_[3] * vecX_ / d -
 178                 coefficients_[4] * vecY_ / d
 179             );
 180
 181             normal_ = newNormal;
 182         }
 183         else
 184         {
 185             break;
 186         }
 187     }
 188
 189     # ifdef DEBUGQuadric
 190     Info << "Other points " << otherPoints_ << endl;
 191     Info << "Transformed points " << transformedPoints_ << endl;
 192     Info << "normal " << normal_ << endl;
 193     Info << "Coefficients " << coefficients_ << endl;
 194     # endif
 195 }
 196
 197 // * * * * * * * * * * * * * * * *  Constructors * * * * * * * * * * * * * * //
 198
 199 template<class ListType>
 200 inline quadricFitting::quadricFitting
 201 (
 202     const point& origin,
 203     const vector normal,
 204     const ListType& otherPoints
 205 )
 206 :
 207     origin_(origin),
 208     normal_(normal),
 209     vecX_(vector::zero),
 210     vecY_(vector::zero),
 211     otherPoints_(),
 212     transformedPoints_(),
 213     coefficients_()
 214 {
 215     otherPoints_.setSize(otherPoints.size());
 216     forAll(otherPoints, i)
 217         otherPoints_[i] = otherPoints[i];
 218
 219     if( magSqr(normal_) < VSMALL )
 220     {
 221         FatalErrorIn
 222         (
 223             "template<class ListType>\n"
 224             "inline quadricFitting::quadricFitting"
 225             "(const point& origin, const vector normal,"
 226             " const ListType& otherPoints)"
 227         ) << "Cannot construct the quadric surface for point "
 228             << origin_ << " because the normal does not exist!"
 229             << "\nThis indicates that the input"
 230             << " surface mesh is of poor quality" << exit(FatalError);
 231
 232         normal_ = vector::zero;
 233     }
 234     else
 235     {
 236         normal_ /= mag(normal_);
 237     }
 238
 239     calculateBestFit();
 240 }
 241
 242 //- Construct from point and other points
 243 template<class ListType>
 244 inline quadricFitting::quadricFitting
 245 (
 246     const point& origin,
 247     const ListType& otherPoints
 248 )
 249 :
 250     origin_(origin),
 251     normal_(),
 252     vecX_(),
 253     vecY_(),
 254     otherPoints_(),
 255     transformedPoints_(),
 256     coefficients_()
 257 {
 258     otherPoints_.setSize(otherPoints.size());
 259     forAll(otherPoints, i)
 260         otherPoints_[i] = otherPoints[i];
 261
 262     calculateNormalVector();
 263
 264     calculateBestFit();
 265 }
 266
 267 inline quadricFitting::~quadricFitting()
 268 {
 269 }
 270
 271 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 272
 273 inline const vector& quadricFitting::normal() const
 274 {
 275     return normal_;
 276 }
 277
 278 //- Return min curvature
 279 inline scalar quadricFitting::minCurvature() const
 280 {
 281     return
 282     (
 283         coefficients_[0] + coefficients_[1] -
 284         sqrt(sqr(coefficients_[0] - coefficients_[1]) + sqr(coefficients_[2]))
 285     );
 286 }
 287
 288 //- Return max curvature
 289 inline scalar quadricFitting::maxCurvature() const
 290 {
 291     return
 292     (
 293         coefficients_[0] + coefficients_[1] +
 294         sqrt(sqr(coefficients_[0] - coefficients_[1]) + sqr(coefficients_[2]))
 295     );
 296 }
 297
 298 //- Return mean curvature
 299 inline scalar quadricFitting::meanCurvature() const
 300 {
 301     return 0.5 * (minCurvature() + maxCurvature());
 302 }
 303
 304 //- Return Gaussian curvature
 305 inline scalar quadricFitting::gaussianCurvature() const
 306 {
 307     return minCurvature() * maxCurvature();
 308 }
 309
 310 //- Return min curvature vector
 311 inline vector quadricFitting::minCurvatureVector() const
 312 {
 313     const scalar theta =
 314     0.5 * Foam::atan
 315     (
 316         coefficients_[2] /
 317         stabilise(coefficients_[0] - coefficients_[1], VSMALL)
 318     );
 319
 320     return Foam::cos(theta) * vecX_ + Foam::sin(theta) * vecY_;
 321 }
 322
 323 //- Return max curvature vector
 324 inline vector quadricFitting::maxCurvatureVector() const
 325 {
 326     const scalar theta =
 327         0.5 * Foam::atan
 328         (
 329             coefficients_[2] /
 330             stabilise(coefficients_[0] - coefficients_[1], VSMALL)
 331         );
 332
 333     return Foam::sin(theta) * vecX_ + Foam::cos(theta) * vecY_;
 334 }
 335
 336 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 337
 338 } // End namespace Foam
 339
 340 // ************************************************************************* //