simulate_points_dem_6

% simulate_points_dem_6: generate test image of Uwe 1994

 [points,BB,delta_x,sigma_k,sigma_s,sigma,out_in,dem] = ...
        simulate_points_dem_6(N,g_min,delta_g,sigma,outlier_percentage)

 Input:
    N:                   int, mesh size in x- and y-direction
    g_min:               double, min high level
    delta_g:             double, hight steps
    sigma:               double, std noise
    outlier_percentage:  double \in [0,1], fraction of outlies

 Output
    points:    double (N*N)x4, coordinates of dem grid points, [i, j, height(i,j), std]
    BB:        int 1x4, bounding box of dem
    delta_x:   double, gridsize/spacing
    sigma_k:   double, std curvature   
    sigma:     double, std noise
    out_in:    boolean (N*N)x1, outlier mask
    dem:       double NxN, DEM

 Wolfgang Förstner 07/14
 last changes: Susanne Wenzel 09/16, wf 4/2018
 wfoerstn@uni-bonn.de, wenzel@igg.uni-bonn.de

0001 %% simulate_points_dem_6: generate test image of Uwe 1994
0002 %
0003 % [points,BB,delta_x,sigma_k,sigma_s,sigma,out_in,dem] = ...
0004 %        simulate_points_dem_6(N,g_min,delta_g,sigma,outlier_percentage)
0005 %
0006 % Input:
0007 %    N:                   int, mesh size in x- and y-direction
0008 %    g_min:               double, min high level
0009 %    delta_g:             double, hight steps
0010 %    sigma:               double, std noise
0011 %    outlier_percentage:  double \in [0,1], fraction of outlies
0012 %
0013 % Output
0014 %    points:    double (N*N)x4, coordinates of dem grid points, [i, j, height(i,j), std]
0015 %    BB:        int 1x4, bounding box of dem
0016 %    delta_x:   double, gridsize/spacing
0017 %    sigma_k:   double, std curvature
0018 %    sigma:     double, std noise
0019 %    out_in:    boolean (N*N)x1, outlier mask
0020 %    dem:       double NxN, DEM
0021 %
0022 % Wolfgang Förstner 07/14
0023 % last changes: Susanne Wenzel 09/16, wf 4/2018
0024 % wfoerstn@uni-bonn.de, wenzel@igg.uni-bonn.de
0025 %
0026 function [points,BB,delta_x,sigma_k,out_in,dem] = ...
0027         simulate_points_dem_6(N,g_min,delta_g,sigma,outlier_percentage)
0028 
0029 %% plot settings
0030 ss = plot_init;
0031 
0032 BB=[1,1,N,N];
0033 
0034 %% set sigma of curvatures
0035 sigma_k = 1.0*sigma;
0036 
0037 %% start simulation
0038 A=ones(N,N)*g_min;
0039 
0040 %% squares: center (16,16)
0041 cent_x=round(2/7*N);
0042 cent_y=round(2/7*N);
0043 % lowest: half width = 12, z = 40
0044 h = round(12/70*N);
0045 z = g_min+delta_g;
0046 for i=cent_x-h:cent_x+h
0047     for j=cent_y-h:cent_y+h
0048         A(i,j)=z;
0049     end
0050 end
0051 % lowest: half width = 8, z = 60
0052 h = round(8/70*N);
0053 z = g_min+2*g_min;
0054 for i=cent_x-h:cent_x+h
0055     for j=cent_y-h:cent_y+h
0056         A(i,j)=z;
0057     end
0058 end
0059 % lowest: half width = 4, z = 80
0060 h = round(4/70*N);
0061 z = g_min+3*delta_g;
0062 for i=cent_x-h:cent_x+h
0063     for j=cent_y-h:cent_y+h
0064         A(i,j)=z;
0065     end
0066 end
0067 
0068 %% cylinder, sphere, block, peak
0069 cent_x = round(2/7*N);
0070 cent_y = round(5/7*N);
0071 r      = round(13/70*N);
0072 
0073 %% sphere
0074 for i = cent_x-r:cent_x
0075     for j = cent_y-r:cent_y
0076       dx = (i-cent_x)/r;
0077       dy = (j-cent_y)/r;
0078       dr = sqrt(dx^2+dy^2);
0079       if dr < 1;
0080           A(i,j) = g_min + 3*delta_g*(1-sqrt(dx^2+dy^2));
0081       end
0082     end
0083 end
0084 
0085 %% box with peak
0086 for i = cent_x-r:cent_x
0087     for j = cent_y:cent_y+r
0088       if dr < 1;
0089           A(i,j)=g_min + 3*delta_g;
0090       end
0091     end
0092 end
0093 A(cent_x-round(r/2),cent_y+round(r/2)) = g_min+4*delta_g;
0094 
0095 %% cylinder
0096 for i = cent_x:cent_x+r
0097     for j = cent_y-r:cent_y+r
0098       dx = (i-cent_x)/r;
0099       if dr < 1;
0100           A(i,j) = g_min + 3*delta_g*sqrt(1-dx^2);
0101       end
0102     end
0103 end
0104 
0105 %% diagonal cross
0106 cent_x = round(5/7*N);
0107 cent_y = round(2/7*N);
0108 r      = round(13/70*N);
0109 for i = cent_x-r:cent_x+r
0110     for j = cent_y-r:cent_y+r
0111         di = (i-cent_x)/r;
0112         dj = (j-cent_y)/r;
0113       A(i,j) = max(A(i,j),g_min+min([(1-di-dj),(1+di+dj)])*3*delta_g);
0114       A(i,j) = max(A(i,j),g_min+min([(1-di+dj),(1+di-dj)])*3*delta_g);
0115     end
0116 end
0117 
0118 %% torus and box
0119 cent_x = round(5/7*N);
0120 cent_y = round(5/7*N);
0121 b      = round(12/70*N);
0122 R      = round(12/70*N);
0123 r      = round(03/70*N);
0124 rm     = (R+r)/2;
0125 for i = cent_x-b:cent_x+b
0126     for j = cent_y-b:cent_y+b
0127       dx = (i-cent_x);
0128       dy = (j-cent_y);
0129       dr = sqrt(dx^2+dy^2);
0130       ddr = 2*(dr - rm)/(R-r);
0131       if dr >= r && dr <= R
0132           A(i,j) = g_min + 3*delta_g*(1-sqrt(ddr^2));
0133       end
0134     end
0135 end
0136 for i = cent_x-b:cent_x-rm
0137     for j = cent_x-b:cent_y+b
0138         A(i,j) = g_min+3*delta_g;
0139     end
0140 end
0141 
0142 %% add noise
0143 M = (rand(N,N) < outlier_percentage);
0144 A = A + randn(N,N)*sigma + M .*(rand(N,N)-0.5)*8*delta_g;
0145 
0146 %% Check
0147 figure('name','original dem as image','color','w',...
0148     'Position',[0.02*ss(1),0.02*ss(2),0.3*ss(1),0.4*ss(2)]);
0149 imagesc(A)
0150 axis equal;axis off
0151 colormap(gray)
0152 title('original dem as image','FontSize',16)
0153 
0154 %% output
0155 dem = A;
0156 delta_x = 1;
0157 k = 0;
0158 points = zeros(BB(3)*BB(4),4);
0159 out_in = ones(BB(3)*BB(4),1);
0160 for i=BB(1):BB(3)
0161     for j=BB(2):BB(4)
0162         k=k+1;
0163         points(k,:) = [i,j,A(i,j),sigma];
0164         out_in(k) = 1-M(i,j);
0165     end
0166 end
0167 dem = dem(BB(1):BB(3),BB(2):BB(4));
0168 dem = dem(BB(1):BB(3),BB(2):BB(4));
0169 
0170 figure('name','original dem','color','w',...
0171     'Position',[0.02*ss(1),0.52*ss(2),0.3*ss(1),0.4*ss(2)]);
0172 plot_surface(dem,BB,delta_x,'plotfun',@mesh,'view',[-65,29],'alpha',1);
0173 axis equal
0174 title('original dem as mesh','FontSize',16)
0175 return