version 1.2.0.0 (257 KB) by
Vipin Vijayan

Ray casting or collision detection for deformable triangular 3D meshes.

This is a Matlab wrapper for OPCODE, which is a collision detection

or ray casting library for triangular 3D meshes.

OPCODE uses a couple of different aabb trees to store the mesh and this is

a pretty simple wrapper for one of the trees.

Nice thing about opcode is that it allows for deformable meshes,

meaning that you can update the mesh while it is stored in the tree,

which is much faster than what it takes to rebuild the aabb tree.

Input and output:

To make the tree:

tree = opcodemesh(v,f);

where

vertices v : 3 x nv

faces f : 3 x nf

To intersect:

[hit,d,trix,bary,Q] = tree.intersect(orig,dir);

where

starting points orig : 3 x nc

direction dir : 3 x nc

whether hit or not hit : nc x 1 logical

distance from orig to intersection point d : nc x 1

index into f of the intersection triangle trix : nc x 1

barycentric coordinates of the triangle that the rays intersected bary : 2 x nc

actual 3D coordinates of the intersection poitns Q : 3 x nc

If a ray misses, then you have 0's for trix and nan's for d,bary,Q at that index.

To get the actual 3D coords, you can also do

Qhit = repmat(1-B(2,:)-B(1,:),3,1).*v(:,f(1,T)) + ...

repmat(B(1,:),3,1).*v(:,f(2,T)) + ...

repmat(B(2,:),3,1).*v(:,f(3,T));

where B = bary(:,hit), T = trix(hit),

which gives you the coordinates at the intersected points.

To update the mesh with a new set of vertices (as in deform the mesh):

tree.update(vnew);

vnew : 3 x nv

vnew must be the same size and the original set of vertices.

The topology of the mesh cannot change.

To delete the tree (which is actually done automatically by Matlab when the variable is cleared):

tree.delete();

To compile, go to ./src_matlab

and run mexall.m

which compiles the mex code and copies it to ./matlab

(I compile it against everything in opcode, so it's a bit slow.)

To run the demo, go to ./matlab

and run opcodemeshdemo.m

OPCODE is in

http://www.lia.ufc.br/~gilvan/cd/

(which is a more portable version)

and

http://www.codercorner.com/Opcode.htm

(the original site)

Contains much appreciated code from

http://www.mathworks.com/matlabcentral/fileexchange/38964

for nicely wrapping C++ functions in a Matlab class.

Fixes from the previous version:

Fixed compile issues for Win64 by removing all of the assembly code.

You don't need to normalize the direction vectors anymore.

It also returns the actual points now instead of just the barycentric points.

Vipin Vijayan (2021). Ray casting for deformable triangular 3D meshes (https://www.mathworks.com/matlabcentral/fileexchange/41504-ray-casting-for-deformable-triangular-3d-meshes), MATLAB Central File Exchange. Retrieved .

Created with
R2013a

Compatible with any release

**Inspired by:**
Example MATLAB class wrapper for a C++ class

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!Create scripts with code, output, and formatted text in a single executable document.

Carmelo MineoHello @Vipin Vijayan and all.

I found this code very useful. I am using it to find the closest intersections between rays and a triangular mesh.

The code works well and I am impressed by the computation speed. However, I found approximation problems that cause some intersections to be missed. I had a look at the opcode source files and realized it converts the mesh vertices and the ray vectors into float numbers and does the computation in single precision. I am wondering if you anyone has realized this and has edited the code to create a double-precision version of opcode.

Cheers,

CarMineo

Austin FiteJun WWhat is 'direction dir : 3 x nc'?

Chao WangTo help others, I put my building experience here. I just compiled it in matlab 2018b 64bit in win10, and only got one error:

------

C:\Program Files (x86)\Microsoft Visual Studio 14.0\VC\INCLUDE\xkeycheck.h(246): warning C4005: 'override': macro redefinition

c:\dev\kinectir\opcodemesh\opcode\Ice/IcePreprocessor.h(91): note: see previous definition of 'override'

C:\Program Files (x86)\Microsoft Visual Studio 14.0\VC\INCLUDE\xkeycheck.h(250): fatal error C1189: #error: The C++ Standard

Library forbids macroizing keywords. Enable warning C4005 to find the forbidden macro

------

This is actually because the code redefines C++ keyword 'override' in line 91 in IcePreprocessor.h. To fix it, I simply replace all the 'override' in ALL source code files in folder 'opcode' to some other words such as 'overrides'. I can build the code well after fixing this.

mohamed dahymohamed dahyThanks @ThT for your comment it helped me so much, I just want to add additional information to it as I faced other problems in Win10:

1) Install SDK 7.1 from http://www.microsoft.com/en-us/download/details.aspx?id=8279

>However, it fails to detect .NET framework 4. >> so continue the installation process as it is without the c++ compiler.

2) Then you can download c++ compiler "VC-Compiler-KB2519277.exe" manually at the following link: http://www.microsoft.com/en-us/download/details.aspx?id=4422 .

3)configure matlab to use this compiler in matlab by typing

"mex -setup C++" in matlab command window

then press the "Microsoft Windows SDK 7.1 (C++) "

hope it is useful for someone as well.

ThT@Mohammad Alomari try to compile the lib with the winsdk-7.1 (windows sdk 7.1) compiler instead. That worked for me. The point is that there are some declarations in the opcode source code that in the new compilers are existing and that causes a conflict. Unfortunately the opcode library is not updated quite some years now.

Mohammad Alomari@Ngo Nguyen @ThT

I have the same error

fatal error C1189: #error:

The C++ Standard Library forbids macroizing keywords. Enable

warning C4005 to find the forbidden macro.

can anyone help with compiling.

Thanks

ThT@Ngo Nguyen I am getting the same error, did you manage to find a solution for it? @Philipp Ertl can I ask you how did you compile the lib in your case.

Thanks.

Philipp Ertl@Stefan Spelitz, @Vipin Vijayan. I think i found a solution for the bug. Although it might have created other bugs :D

use the following code in OPC_RayTriOverlap.h

https://pastebin.com/VeFhns4W

using the Example code of Stefan Spelitz the following results are calculated:

hit1 = true

hit2 = true

cheers.

ThT@Stefan Spelitz, since the bug is on the original opcode you could contact the developers there given the links from Vipin Vijayan above.

Gaoli Sanghello,

I am really appreciate your good job. I will help me a lot.

now, the first link you provided is invalid. so could you kindly send the code to my email? g.sang@utwente.nl. Thank you very much.

Stefan SpelitzThank you for looking into it anyway!

Vipin VijayanThanks for describing the bug. I looked into it and it seems to be a bug in the original opcode. Unfortunately I won't be fixing it since it is too much effort.

Stefan SpelitzI found a bug! Not all intersections are correctly detected. It seems to occur when shooting a ray orthogonal to the edge of a triangle. It depends on the winding of the triangles and/or side of the ray.

I cannot tell if the bug is already in the original opcode or in the Matlab bridge.

Example code:

f = [2 4 1 ; 2 3 4]';

v = [-1 -1 0 ; -1 1 0 ; 1 1 0 ; 1 -1 0]';

t = opcodemesh(v,f);

p1 = [0.0 0.0 -1]';

p2 = [0.0 0.0 1]';

[hit1,d1] = t.intersect(p1,(p2-p1)./norm(p2-p1)); % ok

[hit2,d2] = t.intersect(p2,(p1-p2)./norm(p1-p2)); % NOK

% Visualization:

trimesh(f',v(1,:),v(2,:),v(3,:));

hold on;

xlabel('x');

ylabel('y');

zlabel('z');

axis equal;

xlim([-2 2]);

ylim([-2 2]);

zlim([-2 2]);

plot3([p1(1) p2(1)], [p1(2) p2(2)], [p1(3) p2(3)], '-ro');

hold off;

Vladislav KramarevAllaNgo NguyenHi Vipin, I got the error:

C:\Program Files (x86)\Microsoft Visual Studio

14.0\VC\INCLUDE\xkeycheck.h(250): fatal error C1189: #error:

The C++ Standard Library forbids macroizing keywords. Enable

warning C4005 to find the forbidden macro.

Can you tell me how to fix it?

Thank you.

Johannes KorsaweReally good.

One small thing to improve (could be much effort, though):

Also return multiple intersections, not only the first.

TCHHere I put a test code but I cannot get the hit result that obviously shown in the figure(1). I am totally confused with the result, is there a bug in the original opcode?

@@@ code start

clc

clear

faces = [ ...

10 , 6 , 4 ; ...

11 , 6 , 10 ; ...

6 , 7 , 4 ; ...

10 , 4 , 13 ; ...

11 , 8 , 6 ; ...

16 , 11 , 10 ; ...

6 , 3 , 7 ; ...

4 , 7 , 2 ; ...

13 , 4 , 5 ; ...

14 , 10 , 13 ; ...

12 , 8 , 11 ; ...

8 , 3 , 6 ; ...

17 , 11 , 16 ; ...

16 , 10 , 14 ; ...

3 , 1 , 7 ; ...

7 , 0 , 2 ; ...

4 , 2 , 5 ; ...

13 , 5 , 9 ; ...

14 , 13 , 15 ...

] + 1;

vertices = [ ...

4.231401060000e-001 , 1.289017940000e-001 , 8.522257000000e-003 ; ...

4.234843750000e-001 , 1.295909270000e-001 , 8.482499000000e-003 ; ...

4.203841550000e-001 , 1.286728970000e-001 , 9.441268000000e-003 ; ...

4.212358700000e-001 , 1.298959500000e-001 , 9.275906000000e-003 ; ...

4.181593630000e-001 , 1.289744720000e-001 , 1.030612100000e-002 ; ...

4.176608580000e-001 , 1.284392850000e-001 , 1.045896700000e-002 ; ...

4.186631770000e-001 , 1.295075680000e-001 , 1.016991900000e-002 ; ...

4.208076170000e-001 , 1.292857210000e-001 , 9.348370000000e-003 ; ...

4.191720280000e-001 , 1.300382540000e-001 , 1.004999700000e-002 ; ...

4.151989140000e-001 , 1.282018740000e-001 , 1.154331100000e-002 ; ...

4.163471070000e-001 , 1.291179810000e-001 , 1.115785700000e-002 ; ...

4.169212040000e-001 , 1.295760350000e-001 , 1.096513000000e-002 ; ...

4.174953000000e-001 , 1.300340880000e-001 , 1.077240300000e-002 ; ...

4.157730100000e-001 , 1.286599270000e-001 , 1.135058400000e-002 ; ...

4.162002870000e-001 , 1.305507970000e-001 , 1.134597100000e-002 ; ...

4.156481020000e-001 , 1.302543030000e-001 , 1.150225100000e-002 ; ...

4.167558590000e-001 , 1.308455510000e-001 , 1.119693700000e-002 ; ...

4.173129270000e-001 , 1.311391140000e-001 , 1.105040500000e-002 ...

];

x = vertices(:,1);

y = vertices(:,2);

z = vertices(:,3);

orig0 = [0.41749369461334546, 0.12947887360062926, 0.010681338615420484]; % ray's origin

dir0 = [-0.40977066483612395, 0.15771530547732460, -0.89845082484126815]; % ray's direction

orig = orig0 - dir0 * 0.01;

orig0 = orig;

orig = [orig;orig];

dir = [dir0;-dir0];

endp = orig + dir * 0.1;

linex = [orig(:,1), endp(:,1)];

liney = [orig(:,2), endp(:,2)];

linez = [orig(:,3), endp(:,3)];

t = opcodemesh(vertices',faces');

tic;

[hit,d,trix,bary,Q] = t.intersect(orig',dir');

% intersect = TriangleRayIntersection(orig, dir, vert1, vert2, vert3);

fprintf('Number of: faces=%i, points=%i, intresections=%i; time=%f sec\n', ...

size(faces,1), size(vertices,1), sum(hit), toc);

%%

% *Display the results: Surface in blue, line in light read and intersected

% triangles in dark red*

figure(1); clf;

hold on;

trisurf(faces,x,y,z,'FaceAlpha', 0.9)

hold on;

plot3(linex,liney,linez,'-k', 'linewidth', 2);

hold on

plot3(orig0(1),orig0(2),orig0(3), '-bo');

% set(gca, 'CameraPosition', [106.2478 -35.9079 136.4875])

%set(gco,'EdgeColor','none');

xlabel('x');

ylabel('y');

zlabel('z');

axis tight

@@@ code end

Vipin VijayanI doubt it. I've never used the mesh to mesh collision function and I'm not really using this code much anymore.

AndreAny chance to get the mesh to mesh collision function wrapped?

EricGreat submission. Very useful.

Thank you Germano Gomes for mex64 advice.

Vipin VijayanYeah, sorry about the confusion with norm(to-from). You're right, normc should be used for the input.

Germano GomesNote(ray casting): the correct value of 'd' is obtained by normalizing each column of the (to-from) 3xN matrix. I use Matlab normc for this:

[hit,d,trix,bary] = t.intersect(from,normc(to-from));

do not use:

[hit,d,trix,bary] = t.intersect(from,(to-from)./norm(to-from));

if you have more than one vector/ray

Germano GomesVery useful Matlab wrapper for fast collision detection and ray tracing.

For those using Matlab win64: look at the error information and comment out (//) the machine code (the __asm keyword invokes the inline assembler).

Example solution: error fix in IceFPU.h line 47:

//! Fast square root for floating-point values.

inline_ float FastSqrt(float square)

{

// #ifdef _MSC_VER

// float retval;

//

// __asm {

// mov eax, square

// sub eax, 0x3F800000

// sar eax, 1

// add eax, 0x3F800000

// mov [retval], eax

// }

// return retval;

// #else

return sqrt(square);

// #endif

In all cases just leave the code below #else (this is in vanilla C++) that does the same work without any incompatibility with win64.

Finally, unfortunately the returned d (matrix of distances) returns wrong values that seem to be scaled or shifted by the number of rays (vectors from-to), even if this vector is normalized as suggested by Vipin.

Vipin can you check and fix this please.

For it to be perfect I suggest returning the correct 'd' matrix, not needing the normalization step in Matlab and returning the intersection points(Q matrix) from the Mex function. That would be nice.

Thank you Vipin

Germano GomesVipin VijayanI don't have an Matlab win64 install now. I tried it on win 32, os x 64, and linux 64. Maybe you can tell me what the error looks like and I will see what I can do.

omarThanks Vipin..now i've the same problem with K G...i get errors when compiling it under win 64 !

K Ghow can I use it in Matlab(win64), it can not be mex in a 64bit system

Vipin VijayanThe compiled file is not in your path or you didn't compile it.

omarHi...i get an error when i try to run the opcodemeshdemo.m in Matlab(2012b-64).

//

Undefined function 'opcodemeshmex' for input

arguments of type 'char'.

Error in opcodemesh (line 8)

this.objectHandle =

opcodemeshmex('create',

varargin{:});

Error in opcodemeshdemo (line 19)

t = opcodemesh(v',f');

//

any help ?

Thanks,

omar

Vipin VijayanMade a mistake in the previous comment. Calculation of Q should actually be like the following example.

f = [1 2 3]'; v = [-1 0 0.1 ; 1 1 0.1 ; 1 -1 0.1]';

t = opcodemesh(v,f);

from = [0.17 0.17 0.5]';

to = [0.17 0.17 -1]';

[hit,d,trix,bary] = t.intersect(from,(to-from)./norm(to-from));

B = bary(:,hit); T = trix(hit);

Q = repmat(1-B(2,:)-B(1,:),3,1).*v(:,f(1,T)) + repmat(B(1,:),3,1).*v(:,f(2,T)) + repmat(B(2,:),3,1).*v(:,f(3,T));

Vipin VijayanYou had to normalize the direction (to - from) and I had gotten the calculation for the points Q wrong. I will upload a version eventually so you don't have to normalize the direction. Btw, here is an example that works:

f = [1 2 3]'; v = [-1 0 0.1 ; 1 1 0.1 ; 1 -1 0.1]';

t = opcodemesh(v,f);

from = [0 0 0.5]';

to = [0 0 -1]';

[hit,d,trix,bary] = t.intersect(from,(to-from)./norm(to-from));

B = bary(:,hit); T = trix(hit);

Q = repmat(B(1,:),3,1).*v(:,f(3,T)) + repmat(B(2,:),3,1).*v(:,f(2,T)) + repmat(1-B(2,:)-B(1,:),3,1).*v(:,f(1,T));

ZhouI used the code following the descriptions.

i have one questions:

(1) why d is not equal to norm(orig-Q)?

(2) The following simple example can not obtian right result (d=0.5, Q=[0 0 0]). (i get a wrong result: d =0.3333, Q = [0.5 -0.25 0]

==================

v = [-1 0 0 ; 1 1 0 ; 1 -1 0];

f = [1 2 3];

t = opcodemesh(v',f');

from = [0 0 0.5]';

to = [0 0 -1]';

[hit,d,trix,bary] = t.intersect(from,to-from);

Q = repmat(bary(1,:),3,1).*v(f(trix,1),:)' + ...

repmat(bary(2,:),3,1).*v(f(trix,2),:)' + ...

repmat(1-bary(1,:)-bary(2,:),3,1).*v(f(trix,3),:)'

===========

I do not know why?

Thanks!

Vipin VijayanPS. The getting the points part should actually look like

Q = repmat(B(1,:),3,1).*v(:,f(1,T)) + ...

repmat(B(2,:),3,1).*v(:,f(2,T)) + ...

repmat(1-B(2,:)-B(1,:),3,1).*v(:,f(3,T));

where B = bary(:,hit), T = trix(hit).