Fixed a bug in SilhouetteGeomEngine::ImageToWorldParameter() that caused
authorTamito Kajiyama <rd6t-kjym@asahi-net.or.jp>
Sun, 10 Jan 2010 14:08:59 +0000 (14:08 +0000)
committerTamito Kajiyama <rd6t-kjym@asahi-net.or.jp>
Sun, 10 Jan 2010 14:08:59 +0000 (14:08 +0000)
instability issues regarding the view map creation.  A new iterative
solver of the 2D-to-3D inverse projection transformation problem was
implemented.  Instead of directly solving the problem in the direction
from the 2D to 3D space, the new solver starts with an initial guess of
an approximated solution and asymptotically approaches to the true
solution by iteratively performing the forward 3D-to-2D projection
transformation and improving the approximation.  Preliminary tests with
one simple and another complex scenes showed that the solver converges
quickly (more and less 20 iterations in many cases, with a stopping
criterion of a residual distance between the true and approximated
solutions less than 1e-6 Blender Unit).

source/blender/freestyle/intern/blender_interface/BlenderFileLoader.cpp
source/blender/freestyle/intern/blender_interface/BlenderFileLoader.h
source/blender/freestyle/intern/blender_interface/FRS_freestyle.cpp
source/blender/freestyle/intern/view_map/SilhouetteGeomEngine.cpp
source/blender/freestyle/intern/view_map/ViewMapBuilder.cpp

index 114dd1bc7e2265f7d13f1c62c918a2e566239e35..6efa80431fc4895bbdd07195b313b22a4b9ea740 100644 (file)
@@ -21,8 +21,20 @@ NodeGroup* BlenderFileLoader::Load()
 
        cout << "\n===  Importing triangular meshes into Blender  ===" << endl;
 
-  // creation of the scene root node
-  _Scene = new NodeGroup;
+       // creation of the scene root node
+       _Scene = new NodeGroup;
+
+       _viewplane_left=   _re->viewplane.xmin;
+       _viewplane_right=  _re->viewplane.xmax;
+       _viewplane_bottom= _re->viewplane.ymin;
+       _viewplane_top=    _re->viewplane.ymax;
+       _z_near= _re->clipsta;
+       _z_far=  _re->clipend;
+#if 0
+       cout << "frustrum: l " << _viewplane_left << " r " << _viewplane_right
+               << " b " << _viewplane_bottom << " t " << _viewplane_top
+               << " n " << _z_near << " f " << _z_far << endl;
+#endif
 
        int id = 0;
        for(obi= (ObjectInstanceRen *) _re->instancetable.first; obi; obi=obi->next) {
@@ -37,8 +49,8 @@ NodeGroup* BlenderFileLoader::Load()
                        cout << "  Sorry, only vlak-based shapes are supported." << endl;
        }
 
-  //Returns the built scene.
-  return _Scene;
+       //Returns the built scene.
+       return _Scene;
 }
 
 void BlenderFileLoader::insertShapeNode(ObjectRen *obr, int id)
index 8068eaf1bcf37db0d3695ea7dad036f11298e746..29acbfc203ffaeed4b9b06e2408777c4da3ff87b 100644 (file)
@@ -55,9 +55,14 @@ protected:
 protected:
        Render* _re;
        SceneRenderLayer* _srl;
-  NodeGroup* _Scene;
-  unsigned _numFacesRead;
-  real _minEdgeSize;
+       NodeGroup* _Scene;
+       unsigned _numFacesRead;
+       real _minEdgeSize;
+       float _viewplane_left;
+       float _viewplane_right;
+       float _viewplane_bottom;
+       float _viewplane_top;
+       float _z_near, _z_far;
 };
 
 #endif // BLENDER_FILE_LOADER_H
index 393cd1b2dd9f32f442932bc8c25e924fdd1f3538..4ef3abe54d851d8f58e72ea4800b01b0fff9a394 100644 (file)
@@ -19,6 +19,7 @@ extern "C" {
 
 #include "BKE_main.h"
 #include "BLI_blenlib.h"
+#include "BLI_math.h"
 #include "BPY_extern.h"
 
 #include "renderpipeline.h"
@@ -106,6 +107,9 @@ extern "C" {
                for( int i = 0; i < 4; i++ )
                   for( int j = 0; j < 4; j++ )
                        freestyle_proj[i][j] = re->winmat[i][j];
+
+               //print_m4("mv", freestyle_mv);
+               //print_m4("proj", freestyle_proj);
        }
 
        
index ecb0a4979f0b7a7832417e967c8ba6870862d679..bb9ef5d574b50868ffd4c9ce84ae1c2f51fd7b89 100755 (executable)
@@ -155,6 +155,7 @@ real SilhouetteGeomEngine::ImageToWorldParameter(FEdge *fe, real t)
        
   // we need to compute for each parameter t the corresponding 
   // parameter T which gives the intersection in 3D.
+#if 0
   //currentEdge = (*fe);
   Vec3r A = (fe)->vertexA()->point3D();
   Vec3r B = (fe)->vertexB()->point3D();
@@ -175,6 +176,46 @@ real SilhouetteGeomEngine::ImageToWorldParameter(FEdge *fe, real t)
   
   real T;
   T = (Ic[2]*Ac[1] - Ic[1]*Ac[2])/(Ic[1]*(Bc[2]-Ac[2])-Ic[2]*(Bc[1]-Ac[1]));
+#else
+       // suffix w for world, i for image
+       Vec3r Aw = (fe)->vertexA()->point3D();
+       Vec3r Bw = (fe)->vertexB()->point3D();
+       Vec3r Ai = (fe)->vertexA()->point2D();
+       Vec3r Bi = (fe)->vertexB()->point2D();
+       Vec3r Ii = Ai + t * (Bi - Ai); // the intersection point in 2D
+       Vec3r Pw, Pi;
+       real T_sta = 0.0;
+       real T_end = 1.0;
+       real T;
+       real delta_x, delta_y, dist, dist_threshold = 1e-6;
+       int i, max_iters = 100;
+       for (i = 0; i < max_iters; i++) {
+        T = T_sta + 0.5 * (T_end - T_sta);
+        Pw = Aw + T * (Bw - Aw);
+               GeomUtils::fromWorldToImage(Pw, Pi, _transform, _viewport);
+               delta_x = Ii[0] - Pi[0];
+               delta_y = Ii[1] - Pi[1];
+        dist = sqrt(delta_x * delta_x + delta_y * delta_y);
+        if (dist < dist_threshold)
+            break;
+               if (Ai[0] < Bi[0]) {
+            if (Pi[0] < Ii[0])
+                T_sta = T;
+            else
+                T_end = T;
+               } else {
+            if (Pi[0] > Ii[0])
+                T_sta = T;
+            else
+                T_end = T;
+               }
+       }
+#if 0
+       printf("SilhouetteGeomEngine::ImageToWorldParameter(): #iters = %d, dist = %e\n", i, dist);
+#endif
+       if (i == max_iters)
+               printf("SilhouetteGeomEngine::ImageToWorldParameter(): reached to max_iters (dist = %e)\n", dist);
+#endif
   
   return T;
 }
index 3a4d86ea16f4900bb9f18ad6159d30430e56d50e..9c7a084c2047ee3bfe226c4b3d9be19ac66990b5 100755 (executable)
@@ -946,14 +946,20 @@ void ViewMapBuilder::ComputeSweepLineIntersections(ViewMap *ioViewMap, real epsi
     real tb = (*i)->tB;
 
     if((ta < -epsilon) || (ta > 1+epsilon))
-        cerr << "Warning: intersection out of range for edge " << fA->vertexA()->getId() << " - " << fA->vertexB()->getId() << endl;
+        cerr << "Warning: 2D intersection out of range for edge " << fA->vertexA()->getId() << " - " << fA->vertexB()->getId() << endl;
     
     if((tb < -epsilon) || (tb > 1+epsilon))
-        cerr << "Warning: intersection out of range for edge " << fB->vertexA()->getId() << " - " << fB->vertexB()->getId() << endl;
+        cerr << "Warning: 2D intersection out of range for edge " << fB->vertexA()->getId() << " - " << fB->vertexB()->getId() << endl;
     
     real Ta = SilhouetteGeomEngine::ImageToWorldParameter(fA, ta);
     real Tb = SilhouetteGeomEngine::ImageToWorldParameter(fB, tb);
 
+    if((Ta < -epsilon) || (Ta > 1+epsilon))
+        cerr << "Warning: 3D intersection out of range for edge " << fA->vertexA()->getId() << " - " << fA->vertexB()->getId() << endl;
+    
+    if((Tb < -epsilon) || (Tb > 1+epsilon))
+        cerr << "Warning: 3D intersection out of range for edge " << fB->vertexA()->getId() << " - " << fB->vertexB()->getId() << endl;
+
     TVertex * tvertex = ioViewMap->CreateTVertex(Vec3r(A1 + Ta*(A2-A1)), Vec3r(a1 + ta*(a2-a1)), fA, 
                                                  Vec3r(B1 + Tb*(B2-B1)), Vec3r(b1 + tb*(b2-b1)), fB, id);