% Script to analyze shot counts from different LiDAR systems
% Ramon Arrowsmith May and June 2010
clear all
close all

figure(1) %Hillshades
Azimuth=315;
Zenith=45;
Zfact=1;

subplot(2,1,1)
%load up the precomputed 4 m DEM and plot as a hillshade
[X, Y, Easting, Northing, Elevation, dElev_dx, dElev_dy, nulls_count] = LoadFile('4moutput.idw.asc');
Hillshds=hillshademe(Azimuth, Zenith, Zfact, dElev_dy, dElev_dx);
plothillshade(Easting, Northing, Hillshds)
xlabel('Easting (m)')
ylabel('Northing (m)')

subplot(2,1,2)
%load up the precomputed 1 m DEM and plot as a hillshade
[X, Y, Easting, Northing, Elevation, dElev_dx, dElev_dy, nulls_count] = LoadFile('1moutput.idw.asc');
Hillshds=hillshademe(Azimuth, Zenith, Zfact, dElev_dy, dElev_dx);
plothillshade(Easting, Northing, Hillshds)
xlabel('Easting (m)')
ylabel('Northing (m)')

subplot(2,1,1) %show the zoomed area on the overview map
hold  on
xmin = min(min(X));
xmax = max(max(X));
ymin = min(min(Y));
ymax = max(max(Y));
xs = [xmin xmax xmax xmin xmin];
ys = [ymin ymin ymax ymax ymin];
plot(xs,ys, 'k-')

%Analyze and visualize the point count maps
figure (2)
clf

CLOW=0; CHIGH=20; %set the color map range 
subplot(2,2,1)
hold on
[X, Y, Easting, Northing, Elevation, dElev_dx, dElev_dy, nulls_count] = LoadFile('4moutput.den.asc');
show_point_count_map(X(1),Y(1), Easting, Northing, Elevation./(pi*16), CLOW, CHIGH) %show overview point count map
colorbar
plot(xs,ys, 'k-')

subplot(2,2,2)
hold on
compute_and_display_shot_count_histogram(Elevation./(pi*16),200, 20) %compute and display overview point count histogram

subplot(2,2,3)
hold on
[X, Y, Easting, Northing, Elevation, dElev_dx, dElev_dy, nulls_count] = LoadFile('1moutput.den.asc'); 
show_point_count_map(X(1),Y(1), Easting, Northing, Elevation./pi, CLOW, CHIGH) %show zoomed point count map
colorbar

subplot(2,2,4)
compute_and_display_shot_count_histogram(Elevation./(pi),100, 20) %compute and display zoomed point count histogram


% %load all the points and concatenate them
load NCALpoints.txt; %output from OpenTopography point cloud (needs % in first character of first line
newnorthing = NCALpoints(:,2);
neweasting = NCALpoints(:,1);
%elevation = NCALpoints(:,3); %Do not need elevation here
clear NCALpoints

%Now we select small areas and pull the points, count them, and plot them.
%For each of the three small areas, the script wants the user to click on
%the lower left of the box and to then push enter. The script uses that
%location and the range specified in widthofinterest to define the box and
%find the points within the box which are then plotted in the adjacent
%graphics
figure(3)
widthofinterest=15; %width in meters. This makes a 15 x 15 m box
expandaxis = 0.1; %scales the size of of the small maps to be about 10% wider than the box
clf
subplot(2,2,1)
hold on
[X, Y, Easting, Northing, Elevation, dElev_dx, dElev_dy, nulls_count] = LoadFile('1moutput.den.asc');
show_point_count_map(X(1),Y(1), Easting, Northing, Elevation./pi, CLOW, CHIGH)
plot([X(1) X(1)+100], [Y(1)+100 Y(1)+100], 'k-')
axis off
colorbar

[xmin xmax ymin ymax] = location_for_count('1',widthofinterest);

%let's zoom in to box 1
subplot(2,2,2)
hold on
show_point_count_map(xmin,ymin, Easting, Northing, Elevation./pi, CLOW, CHIGH)
find_and_show_points(expandaxis, neweasting, newnorthing, xmin, xmax, ymin, ymax)
axis off

%Let's do another (2)
subplot(2,2,1)
[xmin xmax ymin ymax] = location_for_count('2',widthofinterest);

subplot(2,2,3)
hold on
show_point_count_map(xmin,ymin, Easting, Northing, Elevation./pi, CLOW, CHIGH)
find_and_show_points(expandaxis, neweasting, newnorthing, xmin, xmax, ymin, ymax)
axis off

%Let's do another (3)
subplot(2,2,1)
[xmin xmax ymin ymax] = location_for_count('3',widthofinterest);

subplot(2,2,4)
hold on
show_point_count_map(xmin,ymin, Easting, Northing, Elevation./pi, CLOW, CHIGH)
find_and_show_points(expandaxis, neweasting, newnorthing, xmin, xmax, ymin, ymax)
axis off

