4

I'm trying to make a surf plot that looks like:

enter image description here

So far I have:

x = [-1:1/100:1];
y = [-1:1/100:1];

[X,Y] = meshgrid(x,y);

Triangle1 = -abs(X) + 1.5;
Triangle2 = -abs(Y) + 1.5;

Z = min(Triangle1, Triangle2);

surf(X,Y,Z);
shading flat
colormap winter;
hold on;

[X,Y,Z] = sphere();
Sphere = surf(X, Y, Z + 1.5 );% sphere with radius 1 centred at (0,0,1.5)
hold off;

This code produces a graph that looks like :

enter image description here

  • A pyramid with square base ([-1,1]x[-1,1]) and vertex at height c = 1.5 above the origin (0,0) is erected.

  • The top of the pyramid is hollowed out by removing the portion of it that falls within a sphere of radius r=1 centered at the vertex.

So I need to keep the part of the surface of the sphere that is inside the pyramid and delete the rest. Note that the y axis in each plot is different, that's why the second plot looks condensed a bit. Yes there is a pyramid going into the sphere which is hard to see from that angle.

I will use viewing angles of 70 (azimuth) and 35 (elevation). And make sure the axes are properly scaled (as shown). I will use the AXIS TIGHT option to get the proper dimensions after the removal of the appropriate surface of the sphere.

Ahmed Abdelazim
  • 133
  • 1
  • 7
  • That is really challenging to do in MATLAB, it is not the right tool for this kind of geometry composition. Not saying it is impossible, but there are no primitives to help you do it. You must write the code to do all the calculations to figure out the intersection points and tessellate them appropriately to get a surface. Not for the faint of heart. – Cyb3rFly3r Apr 30 '16 at 01:31
  • There may be an alternative way of creating the plot that you have shown, but I am not sure if the method would be suitable for you. You can use equations for flat surfaces/sphere for the determination of the values of explicit points in the form of `z=f(x,y)`. Then you can use `meshgrid` to generate the data for a `surf` plot that should produce the plot that you have shown. If the solution is suitable for you, I can produce the code - let me know. –  Apr 30 '16 at 08:00

2 Answers2

1

Here is my humble suggestion:

N = 400; % resolution
x = linspace(-1,1,N);
y = linspace(-1,1,N);
[X,Y] = meshgrid(x,y);
Triangle1 = -abs(X)+1.5 ;
Triangle2 = -abs(Y)+1.5 ;
Z = min(Triangle1, Triangle2);
Trig = alphaShape(X(:),Y(:),Z(:),2);
[Xs,Ys,Zs] = sphere(N-1);
Sphere = alphaShape(Xs(:),Ys(:),Zs(:)+2,2);
% get all the points from the pyramid that are within the sphere:
inSphere = inShape(Sphere,X(:),Y(:),Z(:));
Zt = Z;
Zt(inSphere) = nan; % remove the points in the sphere
surf(X,Y,Zt)
shading interp
view(70,35)
axis tight

I use alphaShape object to remove all unwanted points from the pyramid and then plot it without them:

enter image description here

I know, it's not perfect, as you don't see the bottom of the circle within the pyramid, but all my tries to achieve this have failed. My basic idea was plotting them together like this:

hold on;
Zc = Zs;
inTrig = inShape(Trig,Xs(:),Ys(:),Zs(:)+1.5);
Zc(~inTrig) = nan;
surf(Xs,Ys,Zc+1.5)
hold off

But the result is not so good, as you can't really see the circle within the pyramid.

Anyway, I post this here as it might give you a direction to work on.

EBH
  • 10,350
  • 3
  • 34
  • 59
1

An alternative to EBH's method.

A general algorithm from subtracting two shapes in 3d is difficult in MATLAB. If instead you remember that the equation for a sphere with radius r centered at (x0,y0,z0) is

r^2 = (x-x0)^2 + (y-y0)^2 + (z-z0)^2

Then solving for z gives z = z0 +/- sqrt(r^2-(x-x0)^2-(y-y0)^2) where using + in front of the square root gives the top of the sphere and - gives the bottom. In this case we are only interested in the bottom of the sphere. To get the final surface we simply take the minimum z between the pyramid and the half-sphere.

Note that the domain of the half-sphere is defined by the filled circle r^2-(x-x0)^2-(y-y0)^2 >= 0. We define any terms outside the domain as infinity so that they are ignored when the minimum is taken.

N = 400; % resolution
z0 = 1.5;  % sphere z offset
r = 1;     % sphere radius
x = linspace(-1,1,N);
y = linspace(-1,1,N);
[X,Y] = meshgrid(x,y);
% pyramid
Triangle1 = -abs(X)+1.5 ;
Triangle2 = -abs(Y)+1.5 ;
Pyramid = min(Triangle1, Triangle2);    
% half-sphere (hemisphere)
sqrt_term = r^2 - X.^2 - Y.^2;
HalfSphere = -sqrt(sqrt_term) + z0;
HalfSphere(sqrt_term < 0) = inf;
Z = min(HalfSphere, Pyramid);
surf(X,Y,Z)
shading interp
view(70,35)
axis tight

enter image description here

jodag
  • 19,885
  • 5
  • 47
  • 66
  • Nice! This is definitely the right direction for this specific problem, with simple geometric shapes. – EBH Jan 15 '18 at 20:19
  • @EBH I hope I didn't undermine a year of hard work on this question! – jodag Jan 15 '18 at 20:25
  • A year?! more like an hour... and I was aware of your idea, but really wanted to solve this with a more 'general' way, that can treat any cloud of points. – EBH Jan 15 '18 at 20:29
  • I was thinking the same. Was thinking of intersecting vertical lines with sphere at each missing point but that's not completely general either. – jodag Jan 15 '18 at 20:33