Efficient data storage for real-time simulation

I am currently working on a real-time simulation using MATLAB and MEX-files. I have a triangulated mesh of ~2000 vertices and ~4000 faces stored in two matrices v (Nx3, N:number of vertices and 3-dimensions) and f. I am using a mass-spring model to simulate deformations, for which I am constantly pulling and inserting data from the vertex matrix v to calculate displacements, distances, spring-forces etc. The reason why I am storing in Nx3 format is to be able to plot using trisurf(f,v(:,1), v(:,2), v(:,3)).

I am currently not reaching my target updates of 20 iterations/sec (currently around ~11 iterations/sec) and the profiler does not really target any specific area in the code for improvement. I am wondering if maybe a different approach to store the vertices may be more efficient? E.g. storing the vertices in a 3*Nx1 long vector instead. This would not make it suitable for trisurf though. Thoughts or comments?

Current version, pseudo-code:

for t=1:T

     for i = 1:N

        %Pull vertex i
        xi = v(i,:);  %xi = [x_i y_i z_i]

        for k = 1:numConn(i)  %number of connections numConn(i) is precomputed

           %Pull a vertex 
           j=connection(k,i);   %j equals a row in v, i.e a vertex
           xj=v(j,:);

           %Distance
           xij=xi-xj;
           d = sqrt(xij*xij');

        end
           ...
           %Update 
           v(i,:)=v(i,)+f(d,xi,vj,t);
     end

     %Plot
     trisurf(f,v(:,1),v(:,2),v(:,3));
     drawnow;
  end

[Cedric Wannaz]

You are mentioning in a comment below that it's not possible to vectorize inner loops completely; the only reason that I can see (but I spent 2 minutes on it, so I am probably wrong), it the way you are using an array of connections that is extracted from somewhere as a function of the vertex ID. Have you tried to use a matrix of connections instead and get v at t+1 as a product of matrices (force, connections_flag, interVertex_dist, v@t) ..?

[Linus]

The connections are not pulled from a function but from a matrix, say B. The rows of column i are indices (rows) in v connected to vertex i. Thus B(:,10) will return all vertices connected to vertex 10. I have tried vectorizing the innermost loop, but the result is slower.

[Wouter]

If you have the parallel processing toolbox, you could perhaps try to make the code run parallel on all you processors; this could theoretically speed up the execution (if processing power is actually the limiting factor).

If you are lucky, you could change the inner for loop into this:

        parfor k = 1:numConn(i)

however I suspect that your iterations depend on previous iterations (line v(i,:) = v(i,:)+f(d,xi,vj,t). In that case ignore my comment :)

[Linus]

It is true that the 1:numConn(i) loop contains the most expensive calculations, however, the loop is short as the numConn(i) is ususally less than 20 (meaning that each vertex is connected to roughly 20 other vertices) and I am not sure if threading this part may actually produce a speedup. As of now I am wondering if I should let MATLAB do all the pre-processing and drawing, leaving the actual mass-spring calculations and post-processing to MEX-files. I think that rearranging v into a vector instead of a matrix would improve cache-performance as lookups become more local.

Answer 1

[Wouter]

You could try to update these lines:

     trisurf(f,v(:,1),v(:,2),v(:,3));
     drawnow;

Into:

     if ~exist('htri','var')
         htri = trisurf(f,v(:,1),v(:,2),v(:,3));
     else
         set(htri,'vertices',v); % use if faces (f) do not change
         % set(htri,'vertices',v,'faces',f); % use if both vertices
         %                                   % and faces change
     end
     drawnow

This might be faster as you then do not run the entire trisurf function again.

[Wouter]

of course this only works when the drawing is the limiting factor :) and not the calculation.

[Linus]

Thank you for your suggestion, much appreciated. Unfortunately the drawing is not an limiting factor. For 500 time-steps your solution is a couple of seconds faster,

htri=trisurf(f,v(:,1),v(:,2),v(:,3));
for t=1:T
   mass_spring_function(v,...);
   set(htri,'vertices',v);
   drawnow;
end

[Wouter]

what is j exactly? all vertices connected to v(i,:)?

[Linus]

Yes, but in the above example it's just a single vertex to illustrate how I pull data from v. In the actual implementation I loop through all vertices sharing an edge with vertex i and calculate the distance (to determine the resulting force from vertex j to i).

[Wouter]

If the connections between the vertices remain the same you could try to find the connected vertices before the for-loops (e.g. using the ismember function) and then use these in the for loop in stead of calculating them each iteration.

Another option would be to vectorize you innermost for-loop (for i = 1:N). If this can be done depends however on the fact that subsequent iterations do not depend on each other.

[Linus]

You are correct, the connections stay the same and are pre-computed. It is not possible to vectorize the inner-loop completely (and vectorizing smaller parts actually makes it slower, probably because of JIT). The connections between vertices are not "well-structured" i.e. there does not exist a "connection pattern" between vertices. I have updated the pseudo-code.

Answer 2

[Linus]

The solution I went with was rewriting the most time-consuming parts into MEX-files. This gave me a 2x increase which put the bottleneck on the drawnow command (this matter is thus solved). I will accept Wouter's answer since it was indeed helpful.