[renderer] mouseWheel zoom now centered on pointer location
[scilab.git] / scilab / modules / graphic_objects / src / cpp / pickSurface.cpp
1 /*
2  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3  * Copyright (C) 2012 - Pedro Arthur dos S. Souza
4  * Copyright (C) 2012 - Caio Lucas dos S. Souza
5  *
6  * Copyright (C) 2012 - 2016 - Scilab Enterprises
7  *
8  * This file is hereby licensed under the terms of the GNU GPL v2.0,
9  * pursuant to article 5.3.4 of the CeCILL v.2.1.
10  * This file was originally licensed under the terms of the CeCILL v2.1,
11  * and continues to be available under such terms.
12  * For more information, see the COPYING file which you should have received
13  * along with this program.
14  *
15  */
16
17 extern "C"
18 {
19 #include <stdio.h>
20 #include <math.h>
21 #include <string.h>
22
23 #include "getGraphicObjectProperty.h"
24 #include "graphicObjectProperties.h"
25 #include "sciprint.h"
26
27 double pickSurface(int uid, double x, double y,  double z, double dx, double dy, double dz);
28
29 }
30
31 #define EPS 1e-8
32
33 class Vec3
34 {
35 public:
36     double x, y, z;
37
38     Vec3(): x(0), y(0), z(0) {}
39     Vec3(double _x, double _y, double _z): x(_x), y(_y), z(_z) {}
40
41     Vec3 operator - (Vec3 v)
42     {
43         return Vec3( x - v.x, y - v.y, z - v.z );
44     }
45
46     Vec3 operator + (Vec3 v)
47     {
48         return Vec3( x + v.x, y + v.y, z + v.z );
49     }
50
51     Vec3 operator * (double s)
52     {
53         return Vec3( x * s, y * s, z * s );
54     }
55
56     Vec3 operator / (double s)
57     {
58         return Vec3( x / s, y / s, z / s );
59     }
60
61     double dot(Vec3 v)
62     {
63         return x * v.x + y * v.y + z * v.z;
64     }
65
66     Vec3 cross(Vec3 v)
67     {
68         return Vec3(y * v.z - z * v.y, z * v.x - x * v.z, x * v.y - y * v.x);
69     }
70
71     Vec3& normalize()
72     {
73         double d = sqrt(this->dot(*this));
74         if (d < EPS)
75         {
76             x = y = z = 0;
77         }
78         else
79         {
80             x /= d;
81             y /= d;
82             z /= d;
83         }
84         return *this;
85     }
86
87     Vec3 getNormalized()
88     {
89         Vec3 n = Vec3(x, y, z);
90         return n.normalize();
91     }
92
93     void print()
94     {
95         printf("\nv = %.8f, %.8f, %.8f", x, y, z);
96     }
97 };
98
99 int test_tri(Vec3 V1, Vec3 V2, Vec3 V3, Vec3 Dir, Vec3 P0, Vec3 &ret);
100 void QuadTestAndSaveT(double *bounds, Vec3 P0, Vec3 P1, Vec3 P2, Vec3 P3, Vec3 direction, Vec3 point, double &retT);
101
102 /*
103  * Given a ray (point(x, y,z) and direction(dx, dy, dz))
104  * check if the ray intersects any triangle from the given surface.
105  * returns t such that the intersection point is (x, y,z) + t * (dx, dy, dz)
106  * or -INFINITY if there isn't intersection.
107  */
108
109 double pickSurface(int uid, double x, double y,  double z, double dx, double dy, double dz)
110 {
111     double* X = NULL;
112     double* Y = NULL;
113     double* Z = NULL;
114
115     int type;
116     int * pType = &type;
117     double lastT = -INFINITY;
118
119     Vec3 direction = Vec3(dx, dy, dz);
120     Vec3 point = Vec3(x, y, z);
121
122     getGraphicObjectProperty(uid, __GO_DATA_MODEL_X__, jni_double_vector, (void**) &X);
123     getGraphicObjectProperty(uid, __GO_DATA_MODEL_Y__, jni_double_vector, (void**) &Y);
124     getGraphicObjectProperty(uid, __GO_DATA_MODEL_Z__, jni_double_vector, (void**) &Z);
125
126
127
128     int axes_uid = 0;
129     int * paxes_uid = &axes_uid;
130     int zoom_enabled = 0;
131     int *ptr = &zoom_enabled;
132     double *bounds;
133     getGraphicObjectProperty(uid, __GO_PARENT_AXES__, jni_int, (void**) &paxes_uid);
134     getGraphicObjectProperty(axes_uid, __GO_ZOOM_ENABLED__, jni_bool, (void**) &ptr);
135
136     if (zoom_enabled)
137     {
138         getGraphicObjectProperty(axes_uid, __GO_ZOOM_BOX__, jni_double_vector, (void**) &bounds);
139     }
140     else
141     {
142         getGraphicObjectProperty(axes_uid, __GO_DATA_BOUNDS__, jni_double_vector, (void**) &bounds);
143     }
144
145
146     getGraphicObjectProperty(uid, __GO_TYPE__, jni_int, (void**) &pType);
147     if (type == __GO_PLOT3D__)
148     {
149         int numX = 0;
150         int* piNumX = &numX;
151         int numY = 0;
152         int* piNumY = &numY;
153
154         getGraphicObjectProperty(uid, __GO_DATA_MODEL_NUM_X__, jni_int, (void**) &piNumX);
155         getGraphicObjectProperty(uid, __GO_DATA_MODEL_NUM_Y__, jni_int, (void**) &piNumY);
156
157         /* for each quad in the mesh separate it in 2 triangles
158          * and use test_tri function to test intersection from
159          * mouse click ray
160          * A point (x, y, z)  at (n,m) is given by
161          * (X[n], Y[m], Z[n][m]) where X, Y are vectors and Z a matrix.
162          */
163         for (int i = 0; i < (numX - 1); ++i)
164         {
165             for (int j = 0; j < (numY - 1); ++j)
166             {
167                 Vec3 P0 = Vec3(X[i],   Y[j],   Z[i + j * numX]);
168                 Vec3 P1 = Vec3(X[i + 1], Y[j],   Z[(i + 1) + j * numX]);
169                 Vec3 P2 = Vec3(X[i + 1], Y[j + 1], Z[(i + 1) + (j + 1) * numX]);
170                 Vec3 P3 = Vec3(X[i],   Y[j + 1], Z[i + (j + 1) * numX]);
171
172                 QuadTestAndSaveT(bounds, P0, P1, P2, P3, direction, point, lastT);
173             }
174         }
175     }
176     else if (type == __GO_FAC3D__)
177     {
178         int ng = 0, nvg = 0;
179         int *png = &ng, *pnvg = &nvg;
180         getGraphicObjectProperty(uid, __GO_DATA_MODEL_NUM_GONS__, jni_int, (void**) &png);
181         getGraphicObjectProperty(uid, __GO_DATA_MODEL_NUM_VERTICES_PER_GON__, jni_int, (void**) &pnvg);
182
183         if (nvg != 4)
184         {
185             return 2.0;
186         }
187
188         /*
189          * Fac3d data model is made by gons
190          * each gon should be a quad
191          * ordered in the vector
192          * X = [ p1, p2, p3, p4, p1, p2, p3, p4, ...]
193          * Y = [ p1, p2, p3, p4, p1, p2, p3, p4, ...]
194          * Z = [ p1, p2, p3, p4, p1, p2, p3, p4, ...]
195          * where a point is given by (x, y, z)
196          */
197         for (int i = 0; i < ng * nvg; i += nvg)
198         {
199             Vec3 P0 = Vec3(X[i],   Y[i],   Z[i]);
200             Vec3 P1 = Vec3(X[i + 1], Y[i + 1], Z[i + 1]);
201             Vec3 P2 = Vec3(X[i + 2], Y[i + 2], Z[i + 2]);
202             Vec3 P3 = Vec3(X[i + 3], Y[i + 3], Z[i + 3]);
203
204             QuadTestAndSaveT(bounds, P0, P1, P2, P3, direction, point, lastT);
205         }
206     }
207     if (type == __GO_GRAYPLOT__)
208     {
209
210         int numX = 0;
211         int* piNumX = &numX;
212         int numY = 0;
213         int* piNumY = &numY;
214
215         getGraphicObjectProperty(uid, __GO_DATA_MODEL_NUM_X__, jni_int, (void**) &piNumX);
216         getGraphicObjectProperty(uid, __GO_DATA_MODEL_NUM_Y__, jni_int, (void**) &piNumY);
217
218         /* Gray plot is a plane with Z = 0 and bounds = {x[0], x[n-1], y[0], y[m-1]}
219          * where n = size of vector x and m = size of vector y
220          */
221
222         Vec3 P0 = Vec3(X[0],      Y[0],      0);
223         Vec3 P1 = Vec3(X[numX - 1], Y[0],      0);
224         Vec3 P2 = Vec3(X[numX - 1], Y[numY - 1], 0);
225         Vec3 P3 = Vec3(X[0],      Y[numY - 1], 0);
226
227         QuadTestAndSaveT(bounds, P0, P1, P2, P3, direction, point, lastT);
228     }
229     if (type == __GO_MATPLOT__)
230     {
231         double* scale = NULL;
232         double* translate = NULL;
233         double zShift = 0.;
234         double* pdZShift = &zShift;
235         double mbounds[4];
236         int numX = 0;
237         int* piNumX = &numX;
238         int numY = 0;
239         int* piNumY = &numY;
240
241         getGraphicObjectProperty(uid, __GO_MATPLOT_SCALE__, jni_double_vector, (void**) &scale);
242         getGraphicObjectProperty(uid, __GO_MATPLOT_TRANSLATE__, jni_double_vector, (void**) &translate);
243         getGraphicObjectProperty(uid, __GO_DATA_MODEL_Z_COORDINATES_SHIFT__, jni_double, (void**) &pdZShift);
244         getGraphicObjectProperty(uid, __GO_DATA_MODEL_NUM_X__, jni_int, (void**) &piNumX);
245         getGraphicObjectProperty(uid, __GO_DATA_MODEL_NUM_Y__, jni_int, (void**) &piNumY);
246
247         mbounds[0] = translate[0];
248         mbounds[1] = translate[1];
249         mbounds[2] = mbounds[0] + (numX - 1) * scale[0];
250         mbounds[3] = mbounds[1] + (numY - 1) * scale[1];
251
252         Vec3 P0 = Vec3(mbounds[0],      mbounds[1],      zShift);
253         Vec3 P1 = Vec3(mbounds[2],              mbounds[1],      zShift);
254         Vec3 P2 = Vec3(mbounds[2],      mbounds[3],      zShift);
255         Vec3 P3 = Vec3(mbounds[0],      mbounds[3],      zShift);
256
257         QuadTestAndSaveT(bounds, P0, P1, P2, P3, direction, point, lastT);
258     }
259     return lastT;
260
261 }
262
263 /*
264  * Test if the ray intersects the given triangle (V1, V2, V3)
265  * Fast, minimum storage ray/triangle intersection
266  * algorithm propose by Tomas Möller and Ben Trumbore.
267  * Calculate barycentric cordinates (u, v) and test if
268  * 0 <= u <= 1 && 0 <= v <= 1 && (u + v) <= 1, then the
269  * intersection point is inside the triangle.
270  */
271 int test_tri(Vec3 V1, Vec3 V2, Vec3 V3, Vec3 Dir, Vec3 P0, double &t)
272 {
273     Vec3 Edge1, Edge2, tVec, pVec, qVec;
274     double det, inv_det, u, v;
275
276     Edge1 = V2 - V1;
277     Edge2 = V3 - V1;
278
279     pVec = Dir.cross(Edge2);
280     det = Edge1.dot(pVec);
281
282     if (det > -EPS && det < EPS)
283     {
284         return 0;
285     }
286
287     inv_det = 1 / det;
288
289     tVec = P0 - V1;
290     u = tVec.dot(pVec) * inv_det;
291
292     if (u < 0.0 || u > 1.0)
293     {
294         return 0;
295     }
296
297     qVec = tVec.cross(Edge1);
298     v = Dir.dot(qVec) * inv_det;
299
300     if (v < 0.0 || (u + v) > 1.0)
301     {
302         return 0;
303     }
304
305     t = Edge2.dot(qVec) * inv_det;
306
307     return 1;
308 }
309
310 bool isInViewBox(double * bounds, Vec3 point)
311 {
312     return (bounds[0] <= point.x && bounds[1] >= point.x &&
313             bounds[2] <= point.y && bounds[3] >= point.y &&
314             bounds[4] <= point.z && bounds[5] >= point.z);
315 }
316
317 void QuadTestAndSaveT(double *bounds, Vec3 P0, Vec3 P1, Vec3 P2, Vec3 P3, Vec3 direction, Vec3 point, double &retT)
318 {
319     Vec3 intersectionPoint;
320     double t;
321     /*test first triangle*/
322     if (test_tri(P0, P1, P2, direction, point, t) == 1)
323     {
324         /*the intersection point can be outside the view box(invisible)*/
325         intersectionPoint = point + direction * t;        
326         if (isInViewBox(bounds, intersectionPoint))
327         {
328             // ray intersects the triangle, we store the greatest t (nearest from viewpoint) 
329             retT = retT > t ? retT : t;
330         }
331     }
332     /*test second triangle*/
333     if (test_tri(P0, P2, P3, direction, point, t) == 1)
334     {
335         intersectionPoint = point + direction * t;        
336         if (isInViewBox(bounds, intersectionPoint))
337         {
338             retT = retT > t ? retT : t;
339         }
340     }
341 }
342
343
344
345