function [ dx, azi, fractureCount ] = linearscanlineorientation( fractureShp, lineCentre, lineLength, lineAzimuth)
% 
% [ dx, azi, fractureCount ] = linearscanlineorientation( fractureShp, lineCentre, length, azimuth)
% Applies the Linear Scanline fracture sampling method to a
% vectorized fracture map. It retrieves the intersections (dx) and orientations of
% intersected fractures (azi). In addition fractureCount is returned, which
% can be different from length(dx) if a single fracture intersects the
% scanline more than once (ie curved fracture).
% 
% Input:
% The fracture map is specified by fractureShp,
% which must be shapefile structure as produced by the 'shaperead'
% function. Linecentre is a 2x1 vector specifying the centre coordinates
% of the scanline, lenght is a scalar specifying the length of the
% scanline, Azimuth specifies its orientation. 
% 
% Output:
% dx  : Distances from scanline startpoint for each encountered fracture
% azi : Azimuths of fractures intersecting the scanline, orientations
%       are determined at the intersection, assuming the fractureShp is
%       horizontal, and the top corresponds to geographic North.
% fractureCount : Number of unique fractures encountered along the
%       scanline, this will usually equal length( dx ) or length( azi ), except
%       when curved fractures intersect the scanline more than once.
% 
% Created by Yannick Kremer May 2009


index = 1;

% Generate endpoint coordinates for scanline
[sla, slb ] = dir2cos( lineAzimuth, 0);
sl1 = lineCentre - 0.5*lineLength*[sla, slb];
sl2 = lineCentre + 0.5*lineLength*[sla, slb];

dx = [];
azi = [];
fractureCount = 0;

% figure
% plot( [sl1(1) sl2(1)], [sl1(2) sl2(2)], '-r')
% hold on;

for fractureIdx = 1 : length( fractureShp),
   l = length( fractureShp( fractureIdx ).X)-1;
   
   if l>1,
     x = [fractureShp( fractureIdx ).X];
     y = [fractureShp( fractureIdx ).Y];
     intersectingscanline = false;
%      plot(x,y,'b-');
     
     for idx = 1 : l,
        a1 = [x(idx), y(idx)];
        a2 = [x(idx + 1), y( idx + 1 )];
        
        [intersection, onboth] = intersecting( a1, a2, sl1, sl2, false);
        if onboth,          
            intersectingscanline = true;
            dx(index) = sqrt(sum((intersection-a1).^2));
            
            alpha = a2(1) - a1(1);
            beta = a2(2) - a1(2);
            azi(index) = cos2dir( alpha, beta, 0 );
            
            index = index + 1;           
        end
     end
     if intersectingscanline, 
         fractureCount = fractureCount + 1;
     end
   end  
    
end

subplot(1, 3, 1:2)
hold off
mapshow( fractureShp )
hold on
plot( [sl1(1) sl2(1)], [sl1(2) sl2(2)], '-r')

if fractureCount > 0,    
    subplot(1, 3, 3)
    
    
    georose( azi, 36, true), hold on
    plot( [ -sla sla], [-slb slb], '-r')
end

