% Schematic25.m %
clf; clear;
nx = 100;ny = 300;
w = imread('bone.jpg');
w = imresize(w,[nx,ny]);
d=double(w);
phi=d(:,:,3);
a = zeros(size(phi));
a(phi<220) = 1;
for it = 1:30
    old_a = a;
    for i = 2:nx-1
        for j = 2:ny-1
            a(i,j) = old_a(i,j) + 0.125 *...
                ( a(i+1,j)+a(i,j+1)+a(i-1,j)+a(i,j-1)-4*a(i,j) );
        end
    end
    
end
mesh(a)
[c,h] = contour(a,[0.5, 0.5]);
x = c(1,2:end);y = c(2,2:end);
plot(x,y,'k-','linewidth',4);
hold on;
patch(x,y,[0.7, 0.7, 0.7]);
axis([0 ny 0 nx]);
axis off
line([0 ny ny 0 0],[0 0 nx nx 0],'color','k','linewidth',2);
h = text(0.5*ny,0.5*nx,'$\Omega$');
set(h,'interpreter','latex','Position',[142.5 49.0 0],'fontsize',35);
set(gcf,'PaperPositionMode','auto')