%
% Function minsurf
%
% If Z is a matrix of heights of a surface over a meshgrid [X,Y]
% minsurf uses the function tarea.m and the MATLAB function
% fmins to find the new matrix W = minsurf(X,Y,Z) that agrees
% with Z on the boundary, but which minimizes the surface area
% as a function of the interior values of Z. W = minsurf(X,Y,Z)
% is the new matrix. The new matrix W may not provide the
% minimal surface. You may need to iterate minsurf several
% times, taking W1 = minsurf(X,Y,W), W2 = minsurf(X,Y,W1), etc.
% The input surface and the modified surface are graphed
% together using tsurf for comparison.
% Depending on the machine, this program can run quite
% slowly. It is best not to use more than 40 interior points.
%
% Example: [X,Y] = meshgrid(0:.2:1)
% Z = 0*X + 1;
% Z(2:4, 2:4) = .5 + zeros(3,3);
% is a surface that is equal to 1 on the boundary of the
% unit square, and 0 at the interior points. The minimal
% surface matrix is obviously W = 0*X +1.
%
function Zout = minsurf(X,Y,Z)
n = size(X,2)-1; m = size(X,1)-1;
nn = n-1;
dx = X(1,2) - X(1,1); dy = Y(2,1) - Y(1,1);
z0 = zeros(1,(n-1)*(m-1));
for i = 2:m
z0((i-2)*nn +1: (i-1)*nn) = Z(i, 2:n);
end
initial_area = tarea(z0,dx,dy,Z)
% calculate the interior points which minimize the surface area.
zmin = fmins('tarea', z0,[],[],dx,dy,Z);
Zout = Z;
% Replace the interior values of Z with the values of zmin.
for i = 2:m
Zout(i,2:n) = zmin((i-2)*nn +1: (i-1)*nn);
end
final_area = tarea(zmin,dx,dy,Zout)
h1 = min(min(Z)); h2 = max(max(Z));
xmax = max(max(X)); xmin = min(min(X));
ymax = max(max(Y)); ymin = min(min(Y));
subplot(1,2,1)
tsurf(X,Y,Z)
title('original surface ')
subplot(1,2,2)
tsurf(X,Y,Zout)
axis([xmin, xmax, ymin, ymax, h1,h2])
title('minimal surface')