Revision: 20622
          
http://projects.blender.org/plugins/scmsvn/viewcvs.php?view=rev&root=bf-blender&revision=20622
Author:   ben2610
Date:     2009-06-04 14:47:59 +0200 (Thu, 04 Jun 2009)

Log Message:
-----------
iTaSC: better auto substep calculation with formula to evaluate the amount of 
singularity.

Modified Paths:
--------------
    branches/ge_dev/intern/itasc/CopyPose.cpp
    branches/ge_dev/intern/itasc/CopyPose.hpp
    branches/ge_dev/intern/itasc/Scene.cpp
    branches/ge_dev/intern/itasc/Scene.hpp
    branches/ge_dev/intern/itasc/Solver.hpp
    branches/ge_dev/intern/itasc/WDLSSolver.cpp
    branches/ge_dev/intern/itasc/WDLSSolver.hpp
    branches/ge_dev/intern/itasc/WSDLSSolver.cpp
    branches/ge_dev/intern/itasc/WSDLSSolver.hpp
    branches/ge_dev/source/blender/ikplugin/intern/itasc_plugin.cpp

Modified: branches/ge_dev/intern/itasc/CopyPose.cpp
===================================================================
--- branches/ge_dev/intern/itasc/CopyPose.cpp   2009-06-04 11:16:56 UTC (rev 
20621)
+++ branches/ge_dev/intern/itasc/CopyPose.cpp   2009-06-04 12:47:59 UTC (rev 
20622)
@@ -457,4 +457,16 @@
        return m_values; 
 }
 
+double CopyPose::getMaxTimestep(double& timestep)
+{
+       // CopyPose should not have any limit on linear velocity: 
+       // in case the target is out of reach, this can be very high.
+       // We will simply limit on rotation
+       e_scalar maxChidot = m_chidot.block(3,0,3,1).cwise().abs().maxCoeff();
+       if (timestep*maxChidot > m_maxDeltaChi) {
+               timestep = m_maxDeltaChi/maxChidot;
+       }
+       return timestep;
 }
+
+}

Modified: branches/ge_dev/intern/itasc/CopyPose.hpp
===================================================================
--- branches/ge_dev/intern/itasc/CopyPose.hpp   2009-06-04 11:16:56 UTC (rev 
20621)
+++ branches/ge_dev/intern/itasc/CopyPose.hpp   2009-06-04 12:47:59 UTC (rev 
20622)
@@ -22,6 +22,7 @@
     virtual void initCache(Cache *_cache);
     virtual void updateControlOutput(const Timestamp& timestamp);
        virtual void CopyPose::modelUpdate(Frame& _external_pose,const 
Timestamp& timestamp);
+       virtual double getMaxTimestep(double& timestep);
 
 public:
     enum ID {          // constraint ID in callback and setControlParameter

Modified: branches/ge_dev/intern/itasc/Scene.cpp
===================================================================
--- branches/ge_dev/intern/itasc/Scene.cpp      2009-06-04 11:16:56 UTC (rev 
20621)
+++ branches/ge_dev/intern/itasc/Scene.cpp      2009-06-04 12:47:59 UTC (rev 
20622)
@@ -29,7 +29,7 @@
        m_ncTotal(0),m_nqTotal(0),m_nuTotal(0),m_nsets(0),
        m_solver(NULL),m_cache(NULL) 
 {
-
+       m_minstep = 0.01;
 }
 
 Scene::~Scene() 
@@ -230,6 +230,7 @@
                numsubstep = 1;
        double timesubstep = timestep/numsubstep;
        double timeleft = timestep;
+       e_scalar nlcoef;
        // initially we keep timestep unchanged so that update function compute 
the velocity over
        while (numsubstep > 0) {
                // get objects
@@ -347,7 +348,7 @@
                }
 
                //Call the solver with A, Wq, Wy, ydot to solver qdot:
-               if(!m_solver->solve(m_A,m_Wy,m_ydot,m_Wq,m_qdot))
+               if(!m_solver->solve(m_A,m_Wy,m_ydot,m_Wq,m_qdot,nlcoef))
                        // this should never happen
                        return false;
                //send result to the objects
@@ -382,6 +383,11 @@
                        // and joint limit gain variation
                        // We will pass the joint velocity to each object and 
they will recommend a maximum timestep
                        timesubstep = timeleft;
+                       double maxsubstep = timestep/nlcoef;
+                       if (maxsubstep < m_minstep)
+                               maxsubstep = m_minstep;
+                       if (timesubstep > maxsubstep)
+                               timesubstep = maxsubstep;
                        for(ObjectMap::iterator 
it=objects.begin();it!=objects.end();++it){
                                Object_struct* os = it->second;
                                if(os->object->getType()==Object::Controlled)
@@ -391,7 +397,7 @@
                                ConstraintSet_struct* cs = it->second;
                                cs->task->getMaxTimestep(timesubstep);
                        }
-                       if (timesubstep >= timeleft-KDL::epsilon) {
+                       if (timesubstep >= timeleft-(m_minstep/2.0)) {
                                timesubstep = timeleft;
                                numsubstep = 1;
                                timeleft = 0.;

Modified: branches/ge_dev/intern/itasc/Scene.hpp
===================================================================
--- branches/ge_dev/intern/itasc/Scene.hpp      2009-06-04 11:16:56 UTC (rev 
20621)
+++ branches/ge_dev/intern/itasc/Scene.hpp      2009-06-04 12:47:59 UTC (rev 
20622)
@@ -39,7 +39,7 @@
        e_matrix6 m_Vf,m_Uf;
     e_vector m_Wy,m_ydot,m_qdot,m_xdot;
        e_vector6 m_Sf,m_tempf;
-
+       double m_minstep;
     unsigned int m_ncTotal,m_nqTotal,m_nuTotal,m_nsets;
        std::vector<bool> m_ytask;
 

Modified: branches/ge_dev/intern/itasc/Solver.hpp
===================================================================
--- branches/ge_dev/intern/itasc/Solver.hpp     2009-06-04 11:16:56 UTC (rev 
20621)
+++ branches/ge_dev/intern/itasc/Solver.hpp     2009-06-04 12:47:59 UTC (rev 
20622)
@@ -20,7 +20,7 @@
        // gc = grouping of constraint output , 
        //      size of vector = nc, alternance of true / false to indicate the 
grouping of output
        virtual bool init(unsigned int nq, unsigned int nc, const 
std::vector<bool>& gc)=0;
-    virtual bool solve(const e_matrix& A, const e_vector& Wy, const e_vector& 
ydot, const e_matrix& Wq, e_vector& qdot)=0;
+    virtual bool solve(const e_matrix& A, const e_vector& Wy, const e_vector& 
ydot, const e_matrix& Wq, e_vector& qdot, e_scalar& nlcoef)=0;
 };
 
 }

Modified: branches/ge_dev/intern/itasc/WDLSSolver.cpp
===================================================================
--- branches/ge_dev/intern/itasc/WDLSSolver.cpp 2009-06-04 11:16:56 UTC (rev 
20621)
+++ branches/ge_dev/intern/itasc/WDLSSolver.cpp 2009-06-04 12:47:59 UTC (rev 
20622)
@@ -32,7 +32,7 @@
     return true;
 }
 
-bool WDLSSolver::solve(const e_matrix& A, const e_vector& Wy, const e_vector& 
ydot, const e_matrix& Wq, e_vector& qdot){
+bool WDLSSolver::solve(const e_matrix& A, const e_vector& Wy, const e_vector& 
ydot, const e_matrix& Wq, e_vector& qdot, e_scalar& nlcoef){
 // Create the Weighted jacobian
     m_AWq = A*Wq;
        for (int i=0; i<Wy.size(); i++)
@@ -52,10 +52,26 @@
     //U'*Wy'*ydot
        m_WyUt_ydot = (m_WyU.transpose()*ydot).lazy();
     //S^-1*U'*Wy'*ydot
-    for(int i=0;i<m_WyUt_ydot.size();++i)
-        m_SinvWyUt_ydot(i) = 
m_WyUt_ydot(i)*m_S(i)/(m_S(i)*m_S(i)+m_lambda*m_lambda);
+       e_scalar maxDeltaS = e_scalar(0.0);
+       e_scalar prevS = e_scalar(0.0);
+       e_scalar maxS = e_scalar(1.0);
+       for(int i=0;i<m_WyUt_ydot.size();++i) {
+               e_scalar S = m_S(i);
+               if (i > 0 && S > KDL::epsilon) {
+                       if ((prevS-S) > maxDeltaS) {
+                               maxDeltaS = (prevS-S);
+                               maxS = prevS;
+                       }
+               }
+        m_SinvWyUt_ydot(i) = m_WyUt_ydot(i)*S/(S*S+m_lambda*m_lambda);
+               prevS = S;
+       }
     //qdot=Wq*V*S^-1*U'*Wy'*ydot
     qdot=(m_WqV*m_SinvWyUt_ydot).lazy();
+       if (maxDeltaS == e_scalar(0.0))
+               nlcoef = e_scalar(1.0/KDL::epsilon);
+       else
+               nlcoef = maxS/(maxS-maxDeltaS)/e_scalar(2.0);
     return true;
 }
 

Modified: branches/ge_dev/intern/itasc/WDLSSolver.hpp
===================================================================
--- branches/ge_dev/intern/itasc/WDLSSolver.hpp 2009-06-04 11:16:56 UTC (rev 
20621)
+++ branches/ge_dev/intern/itasc/WDLSSolver.hpp 2009-06-04 12:47:59 UTC (rev 
20622)
@@ -22,7 +22,7 @@
     virtual ~WDLSSolver();
 
     virtual bool init(unsigned int nq, unsigned int nc, const 
std::vector<bool>& gc);
-    virtual bool solve(const e_matrix& A, const e_vector& Wy, const e_vector& 
ydot, const e_matrix& Wq, e_vector& qdot);
+    virtual bool solve(const e_matrix& A, const e_vector& Wy, const e_vector& 
ydot, const e_matrix& Wq, e_vector& qdot, e_scalar& nlcoef);
 
     void setLambda(double _lambda){m_lambda=_lambda;};
 

Modified: branches/ge_dev/intern/itasc/WSDLSSolver.cpp
===================================================================
--- branches/ge_dev/intern/itasc/WSDLSSolver.cpp        2009-06-04 11:16:56 UTC 
(rev 20621)
+++ branches/ge_dev/intern/itasc/WSDLSSolver.cpp        2009-06-04 12:47:59 UTC 
(rev 20622)
@@ -41,7 +41,7 @@
     return true;
 }
 
-bool WSDLSSolver::solve(const e_matrix& A, const e_vector& Wy, const e_vector& 
ydot, const e_matrix& Wq, e_vector& qdot)
+bool WSDLSSolver::solve(const e_matrix& A, const e_vector& Wy, const e_vector& 
ydot, const e_matrix& Wq, e_vector& qdot, e_scalar& nlcoef)
 {
        unsigned int i, j, l;
        e_scalar N, M;
@@ -59,12 +59,22 @@
        m_Wy_ydot = Wy.cwise() * ydot;
     m_WqV = (Wq*m_V).lazy();
        qdot.fill(e_scalar(0.));
+       e_scalar maxDeltaS = e_scalar(0.0);
+       e_scalar prevS = e_scalar(0.0);
+       e_scalar maxS = e_scalar(1.0);
        for(i=0;i<m_ns;++i) {
                e_scalar norm, mag, alpha, _qmax, Sinv, vmax, damp;
+               e_scalar S = m_S(i);
                bool prev;
-               if (m_S(i) < KDL::epsilon)
+               if (S < KDL::epsilon)
                        break;
-               Sinv = e_scalar(1.)/m_S(i);
+               Sinv = e_scalar(1.)/S;
+               if (i > 0) {
+                       if ((prevS-S) > maxDeltaS) {
+                               maxDeltaS = (prevS-S);
+                               maxS = prevS;
+                       }
+               }
                N = M = e_scalar(0.);
                for (l=0, prev=m_ytask[0], norm=e_scalar(0.); l<m_nc; l++) {
                        if (prev == m_ytask[l]) {
@@ -100,7 +110,12 @@
                        damp = Sinv*alpha;
                }
                qdot += m_WqV.col(i)*damp;
+               prevS = S;
        }
+       if (maxDeltaS == e_scalar(0.0))
+               nlcoef = e_scalar(1.0/KDL::epsilon);
+       else
+               nlcoef = maxS/(maxS-maxDeltaS)/e_scalar(2.0);
     return true;
 }
 

Modified: branches/ge_dev/intern/itasc/WSDLSSolver.hpp
===================================================================
--- branches/ge_dev/intern/itasc/WSDLSSolver.hpp        2009-06-04 11:16:56 UTC 
(rev 20621)
+++ branches/ge_dev/intern/itasc/WSDLSSolver.hpp        2009-06-04 12:47:59 UTC 
(rev 20622)
@@ -24,7 +24,7 @@
     virtual ~WSDLSSolver();
 
     virtual bool init(unsigned int _nq, unsigned int _nc, const 
std::vector<bool>& gc);
-    virtual bool solve(const e_matrix& A, const e_vector& Wy, const e_vector& 
ydot, const e_matrix& Wq, e_vector& qdot);
+    virtual bool solve(const e_matrix& A, const e_vector& Wy, const e_vector& 
ydot, const e_matrix& Wq, e_vector& qdot, e_scalar& nlcoef);
 
     void setQmax(double _qmax){m_qmax=_qmax;};
 

Modified: branches/ge_dev/source/blender/ikplugin/intern/itasc_plugin.cpp
===================================================================
--- branches/ge_dev/source/blender/ikplugin/intern/itasc_plugin.cpp     
2009-06-04 11:16:56 UTC (rev 20621)
+++ branches/ge_dev/source/blender/ikplugin/intern/itasc_plugin.cpp     
2009-06-04 12:47:59 UTC (rev 20622)
@@ -480,6 +480,7 @@
        ikscene->scene = scene;
        ikscene->cache = new iTaSC::Cache();;
        ikscene->solver = new iTaSC::WSDLSSolver();
+       //ikscene->solver->setQmax(10.0);
        ikscene->blArmature = ob;
 
        std::string  joint;
@@ -852,7 +853,7 @@
        if (reiterate) {
                // how many times do we reiterate?
                for (i=0; i<100; i++) {
-                       if (ikscene->armature->getMaxJointChange(timestep) < 
0.001)
+                       if (ikscene->armature->getMaxJointChange(timestep) < 
0.005)
                                break;
                        ikscene->scene->update(timestamp, timestep, 0, true, 
false);
                }
@@ -961,7 +962,7 @@
        arm = get_armature(ob);
        if (arm->ikdata) {
                IK_Data* ikdata = (IK_Data*)arm->ikdata;
-               for (IK_Scene* scene = ikdata->first; scene; scene = 
ikdata->first) {
+               for (IK_Scene* scene = ikdata->first; scene; scene = 
scene->next) {
                        if (scene->channels[0].pchan == pchan) {
                                execute_scene(scene, ctime);
                                break;


_______________________________________________
Bf-blender-cvs mailing list
Bf-blender-cvs@blender.org
http://lists.blender.org/mailman/listinfo/bf-blender-cvs

Reply via email to