%Size of the synthetic image (Our sampled function)
n=100;
%The synthetic takes some values 
f=ones(n,n)*1e100;
f(20,20:70)=0;
f(20:70,20:40)=0;
f(70,70)=0;
%First dimension, for all columns
for c=1:n
    k=1;
    v=1;
    z(1)=-1e100;
    z(2)=1e100;
    for q=2:n
        s=((f(q,c)+q^2)-(f(v(k),c)+v(k)^2))/(2*q-2*v(k));
        while (s<z(k))
            k=k-1;
            s=((f(q,c)+q^2)-(f(v(k),c)+v(k)^2))/(2*q-2*v(k));
        end
        k=k+1;
        v(k)=q;
        z(k)=s;
        z(k+1)=1e100;
    end
    k=0;
    for q=1:n
        while (z(k+1)<q)
            k=k+1;
        end
        g(q,c)=(q-v(k))^2+f(v(k),c);
    end
end

%Second dimension, for all rows
for c=1:n
    k=1;
    v=1;
    z(1)=-1e100;
    z(2)=1e100;
    for q=2:n
        s=((g(c,q)+q^2)-(g(c,v(k))+v(k)^2))/(2*q-2*v(k));
        while (s<z(k))
            k=k-1;
            s=((g(c,q)+q^2)-(g(c,v(k))+v(k)^2))/(2*q-2*v(k));
        end
        k=k+1;
        v(k)=q;
        z(k)=s;
        z(k+1)=1e100;
    end
    k=0;
    for q=1:n
        while (z(k+1)<q)
            k=k+1;
        end
        y(q,c)=(q-v(k))^2+g(c,v(k));
    end
end

subplot(1,2,1)
imagesc(f);
subplot(1,2,2)
%Compute the squared root of the output (not necessary in most
%applications, done here for better visualization
imagesc(y.^.5);
axis tight