Search code examples
matlabvectorizationtriangulation

Retaining color when subsampling a triangulated surface: get indices from reducepatch?

I have a very densely tessellated surface which looks like this: White matter dense

This surface is too densely tessellated for me, so I subsample it to get a coarser surface. To do this, I used Matlab's reducepatch function. This works pretty well: White matter subsampled

Unfortunately, the coloring is based on a variable called sulcal_depth, which is defined for every vertex of my tessellated surface. So I need to retain sulcal depth information only from the vertices which remain after subsampling. Essentially, I need reducepatch to give me not just the subsampled version of the surface, but also the indices of vertex points that it retained. If I know the preserved indices, I can just index my sulcal_depth variable to get the new depth map.

Currently, I'm doing this as follows (this is also how I colored the subsampled version above):

function indices = compute_reduced_indices(before, after)
%% Function to compute the indices of vertices preserved during an operation of
%  reducepatch. This allows you to use reducepatch to subsample a surface and
%  re-compute an original signal on the vertices for the new subsampled mesh

indices = zeros(length(after), 1);
for i = 1:length(after)
    dotprods = (before * after(i, :)') ./ sqrt(sum(before.^2, 2));
    [~, indices(i)] = max(dotprods);
end

But as you might imagine, this is pretty slow, because of the for loop over vertices. I don't have enough memory to vectorize the loop and compute the full dot product matrix in one go.

Is there a smart way to get reducepatch to give me indices, or an alternative approach (with or without reducepatch) that's faster?

Solution

  • If reducepath only delete some vertex but doesn't change the coordinate of the preserved points you can use the function ismember:

    %Load the "flow" matlab's dataset.
[x,y,z,v] = flow(100);
%Patch the isosurface
p = patch(isosurface(x,y,z,v,-3));
%Reducepatch
rp = reducepatch(p,0.15);
%Create an index of the preserved vertex.
[ind,loc] = ismember(p.Vertices,rp.vertices,'rows'); 
%Checksum
sum(find(ind) == sort(indices)) == length(indices) %should be = 1

%if you want to preserve the index order:
locb = loc(ind);
subind = find(ind);
[~,revsor] = sort(locb);
ind = subind(revsor);

    BENCHMARKING

    [x,y,z,v] = flow(100);

p = patch(isosurface(x,y,z,v,-3));
rp = reducepatch(p,0.15);
tic
ind = ismember(p.Vertices,rp.vertices,'rows');
toc

before = p.Vertices;
after = rp.vertices;

tic 
indices = zeros(length(after), 1);
for i = 1:length(after)
    dotprods = (before * after(i, :)') ./ sqrt(sum(before.^2, 2));
    [~, indices(i)] = max(dotprods);
end
toc

    RESULT

    Elapsed time is 0.196078 seconds. %ismember solution 
Elapsed time is 11.280293 seconds. %dotproduct solution