function [p1Clip, p2Clip]=clipLineSegment2Rect(rectMaxMin,p1,p2)
%%CLIPLINESEGMENT2RECT Given a 2D rectangle (axis-aligned) and given the
%           two endpoints of a line segment clip the line segment to the
%           rectangle. If the segment and the rectangle do not overlap,
%           then an empty matrix is returned.
%
%INPUTS: rectMaxMin A length-4 array specifying the maximum and minimum
%                   bounds in x and y. rectMaxMin=[xMin;xMax;yMin;yMax].
%                   The minimum values must be <=the maximum values.
%                p1 A length-2 array specifying the start of the line
%                   segment.
%                p2 A length-2 array specifying the end of the line
%                   segment.
%
%OUTPUTS: p1Clip, p2Clip These are the same p1 and p2 but have been clipped to the
%                start and end of the line segment. If nothing clips, the
%                empty matrices are returned.
%
%The algorithm of [1] is used, without speed improvements that are possible
%from direct substitution, which is mentioned in Section IX of [1]. The
%performance of the algorithm is discussed in [2].
%
%EXAMPLE:
% rect=[-1;1;-1;1];
% numLines=6;
% p1=zeros(2,numLines);
% %Lines entirely inside the box.
% p1(:,1)=[0;0]; p2(:,1)=[0.8;0.8];
% p1(:,2)=[0.6;0]; p2(:,2)=[-0.8;0.8];
% %Lines leaving via each edge of the box.
% p1(:,3)=[0.6;-0.1]; p2(:,3)=[1.6;-0.2];
% p1(:,4)=[-0.6;0.1]; p2(:,4)=[-1.6;0.2];
% p1(:,5)=[-0.1;0.5]; p2(:,5)=[-0.2;2.8];
% p1(:,6)=[0.1;-0.5]; p2(:,6)=[0.2;-1.8];
% 
% figure(1)
% clf
% hold on
% %Draw the bounding box lines, extending them.
% extDist=3;
% plot([-extDist,extDist],[-1,-1],'-k','linewidth',2)
% plot([-extDist,extDist],[1,1],'-k','linewidth',2)
% plot([-1,-1],[-extDist,extDist],'-k','linewidth',2)
% plot([1,1],[-extDist,extDist],'-k','linewidth',2)
% for k=1:numLines
%     %Plot the line segments
%     plot([p1(1,k),p2(1,k)],[p1(2,k),p2(2,k)],'-b','linewidth',4)
%     %Clip the segments
%     [p1Clip,p2Clip]=clipLineSegment2Rect(rect,p1(:,k),p2(:,k));
%     %Plot the clipped segments.
%     if(~isempty(p1Clip))
%         plot([p1Clip(1),p2Clip(1)],[p1Clip(2),p2Clip(2)],'-r','linewidth',2)
%     end
% end
%
%REFERENCES:
%[1] T. M. Nicholl, D. T. Lee, and R. A. Nicholl, "An Efficient  New
%    Algorithm for 2-D Line Clipping: its Development and Analysis", 
%    ACM SIGGRAPH Computer Graphics, vol. 21, no. 4, pp. 253-262, Jul.
%    1987.
%[2] F. Dévai, "Analysis of the Nichioll-Lee-Nicholl Algorithm," in
%    Proceedings of the International Conference on Computational Science
%    and its Applications, Signapore, pp. 726-736, 9-12 May, 2005.
%
%August 2020 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

xMin=rectMaxMin(1);
xMax=rectMaxMin(2);
yMin=rectMaxMin(3);
yMax=rectMaxMin(4);

x1=p1(1);
y1=p1(2);
x2=p2(1);
y2=p2(2);

if(x1<xMin)
    [x1,y1,x2,y2,isVisible]=leftColumn(xMin,yMax,xMax,yMin,x1,y1,x2,y2); 
elseif(x1>xMax)
    [x1,y1]=rotate180c(x1,y1);
    [x2,y2]=rotate180c(x2,y2);
    
    [x1,y1,x2,y2,isVisible]=leftColumn(-xMax,-yMin,-xMin,-yMax,x1,y1,x2,y2);
    
    [x1,y1]=rotate180c(x1,y1);
    [x2,y2]=rotate180c(x2,y2);
else
    [x1,y1,x2,y2,isVisible]=centerColumn(xMin,yMax,xMax,yMin,x1,y1,x2,y2);
end

if(isVisible==false)
    %If the line is not in the box at all.
    p1Clip=[];
    p2Clip=[];
else
    p1Clip=[x1;y1];
    p2Clip=[x2;y2];
end
end

function [x1,y1,x2,y2,isVisible]=leftColumn(xL,yT,xR,yB,x1,y1,x2,y2)
%%LEFTCOLUMN The case where (x1,y1) is to the left of the vertical line
%            x=xL.

if(x2<xL)
    isVisible=false;
elseif(y1>yT)
    [x1,y1,x2,y2,isVisible]=topLeftCorner(xL,yT,xR,yB,x1,y1,x2,y2);
elseif(y1<yB)
    y1=reflectxAxis(y1);
    y2=reflectxAxis(y2);
    [x1,y1,x2,y2,isVisible]=topLeftCorner(xL,-yB,xR,-yT,x1,y1,x2,y2);
    
    y1=reflectxAxis(y1);
    y2=reflectxAxis(y2);
else
    [x1,y1,x2,y2,isVisible]=leftEdge(xL,yT,xR,yB,x1,y1,x2,y2);
end
end

function [x1,y1,x2,y2,isVisible]=topLeftCorner(xL,yT,xR,yB,x1,y1,x2,y2)
%TOPLEFTCORNER The case where (x1,y1) is in the top-left corner.

    if(y2>yT)
        isVisible=false;
    else
        relx2=x2-x1;
        rely2=y2-y1;
        
        topProd=(yT-y1)*relx2;
        leftProd=(xL-x1)*rely2;
        %If (x2,y2) is below the line connecting (x1,y1) and (xL,yT).
        if(topProd>leftProd)
            [x1,y1,x2,y2,isVisible]=leftBottomRegion(xL,xR,yB,x1,y1,x2,y2,relx2,rely2,leftProd);
        else
            [x1,y1]=reflectxMinusy(x1,y1);
            [x2,y2]=reflectxMinusy(x2,y2);
          
            [x1,y1,x2,y2,isVisible]=leftBottomRegion(-yT,-yB,-xR,x1,y1,x2,y2,-rely2,-relx2,topProd);
            
            [x1,y1]=reflectxMinusy(x1,y1);
            [x2,y2]=reflectxMinusy(x2,y2);
        end
    end
end

function [x1,y1,x2,y2,isVisible]=leftBottomRegion(xL,xR,yB,x1,y1,x2,y2,relx2,rely2,leftProd)
%%LEFTBOTTOMREGION The case where (x1,y1) is in the left bottom corner.

if(y2>=yB)
    if(x2>xR)
        %The right-edge intersection.
        y2=y1+(xR-x1)*rely2/relx2;
        x2=xR;
    end
    
    %The left-edge intersection.
    y1=y1+leftProd/relx2;
    x1=xL;
    isVisible=true;
else
    bottomProd=(yB-y1)*relx2;
    
    %If (x2,y2) is below the (x1,y1) to (xR,yB) line.
    if(bottomProd>leftProd)
        isVisible=false;
    else
        if(x2>xR)
            rightProd=(xR-x1)*rely2;
            
            %If (x2,y2) is below the (x1,y1) to (xR,yB) line .
            if(bottomProd>rightProd)
                %The bottom-edge intersection.
                x2=x1+bottomProd/rely2;
                y2=yB;
            else
                %The right-edge intersection.
                y2=y1+rightProd/relx2;
                x2=xR;
            end
        else
            %The bottom-edge intersection.
            x2=x1+bottomProd/rely2;
            y2=yB;
        end
        
        %The left-edge intersection.
        y1=y1+leftProd/relx2;
        x1=xL;
        isVisible=true;
    end
end
end

function [x1,y1,x2,y2,isVisible]=leftEdge(xL,yT,xR,yB,x1,y1,x2,y2)
%%LEFTEDGE (x1,y1) is in the left edge.

    if(x2<xL)
        isVisible=false;
    elseif(y2<yB)
        [x1,y1,x2,y2,isVisible]=p2Bottom(xL,xR,yB,x1,y1,x2,y2);
    elseif(y2>yT)
        y1=reflectxAxis(y1);
        y2=reflectxAxis(y2);
        
        [x1,y1,x2,y2,isVisible]=p2Bottom(xL,xR,-yT,x1,y1,x2,y2);
        
        y1=reflectxAxis(y1);
        y2=reflectxAxis(y2);
    else
        relx2=x2-x1;
        rely2=y2-y1;
        
        if(x2>xR)%The right-edge intersection.
            y2=y1+(xR-x1)*rely2/relx2;
            x2=xR;
        end
        
        %The left-edge intersection.
        y1=y1+(xL-x1)*rely2/relx2;
        x1=xL;
        isVisible=true;
    end
end

function [x1,y1,x2,y2,isVisible]=centerColumn(xL,yT,xR,yB,x1,y1,x2,y2)
%CENTERCOLUMN The case where (x1,y1) is between the left and right
%             boundaries.

    if(y1>yT)
        [x1,y1]=rotate270c(x1,y1);
        [x2,y2]=rotate270c(x2,y2);
        [x1,y1,x2,y2,isVisible]=leftEdge(-yT,xR,-yB,xL,x1,y1,x2,y2);
        [x1,y1]=rotate90c(x1,y1);
        [x2,y2]=rotate90c(x2,y2);
    elseif(y1<yB)
        [x1,y1]=rotate90c(x1,y1);
        [x2,y2]=rotate90c(x2,y2);
        [x1,y1,x2,y2,isVisible]=leftEdge(yB,-xL,yT,-xR,x1,y1,x2,y2);
        [x1,y1]=rotate270c(x1,y1);
        [x2,y2]=rotate270c(x2,y2);
    else
        [x1,y1,x2,y2,isVisible]=inside(xL,yT,xR,yB,x1,y1,x2,y2);
    end
end

function [x1,y1,x2,y2,isVisible]=inside(xL,yT,xR,yB,x1,y1,x2,y2)
%%INSIDE The case where (x1,y1) is inside the rectangle.

if(x2<xL)
    [x1,y1,x2,y2]=p2Left(xL,yT,yB,x1,y1,x2,y2);
elseif(x2>xR)
	[x1,y1]=rotate180c(x1,y1);
	[x2,y2]=rotate180c(x2,y2);
    [x1,y1,x2,y2]=p2Left(-xR,-yB,-yT,x1,y1,x2,y2);
    [x1,y1]=rotate180c(x1,y1);
	[x2,y2]=rotate180c(x2,y2);
elseif(y2>yT)
    x2=x1+(x2-x1)*(yT-y1)/(y2-y1);
    y2=yT;
elseif(y2<yB)
    x2=x1+(x2-x1)*(yB-y1)/(y2-y1);
    y2=yB;
end

isVisible=true;
end

function [x1,y1,x2,y2]=p2Left(xL,yT,yB,x1,y1,x2,y2)
%%P2LEFT

    if(y2>yT)
        [x1,y1,x2,y2]=p2LeftTop(xL,yT,x1,y1,x2,y2);
    elseif(y2<yB)
        [x1,y1]=rotate90c(x1,y1);
        [x2,y2]=rotate90c(x2,y2);
        [x1,y1,x2,y2]=p2LeftTop(yB,-xL,x1,y1,x2,y2);
        [x1,y1]=rotate270c(x1,y1);
        [x2,y2]=rotate270c(x2,y2);
    else
        y2=y1+(y2-y1)*(xL-x1)/(x2-x1); 
        x2=xL;
    end
end

function [x1,y1,x2,y2]=p2LeftTop(xL,yT,x1,y1,x2,y2)
%%P2LEFTTOP

    relx2=x2-x1;
    rely2=y2-y1;
    leftProduct=rely2*(xL-x1);
    topProduct=relx2*(yT-y1);
    
    if(topProduct>leftProduct)
        x2=x1+topProduct/rely2;
        y2=yT;
    else
        y2=y1+leftProduct/relx2;
        x2=xL;
    end
end

function [x1,y1,x2,y2,isVisible]=p2Bottom(xL,xR,yB,x1,y1,x2,y2)
%%P2BOTTOM (x1,y1) is in the left edge; (x2,y2) is not beyond the left
%          edge,  and (x2,y2) is beyond the bottom edge.

    relx2=x2-x1;
    rely2=y2-y1;
    leftProd=(xL-x1)*rely2;
    bottomProd=(yB-y1)*relx2;
    if(bottomProd>leftProd)
        isVisible=false;
    else
       if(x2<=xR)
           x2=x1+bottomProd/rely2;
           y2=yB;
       else
           rightProd=(xR-x1)*rely2;
           if(bottomProd>rightProd)
               x2=x1+bottomProd/rely2;
               y2=yB;
           else
               y2=y1+rightProd/relx2;
               x2=xR;
           end
       end
       y1=y1+leftProd/relx2;
       x1=xL;
       isVisible=true;
    end
end

function [x,y]=reflectxMinusy(x,y)
%%REFLECTXMINUSY This function is from Section VI-2 of [1].
%
%REFERENCES:
%[1] T. M. Nicholl, D. T. Lee, and R. A. Nicholl, "An Efficient  New
%    Algorithm for 2-D Line Clipping: its Development and Analysis", 
%    ACM SIGGRAPH Computer Graphics, vol. 21, no. 4, pp. 253-262, Jul.
%    1987.

t=x;
x=-y;
y=-t;

end

function y=reflectxAxis(y)
%%REFLECTAXIS This function is from Section VI-2 of [1].
%
%REFERENCES:
%[1] T. M. Nicholl, D. T. Lee, and R. A. Nicholl, "An Efficient  New
%    Algorithm for 2-D Line Clipping: its Development and Analysis", 
%    ACM SIGGRAPH Computer Graphics, vol. 21, no. 4, pp. 253-262, Jul.
%    1987.

y=-y;

end

function [x,y]=rotate90c(x,y)
%%ROTATE90C This function is from Section VI-2 of [1].
%
%REFERENCES:
%[1] T. M. Nicholl, D. T. Lee, and R. A. Nicholl, "An Efficient  New
%    Algorithm for 2-D Line Clipping: its Development and Analysis", 
%    ACM SIGGRAPH Computer Graphics, vol. 21, no. 4, pp. 253-262, Jul.
%    1987.

t=x;
x=y;
y=-t;

end

function [x,y]=rotate180c(x,y)
%%ROTATE180C This function is from Section VI-2 of [1].
%
%REFERENCES:
%[1] T. M. Nicholl, D. T. Lee, and R. A. Nicholl, "An Efficient  New
%    Algorithm for 2-D Line Clipping: its Development and Analysis", 
%    ACM SIGGRAPH Computer Graphics, vol. 21, no. 4, pp. 253-262, Jul.
%    1987.

x=-x;
y=-y;
end

function [x,y]=rotate270c(x,y)
%%ROTATE270C This function is from Section VI-2 of [1].
%
%REFERENCES:
%[1] T. M. Nicholl, D. T. Lee, and R. A. Nicholl, "An Efficient  New
%    Algorithm for 2-D Line Clipping: its Development and Analysis", 
%    ACM SIGGRAPH Computer Graphics, vol. 21, no. 4, pp. 253-262, Jul.
%    1987.

t=x;
x=-y;
y=t;

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.
