function [X,Y,Z,Z2,Z3] = RP014(option)
%function remotesense rev 0.13
% Bathymetry sensing with spiking neuron anomaly detection
% © 2021 Johannes Lawen, license: creative commons CC BY-NC-SA
% Email: jl@environment.report or joh.lawen@gmail.com
[X,Y,Z] = RS013(false);
Z2 = excisetri(X,Y,Z); %Excise land
if option == true
    [X,Y,Z3] = RS013(true);
    Z3 = UHI(Z3);
    Z3 = Z3.*ceil(Z2/100); %Excise filtered
end
xmi = 2500; xma = xmi + 450; ymi = 2900; yma = ymi + 600;
imagesc(X(xmi:xma),Y(yma:-1:ymi),-Z3(xmi:xma,ymi:yma)'); view(2);
axis([X(xmi) X(xma) Y(ymi) Y(yma)])
colormap(turbo)
gray = [90 90 90]/255;
set(gca,'Color',gray)
gray = [90 90 90]/255;
pbaspect([450 600 1])
box on
set(gcf, 'Renderer', 'Painters');
set(gca,'xtick',[]);
set(gca,'ytick',[])

end

function [X,Y,Z] = RS013(option)
% Converted to UTM with limited accuracy, e.g. rounding away 6126913
% Zone 39R (converted with https://www.latlong.net/lat-long-utm.html)
xmin = 379185.04;
ymin = 2757585.04;
xmax = 608114.95;
ymax = 2991015.04;
SNN = option; %false without SNN, true with SNN
nz = 1650; %Number of increments in IDed bathymetry depth
th = [6 6 6];
soundings %Load survey measurements into array named 'depth'
[imagename, reflistfile, depthpfile, synapsepfile] = fname;
for t = 1:1
    for n = 1:3 %For band 2 to 4
        i1 = Tiff(strcat(string(t),imagename(n,1:7)),'r');
        imageB = read(i1);
        %Referencing remote sensed with ground soundings
        X = linspace(xmin, xmax, size(imageB,2));
        Y = linspace(ymin, ymax, size(imageB,1));
        Yi = linspace(ymax, ymin, size(imageB,1)); %Inverted as per lcd coordinate
        m = 0;
        for i = 1:size(depth,1) %Remove dry soundings
            if depth(i,3) < 0
                m = m + 1;
                depthf(m,1:3) = depth(i,1:3);
            end
        end
        list(1:size(depthf,1),1:4) = NaN; % northern coordinate inverted
        for i = 1:size(depthf,1)
            rmin = xmax - xmin;
            for j = 2900:3500 %1:length(X)
                for k = 5300:5800 %1:length(Y) Screen coordinate inverted
                    r = sqrt( ( depthf(i,1) - X(j) )^2 ...
                        + ( depthf(i,2) - Yi(k) )^2 );
                    if r < rmin
                        rmin = r; jmin = j; kmin = k; vmin = imageB(k,j);
                    end %Satellite image pixel closest to sounding/min distance
                end
            end
            list(i,1) = jmin; list(i,2) = kmin; %Recording k for inverted coordina.
            list(i,3) = -depthf(i,3); list(i,4) = vmin;
        end
        %Average if soundings share pixel (include dist weighted averaging in rev2)
        %shared pixel list slist
         for i = 1:size(list,1)
             r = 1;
             j1 = list(i,1); k1 = list(i,2);
             for m = 1:size(list,1)
                 if j1 == list(m,1) && k1 == list(m,2)
                     slist(r) = m;
                     r = r + 1;
                 end
             end
             list(slist,3) = mean(list(slist,3));
             clear slist
         end
         %End of averaging soundings that share pixel
         testlist = list; %Store before band dependent mods for error ref.
         %Average if soundings share shading
        %shared pixel list slist
        for i = 1:size(list,1)
            r = 1;
            v1 = list(i,4);
            for m = 1:size(list,1)
                if v1 == list(m,4)
                    slist(r) = m;
                    r = r + 1;
                end
            end
            list(slist,3) = mean(list(slist,3));
            clear slist
        end
        %End of averaging soundings that share shading
        save(strcat(string(t),reflistfile(n,1:14)),'list')

        listA = list(1:3:end,1:4); %Partial list of soundings

        %save(strcat(string(t),refplistfile(n,1:15)),'list')
        %load('reflistB3.mat')
        minv = min(listA(:,4)); %Minimal reflectivity in soundings
        %CALIBRATION
        %Set up sparsely stored synapses for bands sb in perceptron as per survey
        sb(1:2^16) = 0;
        for i = 1:1:size(listA,1) %Go through all soundings
            for j = 0:nz-1 %Go through all depths with 1 cm increments
                if j <= 100*listA(i,3) && j+1 > 100*listA(i,3) %Interval
                    sb(listA(i,4)) = j; %Assign depth/output neuron to shade/input
                end
            end
        end
        %Interpolating for unsampled bandvalues with preceeding deeper depth
        pd = NaN;
        for i = 1:1:length(sb) %Go through all band value input neurons
            if sb(i) == 0 %If no (outgoing) synapse
                if pd > 0 %If a preceding deeper value pd has been recorded
                    sb(i) = pd; %Assign deeper depth
                end
            else
                pd = sb(i); %Set the preceding deeper depth
            end
        end
        %Remote Sensing
        Z(1:length(X),1:length(Yi)) = 0;
        for j = 1:length(X) %Go through each image pixel
            for k = 1:length(Yi) %in each image line
                Z(j,k) = sb(imageB(length(Yi)-k+1,j)+1)/100; %Associate
            end %with screen inverted coordinates
        end
        for j = 1:length(X) %Go through each image pixel
            for k = 1:length(Yi) %in each image line
                if Z(j,k) == 0
                    Z(j,k) = NaN;
                end
            end %with screen inverted coordinates
        end
        Y = Y - 2000000; %False northing
        xmi = 2500; xma = 2950; ymi = 2900; yma = 3500;
        %subplot(2,2,1)
        mesh(X(xmi:xma),Y(ymi:yma),Z(xmi:xma,ymi:yma)')
        axis([X(xmi) X(xma) Y(ymi) Y(yma)])
        view(2)
        pbaspect([1 1 1]) %Even axis scalling
        hAx = gca; %Handle to current axes
        hAx.XAxis.Exponent=0; %Don't use exponent in axis label
        hAx.YAxis.Exponent=0;
        xlabel('Eastern UTM/[$m$]','Interpreter','latex')
        ylabel('Excess Northern UTM $>2\:10^6$/[$m$]','Interpreter','latex')
        set(gca,'TickLabelInterpreter','latex');
        box on
        grid off
        set(gca,'layer','top')
        hold off
        m = 1;
        for i = 2:3:size(list,1)
            e(m) = abs( testlist(i,3) - ...
                Z(list(i,1),length(Yi)-list(i,2)+1) )/ ...
                testlist(i,3);
            m = m + 1;
        end
        mean(e)*100
        %2. SNN outlier filtration
        if SNN == true
        Z = SNNfilter(Z,th(n));
        %Z = fillmissing(Z,'linear');
        m = 1;
        for i = 2:3:size(list,1)
            if isnan(Z(list(i,1),length(Yi)-list(i,2)+1))
            else
                e(m) = abs( testlist(i,3) - ...
                    Z(list(i,1),length(Yi)-list(i,2)+1) )/ ...
                testlist(i,3);
                m = m + 1;
            end
        end
        mean(e)*100
        end
        save(strcat(string(t),depthpfile(n,1:13)),'Z')  %***** before or aft
        save(strcat(string(t),synapsepfile(n,1:14)),'sb') %***** loop
    end
    %3. Combine estimates from different bands
    clearvars -except t testlist Zp % Zi listi
    xmin = 379185.04;
    ymin = 2757585.04;
    xmax = 608114.95;
    ymax = 2991015.04;
    [imagename, reflistfile, depthpfile, synapsepfile] = fname; %Initialize bathymetry function file names
    i1 = Tiff(strcat(string(t),imagename(1,1:7)),'r');
    imageB2 = read(i1);
    i1 = Tiff(strcat(string(t),imagename(2,1:7)),'r');
    imageB3 = read(i1);
    i1 = Tiff(strcat(string(t),imagename(3,1:7)),'r');
    imageB4 = read(i1);
    X = linspace(xmin, xmax, size(imageB2,2));
    Y = linspace(ymin, ymax, size(imageB2,1));
    Yi = linspace(ymax, ymin, size(imageB2,1)); %Screen coordinate inversion
    index = [1 2 3]; ind2 = [NaN NaN]; w(index) = 0;
    load(strcat(string(t),reflistfile(1,1:14))) %Load list referencing soundings to image
    ones(1:size(list,1)) = 1;

    load(strcat(string(t),synapsepfile(1,1:14))) %Load synapse linking
    sb2 = sb;
    load(strcat(string(t),synapsepfile(2,1:14))) %Load synapse linking
    sb3 = sb;
    load(strcat(string(t),synapsepfile(3,1:14))) %Load synapse linking
    sb4 = sb;

    %LOAD ALSO SYNAPSE LINKING and all three images
    %Band 2
    load(strcat(string(t),depthpfile(1,1:13))) %Load depth extracted from band 2
    Z1(:,:) = Z(:,:);
    load(strcat(string(t),reflistfile(1,1:14))) %Load list referencing soundings to image
    list1(:,:) = list(:,:);
    %Band 3
    load(strcat(string(t),depthpfile(2,1:13))) %Load depth extracted from band 3
    Z2(:,:) = Z(:,:);
    load(strcat(string(t),reflistfile(2,1:14))) %Load list referencing soundings to image
    list2(:,:) = list(:,:);
    %Band 4
    load(strcat(string(t),depthpfile(3,1:13))) %Load depth extracted from band 4
    Z3(:,:) = Z(:,:);
    load(strcat(string(t),reflistfile(3,1:14))) %Load list referencing soundings to image
    list3(:,:) = list(:,:);
    clear Z
    ts2(1:2^16) = 0; ts3 = ts2; ts4 = ts2;
    list1B = list1(2:3:size(list,1),1:4);
    list2B = list2(2:3:size(list,1),1:4);
    list3B = list3(2:3:size(list,1),1:4);
    for i = 2:3:size(list,1) %1:size(list1B,1)
        clear d
        d(1) = Z1(list1(i,1),length(Yi)-list1(i,2)+1);
        d(2) = Z2(list1(i,1),length(Yi)-list1(i,2)+1);
        d(3) = Z3(list1(i,1),length(Yi)-list1(i,2)+1);
        d0 = testlist(i,3);
        [mind, imin] = min(d);
        [maxd, imax] = max(d);
        w(index) = 0;
        if isnan(sum(d(:)))
            w(:) = NaN;
        elseif mind >= d0 %If all bands sense deeper than sounding, take shallowest
            w(imin) = 1;
        elseif maxd <= d0 %& if all sense shallower than sounding, take deepest
            w(imax) = 1;
        elseif ( d0 - d(1) )*( d0 - d(2) )*( d0 - d(3) ) > 0 %If 2 bands deeper,
            ind2 = index(index~=imax); %excise deepest
            w(ind2(1)) = ( d0 - d(ind2(2)) )/( d(ind2(1)) - d(ind2(2)) );
            w(ind2(2)) = ( d0 - d(ind2(1)) )/( d(ind2(2)) - d(ind2(1)) );
        elseif ( d0 - d(1) )*( d0 - d(2) )*( d0 - d(3) ) < 0 %two shallower, so excise shallowest
            ind2 = index(index~=imin);
            w(ind2(1)) = ( d0 - d(ind2(2)) )/( d(ind2(1)) - d(ind2(2)) );
            w(ind2(2)) = ( d0 - d(ind2(1)) )/( d(ind2(2)) - d(ind2(1)) );
        else %if one or more band estimates is NaN
            ind2 = index(~isnan(index));
            if length(ind2) == 2
            w(ind2(1)) = ( d0 - d(ind2(2)) )/( d(ind2(1)) - d(ind2(2)) );
            w(ind2(2)) = ( d0 - d(ind2(1)) )/( d(ind2(2)) - d(ind2(1)) );
            elseif length(ind2) == 1
                w(ind2(1)) = 1;
            else
                w(2) = 1;
            end
        end
        %Assign synapse transmissivity/weighing factor, averaged for shade x
        ts2(list1(i,4)) = ts2(list1(i,4)) ...
            + w(1)/sum(ones( list1B(:,4) == list1(i,4) ));
        ts3(list2(i,4)) = ts3(list2(i,4)) ...
            + w(2)/sum(ones( list2B(:,4) == list2(i,4) ));
        ts4(list3(i,4)) = ts4(list3(i,4)) ...
            + w(3)/sum(ones( list3B(:,4) == list3(i,4) ));
    end
    %Interpolating for unsampled bandvalues with preceeding deeper depth
    pt1 = NaN; pt2 = NaN; pt3 = NaN;
    for i = 1:1:length(sb2) %Go through all transmissivities
        if ts2(i) == 0 %If no (outgoing) synapse
            if pt1 > 0 %If a preceding value pt has been recorded,
                ts2(i) = pt1; %then assign it
            end
        else
            pt1 = ts2(i); %Set the preceding deeper depth
        end
        if ts3(i) == 0 %Repeat for band 3 neurons
            if pt2 > 0
                ts3(i) = pt2;
            end
        else
            pt2 = ts3(i);
        end
        if ts4(i) == 0 %Repeat for band 4 neurons
            if pt3 > 0
                ts4(i) = pt3;
            end
        else
            pt3 = ts4(i);
        end
    end
    %Remote Sensing
    Z(1:length(X),1:length(Yi)) = 0;
    for j = 1:length(X) %Go through each image pixel
        for k = 1:length(Yi) %in each image line
            ki = length(Yi)-k+1;
            Z(j,k) = ( ts2(imageB2(ki,j)+1)*sb2(imageB2(ki,j)+1) ...
                +  ts3(imageB3(ki,j)+1)*sb3(imageB3(ki,j)+1) ...
                +  ts4(imageB4(ki,j)+1)*sb4(imageB4(ki,j)+1) )/ ...
                ( ts2(imageB2(ki,j)+1) + ts3(imageB3(ki,j)+1) ...
                + ts4(imageB4(ki,j)+1) );
            Z(j,k) = Z(j,k)/100;
        end %with screen inverted coordinates
    end
    save('tridepthf2.mat','Z')
    Y = Y - 2000000; %False northing
    xmi = 2500; xma = 2950; ymi = 2900; yma = 3500;
    mesh(X(xmi:xma),Y(ymi:yma),Z(xmi:xma,ymi:yma)')
    axis([X(xmi) X(xma) Y(ymi) Y(yma)])
    view(2)
    pbaspect([1 1 1]) %Even axis scalling
    hAx = gca; %Handle to current axes
    hAx.XAxis.Exponent=0; %Don't use exponent in axis label
    hAx.YAxis.Exponent=0;
    xlabel('Eastern UTM/[$m$]','Interpreter','latex')
    ylabel('Excess Northern UTM $>2\:10^6$/[$m$]','Interpreter','latex')
    set(gca,'TickLabelInterpreter','latex');
    box on
    grid off
    set(gca,'layer','top')
clear e
        m = 1;
        for i = 3:3:size(list,1)
            if ~isnan(Z(list(i,1),length(Yi)-list(i,2)+1))
            e(m) = abs( testlist(i,3) - ...
                Z(list(i,1),length(Yi)-list(i,2)+1) )/ ...
                testlist(i,3);
            m = m + 1;
            end
        end
    mean(e)*100
end

end

function Z = SNNfilter(Z,thr)

%Pixel neighbor molecule
dx(1:4,1:8) = 0; dy(1:3,1:8) = 0;
%r = 1
%(i-1,j)
%(i+1,j)
%(i,j-1)
%(i,j+1)
dx(1,1:4) = [-1 1 0 0]; dy(1,1:4) = [0 0 -1 1];
%r = sqrt(2)
%(i-1,j-1)
%(i-1,j+1)
%(i+1,j-1)
%(i+1,j+1)
dx(2,1:4) = [-1 -1 1 1]; dy(2,1:4) = [-1 1 -1 1];
%r = 2
%(i-2,j)
%(i+2,j)
%(i,j-2)
%(i,j+2)
dx(3,1:4) = [-2 2 0 0]; dy(3,1:4) = [0 0 -2 2];
%r = sqrt(5)
%(i-2,j-1)
%(i+2,j-1)
%(i-1,j-2)
%(i-1,j+2)
%(i-2,j+1)
%(i+2,j+1)
%(i+1,j-2)
%(i+1,j+2)
dx(4,:) = [-2 2 -1 -1 -2 2 1 1]; dy(4,:) = [-1 -1 -2 2 1 1 -2 2];
ce = [4 4 4 8]; %Circle elements
f = Z; %Initilize
clear out
out = Z;
f(:,:) = 0; out(:,:) = 0;
r(1:5) = [1 sqrt(2)-1 2-sqrt(2) sqrt(5)-2 0];
for x = 3:size(out,1)-2
    for y = 3:size(out,2)-2
        for i = 1:4
            s = 0;
            for j = 1:ce(i)
                s = s + abs( Z(x,y) - Z(x+dx(i,j),y+dy(i,j)) )/ ...
                    abs( Z(x,y) )/sum(r(1:i));
%                s = s + 2*abs( Z(x,y) - Z(x+dx(i,j),y+dy(i,j)) )/ ...
%                    abs( Z(x,y) + Z(x+dx(i,j),y+dy(i,j)) )/sum(r(1:i));
            end
            f(x,y) = f(x,y) + s;
            if f(x,y) > thr %6 %10
                out(x,y) = 1;
            %    Z(x,y) = NaN;
            end
            f(x,y) = f(x,y)*exp(-r(i+1)); %Decay
        end
    end
end
for x = 3:size(out,1)-2
    for y = 3:size(out,2)-2
        if out(x,y) == 1
            Z(x,y) = NaN;
        end
    end
end

end

function [imagename,reflistfile,depthpfile,synapsepfile]=fname
    imagename(1,1:7) = '/B2.TIF';
    imagename(2,1:7) = '/B3.TIF';
    imagename(3,1:7) = '/B4.TIF';
    reflistfile(1,1:14) = '/reflistB2.mat';
    reflistfile(2,1:14) = '/reflistB3.mat';
    reflistfile(3,1:14) = '/reflistB4.mat';
    depthpfile(1,1:13) = '/b2pdepth.mat';
    depthpfile(2,1:13) = '/b3pdepth.mat';
    depthpfile(3,1:13) = '/b4pdepth.mat';
    synapsepfile(1,1:14) = '/synapsep2.mat';
    synapsepfile(2,1:14) = '/synapsep3.mat';
    synapsepfile(3,1:14) = '/synapsep4.mat';
end

function Ze = excisetri(X,Y,Z)
Ze = Z + 100; % mark unprocessed pixels
xmi = 2500; xma = xmi + 450; ymi = 2900; yma = ymi + 600;
xst = xmi;
yst = ymi;
Ze(xst,yst) = Ze(xst,yst) - 100;
er(1) = 0.26; er(2) = 1.0; er(3) = 2.1;
no =1;
runs = 1;
while no > 0
	no = 0;
	for x = xmi:xma
		for y = ymi:yma
if Ze(x-1,y) < 100 | Ze(x+1,y) < 100 | Ze(x,y-1) < 100 | Ze(x,y+1) < 100 % if neighbor < 100, then processed & not NaN/dry
                if Ze(x,y) >= 100 %process if not yet processed
                    if y > (ymi+yma)/2  & x > xmi+(xma-xmi)/18 %high threshold for artificial islands
                        err = er(2);
                        if x > xmi+(xma-xmi)/2.5
                        if y < ymi+(yma-ymi)/1.23 & x < xmi+(xma-xmi)/1.23
                           err = er(3);
                        end
                        end
                    else
                        err = er(1);
                    end
					if Ze(x,y) < 100 + err %mark as dry
						Ze(x,y) = NaN;
					else %or wet, mark as processed & wet
                        Ze(x,y) = Ze(x,y) - 100;
						no = no + 1;
                        runs = runs + 1;
					end
				end
			end
		end
	end
end
Ze(Ze>90) = NaN;
mesh(X(xmi:xma),Y(ymi:yma),-Ze(xmi:xma,ymi:yma)'); view(2);
colormap(turbo)
sand = [194, 178, 128]/255;
gray = [90 90 90]/255;
set(gca,'Color',gray)
axis([X(xmi) X(xma) Y(ymi) Y(yma)])
pbaspect([1 1 1])
hAx = gca; %Handle to current axes
hAx.XAxis.Exponent=0; %Don't use exponent in axis label
hAx.YAxis.Exponent=0;
xlabel('Eastern UTM/[$m$]','Interpreter','latex')
ylabel('Excess Northern UTM $>2\:10^6$/[$m$]','Interpreter','latex')
title('Bathymetry/[$m$]','Interpreter','latex')
cb = colorbar;
box on
set(gca,'TickLabelInterpreter','latex')
set(cb,'TickLabelInterpreter','latex')
end

%Export as follows for Latex as .eps
% ascertain file size below 300 MB in case of overleaf.com use
% Ascertain vector rendering of the .eps with:
% exportgraphics(gcf,'filename.eps','ContentType','vector')
% Zip and ascertain zipped folder below 50 MB instead of overleaf.com
% Upload folder when launching new project in overleaf as it will be
% automatically decompressed.
% Add .tex files to the new project folder that holds the .eps

%Export as follows for MS Word as .emf
% Save file as Matlan figure under filename.fig
% Open with fig = openfig('sensefig.fig')
% Force to use 'painters' instead of 'openGL' as renderer with
% fig.Renderer = 'painters'
% Store as vector emf
% saveas(gcf,'filename.emf')
% Insert .emf as image into word document
% Render word into pdf with print function

%Test with e.g. UHI([1 0 1; 0 0 0; 2 2 6])
function F = UHI(IC)                            %Unstruct. Heat Interp.
    F = ( fillmissing(IC,'linear',1) + ...      %Save iterations w/ initial
        fillmissing(IC,'linear',2) )/2;         %approx. interp. (1D linear
    IC(isnan(IC)) = 0;                          %Initial condition IC
    q = .1^4*max(IC,[],'all');                  %Tolerance q
	DL = 0.^IC;	ICL = 0.^DL;                    %Boolean domain- & IC loci
	while max(DL.*(HC(F)-F),[],'a') > q         %Fill domain/cells valued 0
		F = HC(F);                              %Central diff. heat smooth
		F = DL.*F + ICL.*IC;                    %Applying IC & BC @ IC loci
	end
end
function F = HC(F)                              %Centr. Diff. Heat Conduct.
	n = size(F,1); m = size(F,2);               %Last indices of domain
	NR(1:n) = (1:n) + 1; NL(1:n) = (1:n) - 1;	%Right & left neighbors
	NU(1:m) = (1:m) + 1; ND(1:m) = (1:m) - 1;	%Up & down neighbors
	NL(1) = 1; NR(n) = n; ND(1) = 1; NU(m) = m; %Modded Boundary indices
	dfx = F(NR,:) - 2*F(:,:) + F(NL,:);         %Central differences vs x
	dfy = F(:,NU) - 2*F(:,:) + F(:,ND);         %Central differences vs y
	F = ( dfx + dfy )*0.1 +  F;                 %Solve FDE of PDE
end

Published with MATLAB® R2021a