shape_reconstruction: smooth_field.h Source File

Go to the documentation of this file.
00001 /****************************************************************************
00002 * VCGLib                                                            o o     *
00003 * Visual and Computer Graphics Library                            o     o   *
00004 *                                                                _   O  _   *
00005 * Copyright(C) 2014                                                \/)\/    *
00006 * Visual Computing Lab                                            /\/|      *
00007 * ISTI - Italian National Research Council                           |      *
00008 *                                                                    \      *
00009 * All rights reserved.                                                      *
00010 *                                                                           *
00011 * This program is free software; you can redistribute it and/or modify      *
00012 * it under the terms of the GNU General Public License as published by      *
00013 * the Free Software Foundation; either version 2 of the License, or         *
00014 * (at your option) any later version.                                       *
00015 *                                                                           *
00016 * This program is distributed in the hope that it will be useful,           *
00017 * but WITHOUT ANY WARRANTY; without even the implied warranty of            *
00018 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             *
00019 * GNU General Public License (http://www.gnu.org/licenses/gpl.txt)          *
00020 * for more details.                                                         *
00021 *                                                                           *
00022 ****************************************************************************/
00023 
00024 #ifndef SMOOTHER_FIELD_H
00025 #define SMOOTHER_FIELD_H
00026 
00027 //eigen stuff
00028 #include <eigenlib/Eigen/Sparse>
00029 
00030 //vcg stuff
00031 #include <vcg/complex/algorithms/update/color.h>
00032 #include <vcg/complex/algorithms/update/quality.h>
00033 #include <vcg/complex/algorithms/parametrization/tangent_field_operators.h>
00034 #include <vcg/complex/algorithms/mesh_to_matrix.h>
00035 
00036 //igl related stuff
00037 #include <igl/n_polyvector.h>
00038 #include <igl/principal_curvature.h>
00039 #include <igl/igl_inline.h>
00040 #include <igl/comiso/nrosy.h>
00041 
00042 namespace vcg {
00043 namespace tri {
00044 
00045 enum SmoothMethod{SMMiq,SMNPoly};
00046 
00047 template <class MeshType>
00048 class FieldSmoother
00049 {
00050     typedef typename MeshType::FaceType FaceType;
00051     typedef typename MeshType::VertexType VertexType;
00052     typedef typename MeshType::ScalarType ScalarType;
00053     typedef typename MeshType::CoordType CoordType;
00054 
00055 
00056     static void InitQualityByAnisotropyDir(MeshType &mesh)
00057     {
00058         std::vector<ScalarType> QVal;
00059         for (size_t i=0;i<mesh.vert.size();i++)
00060         {
00061             ScalarType N1=fabs(mesh.vert[i].K1());
00062             ScalarType N2=fabs(mesh.vert[i].K2());
00063             QVal.push_back(N1);
00064             QVal.push_back(N2);
00065         }
00066 
00067         std::sort(QVal.begin(),QVal.end());
00068         int percUp=int(floor(QVal.size()*0.95+0.5));
00069         ScalarType trimUp=QVal[percUp];
00070 
00071         for (size_t i=0;i<mesh.vert.size();i++)
00072         {
00073             ScalarType N1=(mesh.vert[i].K1());
00074             ScalarType N2=(mesh.vert[i].K2());
00075 
00076            ScalarType NMax=std::max(N1,N2)/trimUp;
00077            ScalarType NMin=std::min(N1,N2)/trimUp;
00078 
00079            if (NMax>1)NMax=1;
00080            if (NMax<-1)NMax=-1;
00081 
00082            if (NMin>1)NMin=1;
00083            if (NMin<-1)NMin=-1;
00084 
00085            ScalarType CurvAni=(NMax-NMin)/2;
00086            mesh.vert[i].Q()=CurvAni;
00087         }
00088         vcg::tri::UpdateQuality<MeshType>::FaceFromVertex(mesh);
00089     }
00090 
00091     static void SetEdgeDirection(FaceType *f,int edge)
00092     {
00093         CoordType dir=f->P0(edge)-f->P1(edge);
00094         dir.Normalize();
00095         ScalarType prod1=fabs(dir*f->PD1());
00096         ScalarType prod2=fabs(dir*f->PD2());
00097         if (prod1>prod2)
00098         {
00099             f->PD1()=dir;
00100             f->PD2()=f->N()^dir;
00101         }else
00102         {
00103             f->PD2()=dir;
00104             f->PD1()=f->N()^dir;
00105         }
00106     }
00107 
00108     static void AddSharpEdgesConstraints(MeshType & mesh,
00109                                          const ScalarType &thr=0.2)
00110     {
00111         for (size_t i=0;i<mesh.face.size();i++)
00112             for (int j=0;j<mesh.face[i].VN();j++)
00113             {
00114                 FaceType *f0=&mesh.face[i];
00115                 FaceType *f1=f0->FFp(j);
00116                 if (f0==f1)continue;
00117                 CoordType N0=f0->N();
00118                 CoordType N1=f1->N();
00119                 if ((N0*N1)>thr)continue;
00120                 SetEdgeDirection(f0,j);
00121                 f0->SetS();
00122             }
00123     }
00124 
00125     static void AddBorderConstraints(MeshType & mesh)
00126     {
00127         for (size_t i=0;i<mesh.face.size();i++)
00128             for (int j=0;j<mesh.face[i].VN();j++)
00129             {
00130                 FaceType *f0=&mesh.face[i];
00131                 FaceType *f1=f0->FFp(j);
00132                 assert(f1!=NULL);
00133                 if (f0!=f1)continue;
00134                 SetEdgeDirection(f0,j);
00135                 f0->SetS();
00136             }
00137     }
00138 
00139     static void AddCurvatureConstraints(MeshType & mesh,const ScalarType &thr=0.5)
00140     {
00141         for (size_t i=0;i<mesh.face.size();i++)
00142             if (mesh.face[i].Q()>thr)mesh.face[i].SetS();
00143     }
00144 
00145     //hard constraints have selected face
00146     static void CollectHardConstraints( MeshType & mesh,
00147                                     Eigen::VectorXi &HardI,
00148                                     Eigen::MatrixXd &HardD,
00149                                     SmoothMethod SMethod,
00150                                     int Ndir)
00151     {
00152         //count number of hard constraints
00153         int numS=vcg::tri::UpdateSelection<MeshType>::FaceCount(mesh);
00154         HardI=Eigen::MatrixXi(numS,1);
00155         if ((Ndir==2)||(SMethod==SMMiq))
00156             HardD=Eigen::MatrixXd(numS,3);
00157         else
00158             HardD=Eigen::MatrixXd(numS,6);
00159         //then update them
00160         int curr_index=0;
00161         for (size_t i=0;i<mesh.face.size();i++)
00162         {
00163             if (!mesh.face[i].IsS())continue;
00164 
00165             CoordType dir=mesh.face[i].PD1();
00166             dir.Normalize();
00167 
00168             HardI(curr_index,0)=i;
00169 
00170             HardD(curr_index,0)=dir.X();
00171             HardD(curr_index,1)=dir.Y();
00172             HardD(curr_index,2)=dir.Z();
00173             if ((Ndir==4)&&(SMethod==SMNPoly))
00174             {
00175                 dir=mesh.face[i].PD2();
00176                 HardD(curr_index,3)=dir.X();
00177                 HardD(curr_index,4)=dir.Y();
00178                 HardD(curr_index,5)=dir.Z();
00179             }
00180             curr_index++;
00181 
00182         }
00183 
00184     }
00185 
00186     //hard constraints have selected face
00187     static void CollectSoftConstraints( MeshType & mesh,
00188                                         Eigen::VectorXi &SoftI,
00189                                         Eigen::MatrixXd &SoftD,
00190                                         Eigen::VectorXd &SoftW)
00191     {
00192         //count number of soft constraints
00193         int numS=vcg::tri::UpdateSelection<MeshType>::FaceCount(mesh);
00194         numS=mesh.fn-numS;
00195         //allocate eigen matrix
00196         SoftI=Eigen::MatrixXi(numS,1);
00197         SoftD=Eigen::MatrixXd(numS,3);
00198         SoftW=Eigen::MatrixXd(numS,1);
00199 
00200         //then update them
00201         int curr_index=0;
00202         for (size_t i=0;i<mesh.face.size();i++)
00203         {
00204             if (mesh.face[i].IsS())continue;
00205 
00206             CoordType dir=mesh.face[i].PD1();
00207             dir.Normalize();
00208 
00209             SoftI(curr_index,0)=i;
00210 
00211             SoftD(curr_index,0)=dir.X();
00212             SoftD(curr_index,1)=dir.Y();
00213             SoftD(curr_index,2)=dir.Z();
00214 
00215             SoftW(curr_index,0)=mesh.face[i].Q();
00216 
00217             curr_index++;
00218 
00219         }
00220     }
00221 
00222     static void SmoothMIQ(MeshType &mesh,
00223                           Eigen::VectorXi &HardI,   //hard constraints index
00224                           Eigen::MatrixXd &HardD,   //hard directions
00225                           Eigen::VectorXi &SoftI,   //soft constraints
00226                           Eigen::MatrixXd &SoftD,   //weight of soft constraints
00227                           Eigen::VectorXd &SoftW,   //soft directions
00228                           ScalarType alpha_soft,
00229                           int Ndir)
00230     {
00231 
00232         assert((Ndir==2)||(Ndir==4));
00233         Eigen::MatrixXi F;
00234         Eigen::MatrixXd V;
00235 
00236         MeshToMatrix<MeshType>::GetTriMeshData(mesh,F,V);
00237 
00238         Eigen::MatrixXd output_field;
00239         Eigen::VectorXd output_sing;
00240 
00241         igl::nrosy(V,F,HardI,HardD,SoftI,SoftW,SoftD,Ndir,alpha_soft,output_field,output_sing);
00242 
00243         //finally update the principal directions
00244         for (size_t i=0;i<mesh.face.size();i++)
00245         {
00246             CoordType dir1;
00247             dir1[0]=output_field(i,0);
00248             dir1[1]=output_field(i,1);
00249             dir1[2]=output_field(i,2);
00250 
00251             dir1.Normalize();
00252             CoordType dir2=mesh.face[i].N()^dir1;
00253             dir2.Normalize();
00254 
00255             ScalarType Norm1=mesh.face[i].PD1().Norm();
00256             ScalarType Norm2=mesh.face[i].PD2().Norm();
00257 
00258             mesh.face[i].PD1()=dir1*Norm1;
00259             mesh.face[i].PD2()=dir2*Norm2;
00260         }
00261     }
00262 
00263     static void SmoothNPoly(MeshType &mesh,
00264                             Eigen::VectorXi &HardI,   //hard constraints index
00265                             Eigen::MatrixXd &HardD,   //hard directions
00266                             int Ndir)
00267     {
00268         assert((Ndir==2)||(Ndir==4));
00269 
00270         Eigen::MatrixXi F;
00271         Eigen::MatrixXd V;
00272 
00273         MeshToMatrix<MeshType>::GetTriMeshData(mesh,F,V);
00274 
00275         Eigen::MatrixXd output_field;
00276         Eigen::VectorXd output_sing;
00277 
00278         igl::n_polyvector(V,F,HardI,HardD,output_field);
00279 
00280         //finally update the principal directions
00281         for (size_t i=0;i<mesh.face.size();i++)
00282         {
00283             CoordType dir1;
00284             dir1[0]=output_field(i,0);
00285             dir1[1]=output_field(i,1);
00286             dir1[2]=output_field(i,2);
00287 
00288             dir1.Normalize();
00289             CoordType dir2=mesh.face[i].N()^dir1;
00290             dir2.Normalize();
00291 
00292             ScalarType Norm1=mesh.face[i].PD1().Norm();
00293             ScalarType Norm2=mesh.face[i].PD2().Norm();
00294 
00295             mesh.face[i].PD1()=dir1*Norm1;
00296             mesh.face[i].PD2()=dir2*Norm2;
00297         }
00298     }
00299 
00300     static void PickRandomDir(MeshType &mesh,
00301                               int &indexF,
00302                               CoordType &Dir)
00303     {
00304         indexF=rand()%mesh.fn;
00305         FaceType *currF=&mesh.face[indexF];
00306         CoordType dirN=currF->N();
00307         dirN.Normalize();
00308         Dir=CoordType(1,0,0);
00309         if (fabs(Dir*dirN)>0.9)
00310             Dir=CoordType(0,1,0);
00311         if (fabs(Dir*dirN)>0.9)
00312             Dir=CoordType(0,0,1);
00313 
00314         Dir=dirN^Dir;
00315         Dir.Normalize();
00316     }
00317 
00318 public:
00319 
00320     struct SmoothParam
00321     {
00322         //the 90° rotation independence while smoothing the direction field
00323         int Ndir;
00324         //the weight of curvature if doing the smoothing keeping the field close to the original one
00325         ScalarType alpha_curv;
00326         //align the field to border or not
00327         bool align_borders;
00328         //threshold to consider some edge as sharp feature and to use as hard constraint (0, not use)
00329         ScalarType sharp_thr;
00330         //threshold to consider some edge as high curvature anisotropyand to use as hard constraint (0, not use)
00331         ScalarType curv_thr;
00332         //the method used to smooth MIQ or "Designing N-PolyVector Fields with Complex Polynomials"
00333         SmoothMethod SmoothM;
00334         //the number of faces of the ring used ot esteem the curvature
00335         int curvRing;
00336         //this are additional hard constraints
00337         std::vector<std::pair<int,CoordType> > AddConstr;
00338 
00339         SmoothParam()
00340         {
00341             Ndir=4;
00342             curvRing=2;
00343             alpha_curv=0.0;
00344 
00345             align_borders=false;
00346 
00347             SmoothM=SMMiq;
00348 
00349             sharp_thr=0.0;
00350             curv_thr=0.4;
00351         }
00352 
00353     };
00354 
00355     static void SelectConstraints(MeshType &mesh,SmoothParam &SParam)
00356     {
00357         //clear all selected faces
00358         vcg::tri::UpdateFlags<MeshType>::FaceClear(mesh);
00359 
00360         //add curvature hard constraints
00361         //ScalarType Ratio=mesh.bbox.Diag()*0.01;
00362 
00363         if (SParam.curv_thr>0)
00364             AddCurvatureConstraints(mesh,SParam.curv_thr);
00365 
00366          //add alignment to sharp features
00367         if (SParam.sharp_thr>0)
00368             AddSharpEdgesConstraints(mesh,SParam.sharp_thr);
00369 
00370         //add border constraints
00371         if (SParam.align_borders)
00372             AddBorderConstraints(mesh);
00373 
00374         //aff final constraints
00375         for (int i=0;i<SParam.AddConstr.size();i++)
00376         {
00377             int indexF=SParam.AddConstr[i].first;
00378             CoordType dir=SParam.AddConstr[i].second;
00379             mesh.face[indexF].PD1()=dir;
00380             mesh.face[indexF].PD2()=mesh.face[indexF].N()^dir;
00381             mesh.face[indexF].PD1().Normalize();
00382             mesh.face[indexF].PD2().Normalize();
00383             mesh.face[indexF].SetS();
00384         }
00385     }
00386 
00387     static void GloballyOrient(MeshType &mesh)
00388     {
00389         vcg::tri::CrossField<MeshType>::MakeDirectionFaceCoherent(mesh,true);
00390     }
00391 
00392 
00393     static void InitByCurvature(MeshType & mesh,
00394                                 int Nring,
00395                                 bool UpdateFaces=true)
00396     {
00397 
00398         tri::RequirePerVertexCurvatureDir(mesh);
00399 
00400         Eigen::MatrixXi F;
00401         Eigen::MatrixXd V;
00402 
00403         Eigen::MatrixXd PD1,PD2,PV1,PV2;
00404         MeshToMatrix<MeshType>::GetTriMeshData(mesh,F,V);
00405         igl::principal_curvature(V,F,PD1,PD2,PV1,PV2,Nring);
00406 
00407         //then copy curvature per vertex
00408         for (size_t i=0;i<mesh.vert.size();i++)
00409         {
00410             mesh.vert[i].PD1()=CoordType(PD1(i,0),PD1(i,1),PD1(i,2));
00411             mesh.vert[i].PD2()=CoordType(PD2(i,0),PD2(i,1),PD2(i,2));
00412             mesh.vert[i].K1()=PV1(i,0);
00413             mesh.vert[i].K2()=PV2(i,0);
00414         }
00415         if (!UpdateFaces)return;
00416         vcg::tri::CrossField<MeshType>::SetFaceCrossVectorFromVert(mesh);
00417         InitQualityByAnisotropyDir(mesh);
00418     }
00419 
00420     static void SmoothDirections(MeshType &mesh,
00421                                  int Ndir,
00422                                  SmoothMethod SMethod=SMNPoly,
00423                                  bool HardAsS=true,
00424                                  ScalarType alphaSoft=0)
00425     {
00426 
00427         Eigen::VectorXi HardI;   //hard constraints
00428         Eigen::MatrixXd HardD;   //hard directions
00429         Eigen::VectorXi SoftI;   //soft constraints
00430         Eigen::VectorXd SoftW;   //weight of soft constraints
00431         Eigen::MatrixXd SoftD;   //soft directions
00432 
00433         if (HardAsS)
00434             CollectHardConstraints(mesh,HardI,HardD,SMethod,Ndir);
00435 
00436         //collect soft constraints , miw only one that allows for soft constraints
00437         if ((alphaSoft>0)&&(SMethod==SMMiq))
00438           CollectSoftConstraints(mesh,SoftI,SoftD,SoftW);
00439 
00440         //add some hard constraints if are not present
00441         int numC=3;
00442         if ((SoftI.rows()==0)&&(HardI.rows()==0))
00443         {
00444             printf("Add Forced Constraints \n");
00445             fflush(stdout);
00446             HardI=Eigen::MatrixXi(numC,1);
00447 
00448             if ((Ndir==4)&&(SMethod==SMNPoly))
00449                 HardD=Eigen::MatrixXd(numC,6);
00450             else
00451                 HardD=Eigen::MatrixXd(numC,3);
00452 
00453             for (int i=0;i<numC;i++)
00454             {
00455                 CoordType Dir;
00456                 int indexF;
00457                 PickRandomDir(mesh,indexF,Dir);
00458 
00459                 HardI(i,0)=indexF;
00460                 HardD(i,0)=Dir.X();
00461                 HardD(i,1)=Dir.Y();
00462                 HardD(i,2)=Dir.Z();
00463 
00464                 if ((Ndir==4)&&(SMethod==SMNPoly))
00465                 {
00466                     CoordType Dir1=mesh.face[indexF].N()^Dir;
00467                     Dir1.Normalize();
00468                     HardD(i,3)=Dir1.X();
00469                     HardD(i,4)=Dir1.Y();
00470                     HardD(i,5)=Dir1.Z();
00471                 }
00472             }
00473         }
00474 
00475          //finally smooth
00476         if (SMethod==SMMiq)
00477          SmoothMIQ(mesh,HardI,HardD,SoftI,SoftD,SoftW,alphaSoft,Ndir);
00478         else
00479         {
00480             assert(SMethod==SMNPoly);
00481             SmoothNPoly(mesh,HardI,HardD,Ndir);
00482         }
00483     }
00484 
00485 
00486     static void SmoothDirections(MeshType &mesh,SmoothParam SParam)
00487     {
00488         //for the moment only cross and line field
00489 
00490         //initialize direction by curvature if needed
00491         if ((SParam.alpha_curv>0)||
00492              (SParam.sharp_thr>0)||
00493              (SParam.curv_thr>0))
00494             InitByCurvature(mesh,SParam.curvRing);
00495 
00496         SelectConstraints(mesh,SParam);
00497         //then do the actual smooth
00498         SmoothDirections(mesh,SParam.Ndir,SParam.SmoothM,true,SParam.alpha_curv);
00499     }
00500 
00501 };
00502 
00503 } // end namespace tri
00504 } // end namespace vcg
00505 #endif // SMOOTHER_FIELD_H