From 0a750b1bc56e81b787f68224ddbe1bfdea282bc3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Tr=C3=B6ger?= Date: Sun, 22 Dec 2013 12:11:01 +0100 Subject: [PATCH] only use scale value for rotation-only solving, not the transformed clusters. This gives the more expected results --- .../App/opendcm/core/imp/kernel_imp.hpp | 27 +++++------ .../App/opendcm/module3d/imp/solver_imp.hpp | 45 ++++++++++++++----- .../Assembly/App/opendcm/module3d/solver.hpp | 3 +- 3 files changed, 45 insertions(+), 30 deletions(-) diff --git a/src/Mod/Assembly/App/opendcm/core/imp/kernel_imp.hpp b/src/Mod/Assembly/App/opendcm/core/imp/kernel_imp.hpp index 9e9f8085d..cd2ae8d5a 100644 --- a/src/Mod/Assembly/App/opendcm/core/imp/kernel_imp.hpp +++ b/src/Mod/Assembly/App/opendcm/core/imp/kernel_imp.hpp @@ -21,6 +21,7 @@ #define DCM_KERNEL_IMP_H #include "../kernel.hpp" +#include #include #include @@ -62,9 +63,7 @@ void Dogleg::calculateStep(const Eigen::MatrixBase& g, const Ei const typename Kernel::Vector h_gn = jacobi.fullPivLu().solve(-residual); const double eigen_error = (jacobi*h_gn + residual).norm(); #ifdef USE_LOGGING - BOOST_LOG(log)<< "Gauss Newton error: "<::calculateStep(const Eigen::MatrixBase& g, const Ei // compute the dogleg step if(h_gn.norm() <= delta) { h_dl = h_gn; -#ifdef USE_LOGGING - BOOST_LOG(log)<< "Gauss Newton Step: "<= delta) { //h_dl = alpha*h_sd; @@ -94,8 +90,6 @@ void Dogleg::calculateStep(const Eigen::MatrixBase& g, const Ei BOOST_LOG(log)<< "Unnormal dogleg descent detected: "<::calculateStep(const Eigen::MatrixBase& g, const Ei BOOST_LOG(log)<< "Unnormal dogleg beta detected: "<::MappedEquationSystem::MappedEquationSystem(int params m_parameter(params), m_residual(equations), m_params(params), m_eqns(equations), Scaling(1.), Jacobi(&m_jacobi(0,0),equations,params,DynStride(equations,1)), - Parameter(&m_parameter(0),params,DynStride(1,1)), + Parameter(&m_parameter(0),params,DynStride(1,1)), Residual(&m_residual(0),equations,DynStride(1,1)) { m_param_rot_offset = 0; @@ -421,21 +413,22 @@ void Kernel::MappedEquationSystem::setAccess(AccessType t) { if(t==complete) { new(&Jacobi) MatrixMap(&m_jacobi(0,0),m_eqns,m_params,DynStride(m_eqns,1)); new(&Parameter) VectorMap(&m_parameter(0),m_params,DynStride(1,1)); - new(&Residual) VectorMap(&m_residual(0), m_eqns, DynStride(1,1)); + new(&Residual) VectorMap(&m_residual(0), m_eqns, DynStride(1,1)); } else if(t==rotation) { int num = m_param_trans_offset; new(&Jacobi) MatrixMap(&m_jacobi(0,0),m_eqn_trans_offset,num,DynStride(m_eqns,1)); new(&Parameter) VectorMap(&m_parameter(0),num,DynStride(1,1)); - new(&Residual) VectorMap(&m_residual(0), m_eqn_trans_offset, DynStride(1,1)); + new(&Residual) VectorMap(&m_residual(0), m_eqn_trans_offset, DynStride(1,1)); } else if(t==general) { int num = m_params - m_param_trans_offset; int eq_num = m_eqns - m_eqn_trans_offset; new(&Jacobi) MatrixMap(&m_jacobi(m_eqn_trans_offset,m_param_trans_offset),eq_num,num,DynStride(m_eqns,1)); new(&Parameter) VectorMap(&m_parameter(m_param_trans_offset),num,DynStride(1,1)); - new(&Residual) VectorMap(&m_residual(m_eqn_trans_offset), eq_num, DynStride(1,1)); + new(&Residual) VectorMap(&m_residual(m_eqn_trans_offset), eq_num, DynStride(1,1)); } + m_access = t; }; @@ -443,11 +436,11 @@ template class Nonlinear> bool Kernel::MappedEquationSystem::hasAccessType(AccessType t) { if(t==rotation) - return (m_param_rot_offset>0); + return (m_param_rot_offset>0 && m_eqn_rot_offset>0); else if(t==general) - return (m_param_trans_offset0); + return (m_params>0 && m_eqns>0); }; template class Nonlinear> diff --git a/src/Mod/Assembly/App/opendcm/module3d/imp/solver_imp.hpp b/src/Mod/Assembly/App/opendcm/module3d/imp/solver_imp.hpp index b4b843c43..77e6e0fc7 100644 --- a/src/Mod/Assembly/App/opendcm/module3d/imp/solver_imp.hpp +++ b/src/Mod/Assembly/App/opendcm/module3d/imp/solver_imp.hpp @@ -107,12 +107,12 @@ SystemSolver::Rescaler::Rescaler(boost::shared_ptr c, Mes& m) : cl template void SystemSolver::Rescaler::operator()() { - mes.Scaling = scaleClusters(); + mes.Scaling = scaleClusters(calculateScale()); rescales++; }; template -typename SystemSolver::Scalar SystemSolver::Rescaler::scaleClusters() { +typename SystemSolver::Scalar SystemSolver::Rescaler::calculateScale() { typedef typename Cluster::cluster_iterator citer; std::pair cit = cluster->clusters(); @@ -134,6 +134,12 @@ typename SystemSolver::Scalar SystemSolver::Rescaler::scaleClusters() sc = (s>sc) ? s : sc; } + return sc; +} + +template +typename SystemSolver::Scalar SystemSolver::Rescaler::scaleClusters(Scalar sc) { + //if no scaling-value returned we can use 1 sc = (Kernel::isSame(sc,0, 1e-10)) ? 1. : sc; @@ -357,20 +363,33 @@ void SystemSolver::solveCluster(boost::shared_ptr cluster, Sys& sy //if(!has_cycle) { #ifdef USE_LOGGING - BOOST_LOG(log)<< "non-cyclic system dedected: solve rotation only"; + BOOST_LOG(log)<< "non-cyclic system dedected: solve rotation only"; #endif - //cool, lets do uncylic. first all rotational constraints with rotational parameters - mes.setAccess(rotation); + //cool, lets do uncylic. first all rotational constraints with rotational parameters + mes.setAccess(rotation); - //solve can be done without catching exceptions, because this only fails if the system is - //unsolvable. We need the rescaler here two as unscaled rotations may be unprecise if the - //parts are big (as ne the normals are always 1) and then distances may not be possible - //to solve correctly with only translations - Rescaler re(cluster, mes); - re(); + //rotations need to be calculated in a scaled manner. thats because the normales used for + //rotation calculation are always 1, no matter how big the part is. This can lead to problems + //when for example two rotated faces have a precision error on the parallel normals but a distance + //at the outer edges is far bigger than the precision as the distance from normal origin to outer edge + //is bigger 1. that would lead to unsolvable translation-only systems. + + //solve need to catch exceptions to reset the mes scaling on failure + Rescaler re(cluster, mes); + mes.Scaling = 1./(re.calculateScale()*SKALEFAKTOR); + + try { sys.kernel().solve(mes, re); + mes.Scaling = 1.; + } + catch(...) { + mes.Scaling = 1.; + throw; + } + + //now let's see if we have to go on with the translations + if(mes.hasAccessType(general)) { - //now let's see if we have to go on with the translations mes.setAccess(general); mes.recalculate(); @@ -392,6 +411,8 @@ void SystemSolver::solveCluster(boost::shared_ptr cluster, Sys& sy done = false; } } + }; + //}; //not done already? try it the hard way! diff --git a/src/Mod/Assembly/App/opendcm/module3d/solver.hpp b/src/Mod/Assembly/App/opendcm/module3d/solver.hpp index 0379b7d14..45e750ac9 100644 --- a/src/Mod/Assembly/App/opendcm/module3d/solver.hpp +++ b/src/Mod/Assembly/App/opendcm/module3d/solver.hpp @@ -84,7 +84,8 @@ struct SystemSolver : public Job { void operator()(); - Scalar scaleClusters(); + Scalar calculateScale(); + Scalar scaleClusters(Scalar sc); void collectPseudoPoints(boost::shared_ptr parent, LocalVertex cluster, std::vector