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
  26     double pickSurface(int uid, double x, double y,  double z, double dx, double dy, double dz, double mx, double my, double mz, double mw);
  27
  28 }
  29
  30 #define EPS 1e-8
  31
  32 class Vec3
  33 {
  34 public:
  35     double x, y, z;
  36
  37     Vec3(): x(0), y(0), z(0) {}
  38     Vec3(double _x, double _y, double _z): x(_x), y(_y), z(_z) {}
  39
  40     Vec3 operator - (Vec3 v)
  41     {
  42         return Vec3( x - v.x, y - v.y, z - v.z );
  43     }
  44
  45     Vec3 operator + (Vec3 v)
  46     {
  47         return Vec3( x + v.x, y + v.y, z + v.z );
  48     }
  49
  50     Vec3 operator * (double s)
  51     {
  52         return Vec3( x * s, y * s, z * s );
  53     }
  54
  55     Vec3 operator / (double s)
  56     {
  57         return Vec3( x / s, y / s, z / s );
  58     }
  59
  60     double dot(Vec3 v)
  61     {
  62         return x * v.x + y * v.y + z * v.z;
  63     }
  64
  65     Vec3 cross(Vec3 v)
  66     {
  67         return Vec3(y * v.z - z * v.y, z * v.x - x * v.z, x * v.y - y * v.x);
  68     }
  69
  70     Vec3& normalize()
  71     {
  72         double d = sqrt(this->dot(*this));
  73         if (d < EPS)
  74         {
  75             x = y = z = 0;
  76         }
  77         else
  78         {
  79             x /= d;
  80             y /= d;
  81             z /= d;
  82         }
  83         return *this;
  84     }
  85
  86     Vec3 getNormalized()
  87     {
  88         Vec3 n = Vec3(x, y, z);
  89         return n.normalize();
  90     }
  91
  92     void print()
  93     {
  94         printf("\nv = %.8f, %.8f, %.8f", x, y, z);
  95     }
  96 };
  97
  98 int test_tri(Vec3 V1, Vec3 V2, Vec3 V3, Vec3 Dir, Vec3 P0, Vec3 &ret);
  99 void QuadTestAndSaveZ(double *bounds, Vec3 P0, Vec3 P1, Vec3 P2, Vec3 P3, Vec3 direction, Vec3 point,
 100                       double mx, double my, double mz, double mw, double &retZ);
 101
 102 /*
 103  * Given a ray (point(x, y,z) + direction(dx, dy, dz))
 104  * check if the ray intersects any triangle from the given surface.
 105  * returns the projected Z from the intersection point (p) (p.x*mx + p.y*my + p.z*mz + mw;)
 106  * or 2.0 if there isn't intersection (projected z vary between -1.0 - 1.0).
 107  */
 108
 109 double pickSurface(int uid, double x, double y,  double z, double dx, double dy, double dz, double mx, double my, double mz, double mw)
 110 {
 111     double* X = NULL;
 112     double* Y = NULL;
 113     double* Z = NULL;
 114
 115     int type;
 116     int * pType = &type;
 117     double lastZ = 2.0;
 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                 QuadTestAndSaveZ(bounds, P0, P1, P2, P3, direction, point, mx, my, mz, mw, lastZ);
 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             QuadTestAndSaveZ(bounds, P0, P1, P2, P3, direction, point, mx, my, mz, mw, lastZ);
 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         QuadTestAndSaveZ(bounds, P0, P1, P2, P3, direction, point, mx, my, mz, mw, lastZ);
 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         QuadTestAndSaveZ(bounds, P0, P1, P2, P3, direction, point, mx, my, mz, mw, lastZ);
 258     }
 259     return lastZ;
 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, Vec3 &ret)
 272 {
 273     Vec3 Edge1, Edge2, tVec, pVec, qVec;
 274     double det, inv_det, t, 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     ret = P0 + Dir * t;
 307
 308
 309     return 1;
 310 }
 311
 312 bool isInViewBox(double * bounds, Vec3 point)
 313 {
 314     return (bounds[0] <= point.x && bounds[1] >= point.x &&
 315             bounds[2] <= point.y && bounds[3] >= point.y &&
 316             bounds[4] <= point.z && bounds[5] >= point.z);
 317 }
 318
 319 void QuadTestAndSaveZ(double *bounds, Vec3 P0, Vec3 P1, Vec3 P2, Vec3 P3, Vec3 direction, Vec3 point,
 320                       double mx, double my, double mz, double mw, double &retZ)
 321 {
 322     Vec3 ret;
 323
 324     /*test first triangle*/
 325     if (test_tri(P0, P1, P2, direction, point, ret) == 1)
 326     {
 327         /*the intersection point can be outside the view box(invisible)*/
 328         if (isInViewBox(bounds, ret))
 329         {
 330             /* ray intersects the triangle, then we project only the Z cordinate
 331              * and store the nearest projected Z.
 332              */
 333             double curZ = ret.x * mx + ret.y * my + ret.z * mz + mw;
 334             retZ = retZ < curZ ? retZ : curZ;
 335         }
 336     }
 337     /*test second triangle*/
 338     if (test_tri(P0, P2, P3, direction, point, ret) == 1)
 339     {
 340         if (isInViewBox(bounds, ret))
 341         {
 342             double curZ = ret.x * mx + ret.y * my + ret.z * mz + mw;
 343             retZ = retZ < curZ ? retZ : curZ;
 344         }
 345     }
 346 }
 347
 348
 349
 350