1 /*
2  * This program is free software; you can redistribute it and/or
3  * modify it under the terms of the GNU General Public License
4  * as published by the Free Software Foundation; either version 2
5  * of the License, or (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program; if not, write to the Free Software Foundation,
14  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
15  */
16 
17 /** \file
18  * \ingroup freestyle
19  * \brief Class to perform all geometric operations dedicated to silhouette. That, for example,
20  * implies that this geom engine has as member data the viewpoint, transformations, projections...
21  */
22 
23 #include <cstdio>
24 #include <cstring>
25 
26 #include "Silhouette.h"
27 #include "SilhouetteGeomEngine.h"
28 
29 #include "../geometry/GeomUtils.h"
30 
31 #include "BKE_global.h"
32 
33 using namespace std;
34 
35 namespace Freestyle {
36 
37 Vec3r SilhouetteGeomEngine::_Viewpoint = Vec3r(0, 0, 0);
38 real SilhouetteGeomEngine::_translation[3] = {0, 0, 0};
39 real SilhouetteGeomEngine::_modelViewMatrix[4][4] = {
40     {1, 0, 0, 0},
41     {0, 1, 0, 0},
42     {0, 0, 1, 0},
43     {0, 0, 0, 1},
44 };
45 real SilhouetteGeomEngine::_projectionMatrix[4][4] = {
46     {1, 0, 0, 0},
47     {0, 1, 0, 0},
48     {0, 0, 1, 0},
49     {0, 0, 0, 1},
50 };
51 real SilhouetteGeomEngine::_transform[4][4] = {
52     {1, 0, 0, 0},
53     {0, 1, 0, 0},
54     {0, 0, 1, 0},
55     {0, 0, 0, 1},
56 };
57 int SilhouetteGeomEngine::_viewport[4] = {1, 1, 1, 1};
58 real SilhouetteGeomEngine::_Focal = 0.0;
59 
60 real SilhouetteGeomEngine::_glProjectionMatrix[4][4] = {
61     {1, 0, 0, 0},
62     {0, 1, 0, 0},
63     {0, 0, 1, 0},
64     {0, 0, 0, 1},
65 };
66 real SilhouetteGeomEngine::_glModelViewMatrix[4][4] = {
67     {1, 0, 0, 0},
68     {0, 1, 0, 0},
69     {0, 0, 1, 0},
70     {0, 0, 0, 1},
71 };
72 real SilhouetteGeomEngine::_znear = 0.0;
73 real SilhouetteGeomEngine::_zfar = 100.0;
74 bool SilhouetteGeomEngine::_isOrthographicProjection = false;
75 
76 SilhouetteGeomEngine *SilhouetteGeomEngine::_pInstance = NULL;
77 
setTransform(const real iModelViewMatrix[4][4],const real iProjectionMatrix[4][4],const int iViewport[4],real iFocal)78 void SilhouetteGeomEngine::setTransform(const real iModelViewMatrix[4][4],
79                                         const real iProjectionMatrix[4][4],
80                                         const int iViewport[4],
81                                         real iFocal)
82 {
83   unsigned int i, j;
84   _translation[0] = iModelViewMatrix[3][0];
85   _translation[1] = iModelViewMatrix[3][1];
86   _translation[2] = iModelViewMatrix[3][2];
87 
88   for (i = 0; i < 4; i++) {
89     for (j = 0; j < 4; j++) {
90       _modelViewMatrix[i][j] = iModelViewMatrix[j][i];
91       _glModelViewMatrix[i][j] = iModelViewMatrix[i][j];
92     }
93   }
94 
95   for (i = 0; i < 4; i++) {
96     for (j = 0; j < 4; j++) {
97       _projectionMatrix[i][j] = iProjectionMatrix[j][i];
98       _glProjectionMatrix[i][j] = iProjectionMatrix[i][j];
99     }
100   }
101 
102   for (i = 0; i < 4; i++) {
103     for (j = 0; j < 4; j++) {
104       _transform[i][j] = 0;
105       for (unsigned int k = 0; k < 4; k++) {
106         _transform[i][j] += _projectionMatrix[i][k] * _modelViewMatrix[k][j];
107       }
108     }
109   }
110 
111   for (i = 0; i < 4; i++) {
112     _viewport[i] = iViewport[i];
113   }
114   _Focal = iFocal;
115 
116   _isOrthographicProjection = (iProjectionMatrix[3][3] != 0.0);
117 }
118 
setFrustum(real iZNear,real iZFar)119 void SilhouetteGeomEngine::setFrustum(real iZNear, real iZFar)
120 {
121   _znear = iZNear;
122   _zfar = iZFar;
123 }
124 
retrieveViewport(int viewport[4])125 void SilhouetteGeomEngine::retrieveViewport(int viewport[4])
126 {
127   memcpy(viewport, _viewport, sizeof(int[4]));
128 }
129 
ProjectSilhouette(vector<SVertex * > & ioVertices)130 void SilhouetteGeomEngine::ProjectSilhouette(vector<SVertex *> &ioVertices)
131 {
132   Vec3r newPoint;
133   vector<SVertex *>::iterator sv, svend;
134   for (sv = ioVertices.begin(), svend = ioVertices.end(); sv != svend; sv++) {
135     GeomUtils::fromWorldToImage(
136         (*sv)->point3D(), newPoint, _modelViewMatrix, _projectionMatrix, _viewport);
137     (*sv)->setPoint2D(newPoint);
138   }
139 }
140 
ProjectSilhouette(SVertex * ioVertex)141 void SilhouetteGeomEngine::ProjectSilhouette(SVertex *ioVertex)
142 {
143   Vec3r newPoint;
144   GeomUtils::fromWorldToImage(
145       ioVertex->point3D(), newPoint, _modelViewMatrix, _projectionMatrix, _viewport);
146   ioVertex->setPoint2D(newPoint);
147 }
148 
ImageToWorldParameter(FEdge * fe,real t)149 real SilhouetteGeomEngine::ImageToWorldParameter(FEdge *fe, real t)
150 {
151   if (_isOrthographicProjection) {
152     return t;
153   }
154 
155   // we need to compute for each parameter t the corresponding parameter T which gives the
156   // intersection in 3D.
157   real T;
158 
159   // suffix w for world, c for camera, r for retina, i for image
160   Vec3r Aw = (fe)->vertexA()->point3D();
161   Vec3r Bw = (fe)->vertexB()->point3D();
162   Vec3r Ac, Bc;
163   GeomUtils::fromWorldToCamera(Aw, Ac, _modelViewMatrix);
164   GeomUtils::fromWorldToCamera(Bw, Bc, _modelViewMatrix);
165   Vec3r ABc = Bc - Ac;
166 #if 0
167   if (G.debug & G_DEBUG_FREESTYLE) {
168     cout << "Ac " << Ac << endl;
169     cout << "Bc " << Bc << endl;
170     cout << "ABc " << ABc << endl;
171   }
172 #endif
173   Vec3r Ai = (fe)->vertexA()->point2D();
174   Vec3r Bi = (fe)->vertexB()->point2D();
175   Vec3r Ii = Ai + t * (Bi - Ai);  // the intersection point in the 2D image space
176   Vec3r Ir, Ic;
177   GeomUtils::fromImageToRetina(Ii, Ir, _viewport);
178 
179   real alpha, beta, denom;
180   real m11 = _projectionMatrix[0][0];
181   real m13 = _projectionMatrix[0][2];
182   real m22 = _projectionMatrix[1][1];
183   real m23 = _projectionMatrix[1][2];
184 
185   if (fabs(ABc[0]) > 1.0e-6) {
186     alpha = ABc[2] / ABc[0];
187     beta = Ac[2] - alpha * Ac[0];
188     denom = alpha * (Ir[0] + m13) + m11;
189     if (fabs(denom) < 1.0e-6) {
190       goto iter;
191     }
192     Ic[0] = -beta * (Ir[0] + m13) / denom;
193 #if 0
194     Ic[1] = -(Ir[1] + m23) * (alpha * Ic[0] + beta) / m22;
195     Ic[2] = alpha * (Ic[0] - Ac[0]) + Ac[2];
196 #endif
197     T = (Ic[0] - Ac[0]) / ABc[0];
198   }
199   else if (fabs(ABc[1]) > 1.0e-6) {
200     alpha = ABc[2] / ABc[1];
201     beta = Ac[2] - alpha * Ac[1];
202     denom = alpha * (Ir[1] + m23) + m22;
203     if (fabs(denom) < 1.0e-6) {
204       goto iter;
205     }
206     Ic[1] = -beta * (Ir[1] + m23) / denom;
207 #if 0
208     Ic[0] = -(Ir[0] + m13) * (alpha * Ic[1] + beta) / m11;
209     Ic[2] = alpha * (Ic[1] - Ac[1]) + Ac[2];
210 #endif
211     T = (Ic[1] - Ac[1]) / ABc[1];
212   }
213   else {
214   iter:
215     bool x_coords, less_than;
216     if (fabs(Bi[0] - Ai[0]) > 1.0e-6) {
217       x_coords = true;
218       less_than = Ai[0] < Bi[0];
219     }
220     else {
221       x_coords = false;
222       less_than = Ai[1] < Bi[1];
223     }
224     Vec3r Pc, Pr, Pi;
225     real T_sta = 0.0;
226     real T_end = 1.0;
227     real delta_x, delta_y, dist, dist_threshold = 1.0e-6;
228     int i, max_iters = 100;
229     for (i = 0; i < max_iters; i++) {
230       T = T_sta + 0.5 * (T_end - T_sta);
231       Pc = Ac + T * ABc;
232       GeomUtils::fromCameraToRetina(Pc, Pr, _projectionMatrix);
233       GeomUtils::fromRetinaToImage(Pr, Pi, _viewport);
234       delta_x = Ii[0] - Pi[0];
235       delta_y = Ii[1] - Pi[1];
236       dist = sqrt(delta_x * delta_x + delta_y * delta_y);
237       if (dist < dist_threshold) {
238         break;
239       }
240       if (x_coords) {
241         if (less_than) {
242           if (Pi[0] < Ii[0]) {
243             T_sta = T;
244           }
245           else {
246             T_end = T;
247           }
248         }
249         else {
250           if (Pi[0] > Ii[0]) {
251             T_sta = T;
252           }
253           else {
254             T_end = T;
255           }
256         }
257       }
258       else {
259         if (less_than) {
260           if (Pi[1] < Ii[1]) {
261             T_sta = T;
262           }
263           else {
264             T_end = T;
265           }
266         }
267         else {
268           if (Pi[1] > Ii[1]) {
269             T_sta = T;
270           }
271           else {
272             T_end = T;
273           }
274         }
275       }
276     }
277 #if 0
278     if (G.debug & G_DEBUG_FREESTYLE) {
279       cout << "SilhouetteGeomEngine::ImageToWorldParameter(): #iters = " << i
280            << ", dist = " << dist << "\n";
281     }
282 #endif
283     if (i == max_iters && G.debug & G_DEBUG_FREESTYLE) {
284       cout << "SilhouetteGeomEngine::ImageToWorldParameter(): reached to max_iters (dist = "
285            << dist << ")\n";
286     }
287   }
288 
289   return T;
290 }
291 
WorldToImage(const Vec3r & M)292 Vec3r SilhouetteGeomEngine::WorldToImage(const Vec3r &M)
293 {
294   Vec3r newPoint;
295   GeomUtils::fromWorldToImage(M, newPoint, _transform, _viewport);
296   return newPoint;
297 }
298 
CameraToImage(const Vec3r & M)299 Vec3r SilhouetteGeomEngine::CameraToImage(const Vec3r &M)
300 {
301   Vec3r newPoint, p;
302   GeomUtils::fromCameraToRetina(M, p, _projectionMatrix);
303   GeomUtils::fromRetinaToImage(p, newPoint, _viewport);
304   return newPoint;
305 }
306 
307 } /* namespace Freestyle */
308