Create new isosurface from a isosurface
5 views (last 30 days)
Show older comments
Hi to everybody.
I have an isosurface, generated with isosurface function, so I have a struct in my workspace as follow
surface
tri: [93372x3 int32]
vert: [47915x3 single]
when I plot it I have the results showed in attacched image.
Now, I want to select a subset of points of isosurface and generate a new isosurface starting from this selected points.
For example, I want to extract from the isosurface only the points included inside the red circle, selecting with datacursor the center of each black dot and creating a new isosurface starting from the coordinates of each black dot.
I hope it's clear and that you can help me.
Thanks in advance
Best regards
Luigi
0 Comments
Answers (7)
Mike Garrity
on 27 Mar 2015
The simplest way is the following.
Let's start with a really simple example:
[x y z v] = flow;
q = z./x.*y.^3;
p = patch(isosurface(x, y, z, q, -.08, v));
p.FaceColor = 'yellow';
view(3)
I want to keep everything between X=2 and X=4.
minx = 2;
maxx = 4;
inverts = p.Vertices(:,1) >= minx & p.Vertices(:,1) <= maxx;
The variable inverts is now a mask that is the size of my Vertices array and has a true for each vertex that I want to keep. Now I need to turn that into a mask of the faces I want to keep. To do that, I need to look at each column of the Faces array:
infaces = inverts(p.Faces(:,1)) & inverts(p.Faces(:,2)) & inverts(p.Faces(:,3));
Now infaces is a mask that is the size of my Faces array and has a true for every face where all three vertices have trues in inverts. Now we can use that to shrink the Faces array.
p.Faces = p.Faces(infaces,:);
The one problem with this approach is that the patch's Vertices array still contains all of those vertices that have been excluded. That's generally OK, but you might want to get rid of them. The problem is that to do that, you need to translate all of the indices in the Faces array from indices in the old, big Vertices array into indices into the new, small Vertices array. That's a messy bit of bookkeeping.
0 Comments
Luigi
on 31 Mar 2015
2 Comments
Mike Garrity
on 31 Mar 2015
I warned you!
It's actually not as bad as it looks at first. There are a couple of ways to do it. One of the simplest involves find:
%%Throw out unused verts & faces
new_verts = p.Vertices(inverts,:);
new_faces = p.Faces(infaces,:);
%%Fix faces
n2o = find(inverts);
for i=1:numel(new_faces)
new_faces(i) = find(n2o==new_faces(i));
end
I think that this can be improved in terms performance and robustness, but it doesn't need to get a lot uglier than this.
Luigi
on 2 Apr 2015
1 Comment
Mike Garrity
on 2 Apr 2015
I'm not sure. If I assemble all the bits from this thread, I get this:
[x y z v] = flow;
q = z./x.*y.^3;
p = patch(isosurface(x, y, z, q, -.08, v));
p.FaceColor = 'yellow';
view(3)
minx = 2;
maxx = 8;
miny = -3;
maxy = 3;
minz = 0;
maxz = 3;
inverts = (p.Vertices(:,1) >= minx & p.Vertices(:,1) <= maxx) ...
& (p.Vertices(:,2) >= miny & p.Vertices(:,2) <= maxy) ...
& (p.Vertices(:,3) >= minz & p.Vertices(:,3) <= maxz);
infaces = inverts(p.Faces(:,1)) & inverts(p.Faces(:,2)) & inverts(p.Faces(:,3));
new_verts = p.Vertices(inverts,:);
new_faces = p.Faces(infaces,:);
n2o = find(inverts);
for i=1:numel(new_faces)
new_faces(i) = find(n2o==new_faces(i));
end
patch('Faces',new_faces,'Vertices',new_verts,'FaceColor','green')
If I run it, I get the following picture:
Is that what you get? How does it differ from what you expected?
Luigi
on 3 Apr 2015
1 Comment
Mike Garrity
on 3 Apr 2015
Sorry, the error you hit running my code is because I'm using a syntax that was introduced in R2014b. I think you're running an older version.
In R2014b, the following are equivalent:
set(p,'FaceColor','yellow')
p.FaceColor = 'yellow'
But only the first worked in earlier versions. It's pretty simple to convert between the two forms, but that line that sets inverts would need a temp variable:
v = get(p,'Vertices');
inverts = (v(:,1) >= minx & v(:,1) <= maxx) ...
& (v(:,2) >= miny & v(:,2) <= maxy) ...
& (v(:,3) >= minz & v(:,3) <= maxz);
But I don't think that's the problem you've hit. At the point where you did this:
surface.vert = surface.vert(inverts,:);
surface.tri = surface.tri(infaces,:);
what do the variables inverts & infaces look like? In my example, they look like this:
>> any(inverts)
ans =
1
>> length(inverts)
ans =
2318
>> sum(inverts)
ans =
688
>> any(infaces)
ans =
1
>> length(infaces)
ans =
4344
>> sum(infaces)
ans =
1271
That's saying that I've selected 688 out of 2318 vertices and 1271 out of 4344 triangles.
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!