function [ handlesToGridLines ] = dlrGrid(axesHandle,varargin)
%DLRGRID plots a nice gray grid with dashed lines
%   DLRGRID places the grid lines at the preset ticks of the axes. It only
%   draws a grid, if the axes grid is active (i.e., turn on grid before
%   executing DLRGRID). Hence, it is also possible to draw grid lines for
%   single directions (e.g., set(gca,'XGrid','on','YGrid','off')).
%
%   Use h = DLRGRID(gca,lineParameter,value)
%
%   DLRGRID returns a column vector of handles to all segments of the grid lines
%
%   Example
%      grid on;
%      h = DLRGRID(gca,'linestyle',':');
%
%   If you want to remove the grid lines again, execute
%      delete(h);
%      set(gca,'XTickMode','auto','YTickMode','auto');
%
%
%   Thomas Hulin 08.08.2014
%   Copyright 2014 DLR - Institute of Robotics and Mechatronics
%   Oberpfaffenhofen, Germany

%% The BSD 2-Clause License
%
% Copyright (c) 2014, DLR.
% All rights reserved.
% 
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
% 
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
% 
% 2. Redistributions in binary form must reproduce the above copyright
% notice, this list of conditions and the following disclaimer in the
% documentation and/or other materials provided with the distribution.
% 
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
% IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
% CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
% EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
% PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
% PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
% NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
% 
% The views and conclusions contained in the software and documentation are
% those of the authors and should not be interpreted as representing
% official policies, either expressed or implied, of the FreeBSD Project.


% define parameters:
segmentLength_cm = 0.25;
minimumGapLength_cm = 0.15;
epsilon_cm=0.01; % should be approx. half of the line-width

if ~exist('axesHandle')
    axesHandle = gca;
end
oldUnits = get(axesHandle,'Units');
set(axesHandle,'Units','centimeters');
hold(axesHandle,'on');

drawXGrid = strcmpi(get(axesHandle,'XGrid'),'on');
drawYGrid = strcmpi(get(axesHandle,'YGrid'),'on');

set(axesHandle,'XGrid','off','YGrid','off');
xt = get(axesHandle,'XTick');
yt = get(axesHandle,'YTick');
xl = get(axesHandle,'XLim');
yl = get(axesHandle,'YLim');
% turn auto lim-mode off:
set(axesHandle,'XLim',xl);
set(axesHandle,'YLim',yl);
tl = get(axesHandle,'TickLength');
pos_cm = get(axesHandle,'Position');
xCm2Unit = (xl(2)-xl(1)) / pos_cm(3);
yCm2Unit = (yl(2)-yl(1)) / pos_cm(4);

% length of ticks (=offset for grid-lines):
if pos_cm(3)>pos_cm(4)
    offset_x = tl(1)*(xl(2)-xl(1));
    offset_y = tl(1)*pos_cm(3)*yCm2Unit;
else
    offset_x = tl(1)*pos_cm(4)*xCm2Unit;
    offset_y = tl(1)*(yl(2)-yl(1));
end

% distance between two ticks:
%xRatio = (max(xt)-min(xt)) / (xl(2)-xl(1));
%xTickDistance_cm = xRatio*pos_cm(3)/(length(xt)-1);
%yRatio = (max(yt)-min(yt)) / (yl(2)-yl(1));
%yTickDistance_cm = yRatio*pos_cm(4)/(length(yt)-1);
xTickDistance_cm = (max(xt)-min(xt))/(length(xt)-1)/xCm2Unit;
yTickDistance_cm = (max(yt)-min(yt))/(length(yt)-1)/yCm2Unit;
xSegmentLength=segmentLength_cm*xCm2Unit;
ySegmentLength=segmentLength_cm*yCm2Unit;

xNrSegsInBetween = floor(xTickDistance_cm/(segmentLength_cm+minimumGapLength_cm))-1;
yNrSegsInBetween = floor(yTickDistance_cm/(segmentLength_cm+minimumGapLength_cm))-1;

%nrOfGridLines = length(xt)+length(yt);
h=[]; % handles to grid-lines
dx = xt(2)-xt(1);
dy = yt(2)-yt(1);
if abs(xt(1)-xl(1))<epsilon_cm*xCm2Unit
    xt=xt(2:end);
end
if abs(yt(1)-yl(1))<epsilon_cm*yCm2Unit
    yt=yt(2:end);
end

ytSegmentsPos = [yt'-(ySegmentLength/2) yt'+(ySegmentLength/2)];
if (yNrSegsInBetween>0)
    o = linspace(0,dy*(1-1/(yNrSegsInBetween+1)),yNrSegsInBetween+1);
    ytSegmentsPos_org = ytSegmentsPos;
    for j=1:yNrSegsInBetween
        p = ytSegmentsPos_org(1:end-1,:)+o(j+1);
        ytSegmentsPos = [ytSegmentsPos; p];
        % upper inbetween grid line
        l = ytSegmentsPos_org(end,:)+o(j+1);
        if l(1)<yl(2)
            if l(2)>yl(2)
                l(2)=yl(2);
            end
            ytSegmentsPos = [ytSegmentsPos; l];
        end
        % lower inbetween grid line
        l = ytSegmentsPos_org(1,:)-o(j+1);
        if l(2)>yl(1)+offset_y
            if l(1)<yl(1)+offset_y
                l(1)=yl(1)+offset_y;
            end
            ytSegmentsPos = [l; ytSegmentsPos];
        end
    end
end
if ytSegmentsPos(end,2)>yl(2)
    ytSegmentsPos(end,2)=yl(2);
end
if ytSegmentsPos(1,1)<yl(1)+offset_y
    ytSegmentsPos(1,1)=yl(1)+offset_y;
end

xtSegmentsPos = [xt'-(xSegmentLength/2) xt'+(xSegmentLength/2)];
if (xNrSegsInBetween>0)
    o = linspace(0,dx*(1-1/(xNrSegsInBetween+1)),xNrSegsInBetween+1);
    xtSegmentsPos_org = xtSegmentsPos;
    for j=1:xNrSegsInBetween
        p = xtSegmentsPos_org(1:end-1,:)+o(j+1);
        xtSegmentsPos = [xtSegmentsPos; p];
        % upper inbetween grid line
        l = xtSegmentsPos_org(end,:)+o(j+1);
        if l(1)<xl(2)
            if l(2)>xl(2)
                l(2)=xl(2);
            end
            xtSegmentsPos = [xtSegmentsPos; l];
        end
        % lower inbetween grid line
        l = xtSegmentsPos_org(1,:)-o(j+1);
        if l(2)>xl(1)+offset_x
            if l(1)<xl(1)+offset_x
                l(1)=xl(1)+offset_x;
            end
            xtSegmentsPos = [l; xtSegmentsPos];
        end
    end
end
if xtSegmentsPos(end,2)>xl(2)
    xtSegmentsPos(end,2)=xl(2);
end
if xtSegmentsPos(1,1)<xl(1)+offset_x
    xtSegmentsPos(1,1)=xl(1)+offset_x;
end
if drawXGrid==true
    for k=1:length(xt)
        hNew=plot( (ones(length(ytSegmentsPos),1)*[xt(k) xt(k)])',ytSegmentsPos','-','color',[0.8 0.8 0.8],'parent',axesHandle,varargin{:});
        h=[h; hNew];
    end
    % extend axes-limits such that the right and top grip-line is fully printed
    if abs(xt(end)-xl(2))<epsilon_cm*xCm2Unit
        xl(2)=xt(end)+epsilon_cm*xCm2Unit;
        set(axesHandle,'XLim',xl);
    end
end
if drawYGrid==true
    for k=1:length(yt)
        hNew=plot( xtSegmentsPos',(ones(length(xtSegmentsPos),1)*[yt(k) yt(k)])','-','color',[0.8 0.8 0.8],'parent',axesHandle,varargin{:});
        h=[h; hNew];
    end
    % extend axes-limits such that the right and top grip-line is fully printed
    if abs(yt(end)-yl(2))<epsilon_cm*yCm2Unit
        yl(2)=yt(end)+epsilon_cm*yCm2Unit;
        set(axesHandle,'YLim',yl);
    end
end

uistack(h,'bottom'); % bring the grid lines behind all other objects
handlesToGridLines=h;

% prevent changing the axis-ticks:
set(axesHandle,'XTickMode','manual','YTickMode','manual');

% reset units:
set(axesHandle,'Units',oldUnits);

end