% Schematic26.m %
clear; clc; clf;
nx = 100; ny = 300;
w = imread('bone.jpg');
w = imresize(w,[nx,ny]);
x1 = linspace(0,300,ny);
y1 = linspace(0,100,nx);
[xx1,yy1] = meshgrid(x1,y1);
d = double(w);phi = d(:,:,3);
a = zeros(size(phi));
a(phi<220) = 1;
for it = 1:50
    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);
mesh(a*0);
view(2);
plot3(x,y,x*0+2,'k-','linewidth',1);
hold on;
ny = 10;nx = ny*3;h = 100/ny;
xt = linspace(0.5*h,300-0.5*h,nx);
yt = linspace(0.5*h,100-0.5*h,ny);
[xx,yy]=meshgrid(xt,yt);
b = interp2(xx1,yy1,a,xx,yy);
c = zeros(size(b));
c (b >= 0.50) = 1;
for j = 1:ny+1;
    line([0 300],[(j-1)*h (j-1)*h],'color','k');
end
for i = 1:nx+1;
    line([(i-1)*h (i-1)*h],[0 100],'color','k');
end
c = c';
for i=1:nx
    for j=1:ny
        if c(i,j) > 0.5
            patch([xt(i)-0.5*h xt(i)+0.5*h xt(i)+0.5*h xt(i)-0.5*h],[yt(j)-0.5*h yt(j)-0.5*h yt(j)+0.5*h yt(j)+0.5*h],[0.7 0.7 0.7]);
        end
    end
end
for i=1:nx
    for j=1:ny
        if c(i,j) < 0.5
            plot(xt(i)+0.5*h,yt(j),'kx','markersize',5);
            plot(xt(i),yt(j)+0.5*h,'kx','markersize',5);
            plot(xt(i)-0.5*h,yt(j),'kx','markersize',5);
            plot(xt(i),yt(j)-0.5*h,'kx','markersize',5);
        end
        if c(i,j) > 0.5
            if c(i+1,j) > 0.5
                plot(xt(i)+0.5*h,yt(j),'k.','markersize',8);
            end
            if c(i-1,j) > 0.5
                plot(xt(i)-0.5*h,yt(j),'k.','markersize',8);
            end
            if c(i,j+1) > 0.5
                plot(xt(i),yt(j)+0.5*h,'k.','markersize',8);
            end
            if c(i,j-1) > 0.5
                plot(xt(i),yt(j)-0.5*h,'k.','markersize',8);
            end
        end
    end
end
for i=1:nx
    for j=1:ny
        plot(xt(i),yt(j),'ko','markersize',2);
    end
end
axis image;axis off;view(2);
axis([0 300 0 100]);
h = text(0.5*ny,0.5*nx,'$\Omega_{in}^h$');
set(h,'interpreter','latex','Position',[142.5 47.8],'fontsize',20);
h = text(0.5*ny,0.5*nx,'$\Omega_{out}^h$');
set(h,'interpreter','latex','Position',[2.0 90.8],'fontsize',20);
h = text(0.5*ny,0.5*nx,'$\Gamma$');
set(h,'interpreter','latex','Position',[280.0 80.8],'fontsize',20);
h = text(0.5*ny,0.5*nx,'$\Gamma^h$');
set(h,'interpreter','latex','Position',[285.0 32.8],'fontsize',20);
annotation('arrow',[0.857 0.84],[0.622 0.61],'Headlength',4);
annotation('arrow',[0.868 0.855],[0.4556 0.4556],'Headlength',4);
set(gcf,'PaperPositionMode','auto')