function[lidarResults sampleVegLidar sampleGrdLidar]=sampleLidar(polyVertex, lidarVeg, lidarGrd, widthThreshold)

%this function samples lidar points along the transect within an x by x
%sampling square
% xv=xvertices of the polygon
% yv = yvertices of the polygon

%might want to go back and clear out some variables in this code later....

% ----------------------------- BEGIN code for testing only -----------

%{

%%

%load lidarSampleLoc lidarSampleLoc
load samplePolyVert_StrmCurrent 
load streamVegLidar
load streamGrdLidar
widthThreshold=8;
%sampleLocation = lidarSampleLoc;
lidarVeg = streamVegLidar;
lidarGrd = streamGrdLidar;
polyVertex=(samplePolyVert_StrmCurrent);

%A is a temporary sample set of points to use

% will create a for loop that will go through each sample polygon in the
% matrix


%load Stream
disp('lidar data loaded')

% ----------------------------- END code for testing only ---------------

%}
% removed sample location as it's in polyVertex now. clean up other code...
% it's not used.
%sampleLocation=(polyVertex(:,5:6);

totalSamples = (length(polyVertex(:,1)));
xv=zeros(totalSamples,5);
yv=zeros(totalSamples,5);
%define matrices even if they may be empty...
vegGr2m=[];
vegGr14m=[];
sampleVegLidar=[];
sampleGrdLidar=[];
closestIndexVeg=[];
sampleVegLidar_Gr2m_temp=[];
totalVeg=[];
totalGrd=[];
%preallocate for lidarResults
lidarResults=NaN((length(polyVertex)),15);

%progressbar('Calculating Lidar Sample Values');

%totalSamples=9;
for s=1:totalSamples 
%for s=1:25
    
    percentiles=[];
    inVeg=[];
    inGrd=[];
    sampleVegLidar_Gr2m_temp=[];
    %progressbar(s/totalSamples);
    %disp('calculating lidar for next polygon...')
    %march 31 -- was drawing this box wrong and fixed it -- but it may not
    %work when the polygon is in different orientations... check this...
    
    xv(s,:)=[polyVertex(s,14);polyVertex(s,12);polyVertex(s,10);polyVertex(s,8);polyVertex(s,14)];  
    yv(s,:)=[polyVertex(s,15);polyVertex(s,13);polyVertex(s,11);polyVertex(s,9);polyVertex(s,15)];

    % if the buffer width is small than a threshold, don't bother calculating anything...
    
% commented out Sept 29: if polyVertex(s,4) == 0
    if polyVertex(s,4) < widthThreshold 
    %assign FC values to zeros across the board if the width of the buffer is less than the threshold... don't calculate lidar
        %values
        FCSample=0;
        percentiles= [0 0 0 0];
        totalGrdPoints=0;
        totalVegPoints=0;
        vegGr2m=0;
        vegGr14m=0;

    else

        % find ground and veg lidar points within the sample box
        inVeg = inpolygon(lidarVeg(:,1),lidarVeg(:,2),xv(s,:),yv(s,:));
        inGrd = inpolygon(lidarGrd(:,1),lidarGrd(:,2),xv(s,:),yv(s,:));

        % find index values of lidar points within the sample box
        closestIndexVeg=find(inVeg == 1);
        closestIndexGrd=find(inGrd == 1);


        %build a matrix of points that are within the sample box -
        %preallocate lidar matrices and clear before loop
        totalVeg=length(closestIndexVeg);
        sampleVegLidar_temp=NaN(totalVeg,6);
        totalGrd=length(closestIndexGrd);
        sampleGrdLidar_temp=NaN(totalGrd,6);

    if isempty(closestIndexVeg)
        percentiles=[0 0 0 0];
        totalVegPoints=0;
        %this shouldn't be commented out - i am recommenting this back
        %because otherwise the vegGr2m value is held from the previous run
        %and assigned mistakenly to this round...
        vegGr2m=0;
        vegGr14m=0;
        
    else
        
            for pointsVeg=1:length(closestIndexVeg)
                %sampleVegLidar_temp(pointsVeg,:)=[lidarVeg(closestIndexVeg(pointsVeg),1) lidarVeg(closestIndexVeg(pointsVeg),2) lidarVeg(closestIndexVeg(pointsVeg),3)];
                sampleVegLidar_temp(pointsVeg,:)=[(lidarVeg(closestIndexVeg(pointsVeg),1:3)),polyVertex(s,2:3),polyVertex(s,7)];
            end
        vegGr2m_index = find(sampleVegLidar_temp(:,3) >= 6.56); %find total points taller than 2m to calculate Percentiles
        vegGr2m = length(vegGr2m_index);
        
        vegGr14m=length(find(sampleVegLidar_temp(:,3) >= 4.59));%find total points taller than 1.4m (DBH) for FC
       % vegGr14m = length(vegGr14m_index);
        %create matrix of lidar veg points taller than 2m to calculate
        %overstory height (percentile values)
        
           % if isempty(vegGr2m)% going to change this to if it =0
            if vegGr2m==0    
                sampleVegLidar_Gr2m_temp=[];
            else
                sampleVegLidar_Gr2m_temp=zeros(vegGr2m,3);
                for l=1:vegGr2m
                 sampleVegLidar_Gr2m_temp(l,:) = (sampleVegLidar_temp(vegGr2m_index(l),1:3));
                end
            end
        
        totalVegPoints=length(sampleVegLidar_temp);
    end

%calculate lidar ground points within sample box
    if isempty(closestIndexGrd)
        totalGrdPoints=0;
    else

        for pointsGrd=1:length(closestIndexGrd)
            %sampleGrdLidar_temp(pointsGrd,:)=[lidarGrd(closestIndexGrd(pointsGrd),1) lidarGrd(closestIndexGrd(pointsGrd),2) lidarGrd(closestIndexGrd(pointsGrd),3)];
            sampleGrdLidar_temp(pointsGrd,:)=[(lidarGrd(closestIndexGrd(pointsGrd),1:3)),polyVertex(s,2:3),polyVertex(s,7)];

        end
        totalGrdPoints=length(sampleGrdLidar_temp);
    end

    %might be able to ignore this in the end as i don't need it for
    %anything but plotting... but i can uncomment it if i want to export
    % sampled lidar points later...
        

    % ----------------------------BEGIN USE FOR PLOTTING SAMPLE POINTS
    if s==1  %i think this is bombing because sampleVegLidar might be empty
        sampleVegLidar=(sampleVegLidar_temp);
        sampleGrdLidar=(sampleGrdLidar_temp);
    else
        sampleVegLidar=[sampleVegLidar; sampleVegLidar_temp];
        sampleGrdLidar=[sampleGrdLidar; sampleGrdLidar_temp];
    end
    % ----------------------------END USE FOR PLOTTING SAMPLE POINTS

  if isempty(sampleVegLidar_Gr2m_temp)
      percentiles=[0 0 0 0]; 
      vegGr2m=0;
      FCSample=0;
      vegGr14m=0;
  else
    heights=sampleVegLidar_Gr2m_temp(:,3);
    percentiles= funPrctile(heights,[10,20,70,95]); % march 31 -- changed this to 95 -- may still have a problem with the orientation of the rectangles...
    %note: this function comes with the Matlab stats package. The code for
    %it is included to allow for testing of the script however please
    %purchase the stats package if you plan to use it.
    FCSample=(vegGr14m / (totalGrdPoints+totalVegPoints)); %fixed this to represent DBH height of FC rather than 2m used in the percentile calculation
    %changed: Oct 28 2011
  end
 
end

% results headings
% Name, Transect Number,??, totalGroundPoints, totalVegPoints, VegPts greater than 2m, Fractional
% Cover, 10, 20, 70, 90 (percentiles)


 lidarResults(s,:)=[polyVertex(s,1:7) totalGrdPoints totalVegPoints vegGr14m FCSample percentiles]; %changed this to points>1.4m ht for FC (was gr2m)


end


%InPointsIndex=find(d(:,xlocation) == fourPoints(closest,1) & d(:,ylocation) == fourPoints(closest,2));

%{
%BEGIN PLOT -------------- remove for function  --------------------

%%





NotIn=[lidarVeg(~inVeg,1) lidarVeg(~inVeg,2)];

     for boxes=1:length(xv)
       plot(xv(boxes,:),yv(boxes,:),'-g+','markerSize',4), hold on
     end

%plot(lidarVeg(1:5:end,1),lidarVeg(1:5:end,2),'ro','MarkerFaceColor','r','markerSize',2), hold on
plot(lidarGrd(1:end,1),lidarGrd(1:end,2),'co','MarkerFaceColor','c','markerSize',2), hold on

plot(sampleVegLidar(:,1),sampleVegLidar(:,2),'bo','markerSize',2,'MarkerFaceColor','b'), hold on
plot(sampleGrdLidar(:,1),sampleGrdLidar(:,2),'k*','markerSize',6), hold on

%plot streams
%%
% plot(Stream(1:10,5),Stream(1:10,6),'--bo','LineWidth',1.25), hold on
% plot(Stream(11:24,5),Stream(11:24,6),'--bo','LineWidth',1.25), hold on
%%
%plot(xv,yv,'-k*','MarkerFaceColor','k','markerSize',10), hold on

%plot sample polygon boundaries
%{
plot(polyVertex(:,7),polyVertex(:,8),'k*','MarkerFaceColor','k','markerSize',5), hold on
plot(polyVertex(:,3),polyVertex(:,4),'k*','MarkerFaceColor','k','markerSize',5), hold on
plot(polyVertex(:,5),polyVertex(:,6),'k*','MarkerFaceColor','k','markerSize',5), hold on
plot(polyVertex(:,9),polyVertex(:,10),'k*','MarkerFaceColor','k','markerSize',5), hold on
%}

%plot(sampleGrdLidar_temp(:,1),sampleGrdLidar_temp(:,2),'ko','markerSize',2,'MarkerFaceColor','k'), hold on
%plot(NotIn(1:10:end,1),NotIn(1:10:end,2),'ro','markerSize',2,'MarkerFaceColor','r')
axis square
 
% END PLOT -------------- remove for function  --------------------

%}
