5.0

5.0 | 4 ratings Rate this file 28 Downloads (last 30 days) File Size: 257 KB File ID: #41504

Ray casting for deformable triangular 3D meshes

by

 

26 Apr 2013 (Updated )

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

| Watch this File

File Information
Description

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.

Acknowledgements

Example Matlab Class Wrapper For A C++ Class inspired this file.

MATLAB release MATLAB 8.1 (R2013a)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (14)
14 Feb 2014 Eric

Great submission. Very useful.

Thank you Germano Gomes for mex64 advice.

08 Oct 2013 Vipin Vijayan

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

04 Oct 2013 Germano Gomes

Note(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

03 Oct 2013 Germano Gomes

Very 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

03 Oct 2013 Germano Gomes  
29 May 2013 Vipin Vijayan

I 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.

20 May 2013 omar

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

18 May 2013 K G

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

13 May 2013 Vipin Vijayan

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

10 May 2013 omar

Hi...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

09 May 2013 Vipin Vijayan

Made 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));

09 May 2013 Vipin Vijayan

You 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));

09 May 2013 Zhou

I 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!

26 Apr 2013 Vipin Vijayan

PS. 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).

Updates
18 Jun 2014

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.

19 Jun 2014

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.

Contact us