curve_offset.m

%% Function Declaration
% Source: https://www.mathworks.com/matlabcentral/fileexchange/52496-offset-curve
function [joinedx, joinedy] = curve_offset(x, y, offset)
haxes = [];
    % Offset a curve by a given distance
    % Inputs:
    %       x, y: input x and y coordinates
    %       offset: offset amount in arbitrary units or in points if haxes
    %           is provided
    %       haxes: handle to parent axis if offset is determined in points
    %       intersectremove: remove self-intersecting portions (default is
    %           true)
    %
    % J. Duchateau, - IHU LIRYC, Bordeaux, France - 2015. 
    %
    % You are free to distribute/modify as you please.
    
    % Now offset...
    % 1) Get unit vector size in points
    if nargin < 4 || isempty(haxes), dirvect = [1 1];
    else
        dirvect = vectorToPoints([1 1], haxes);
    end
    
    % 2) Convert vector directions in points
    compvect = dirvect(1) .* diff(x) + 1i*dirvect(2) .* diff(y);
    directions = angle(compvect);
    
    % 3) Rotate by 90°
    directions = directions + pi/2;
    
    % 4) Offset by input offset value
    offsetvect = offset * exp(1i*directions);
    
    % 5) Convert back to axes units
    offx = real(offsetvect) ./ dirvect(1);
    offy = imag(offsetvect) ./ dirvect(2);
    
    joinedx = offx + x(1:end-1);
    joinedy = offy + y(1:end-1);

    % 3) Remove self-intersections
    iPt = 2;
    while iPt < length(joinedx) - 1
        
        % Find candidates for self intersection
        xbnds = sort(joinedx([iPt-1 iPt]));
        ybnds = sort(joinedy([iPt-1 iPt]));
        
        rside = xbnds(2) < joinedx(iPt+1:end); % Fast scan of candidates
        lside = xbnds(1) > joinedx(iPt+1:end);
        above = ybnds(2) < joinedy(iPt+1:end);
        below = ybnds(1) > joinedy(iPt+1:end);
        
        rside = rside(1:end-1) & rside(2:end); % Both on right side
        lside = lside(1:end-1) & lside(2:end);
        above = above(1:end-1) & above(2:end);
        below = below(1:end-1) & below(2:end);
        
        cands = find(~(rside | lside | above | below));
        if ~isempty(cands)
            [xi, yi, flag] = arrayfun(@(x,y) ...
                intersectSeg(joinedx(iPt-1), joinedy(iPt-1), joinedx(iPt), joinedy(iPt), ...
                joinedx(x), joinedy(x), joinedx(x+1), joinedy(x+1)), (iPt+cands));
            iflag = find(flag == 1, 1, 'last');
            if ~isempty(iflag)
                joinedx(iPt) = xi(iflag);
                joinedy(iPt) = yi(iflag);
                joinedx(iPt+1:iPt+cands(iflag)) = [];
                joinedy(iPt+1:iPt+cands(iflag)) = [];
            end
        end
        iPt = iPt+1;
    end
end
function [xi, yi, flag] = intersectSeg(ax1, ay1, ax2, ay2, bx1, by1, bx2, by2)
    eps = 10^-10;
    adx = ax2-ax1; ady = ay2-ay1;
    bdx = bx2-bx1; bdy = by2-by1;
    
    xprod = @(x1, y1, x2, y2) x1*y2 - x2*y1;
    rxs = xprod(adx, ady, bdx, bdy);
    
    if abs(rxs) < eps % Parallel segments
        xi = NaN; yi = NaN; flag = 0;
    else
        t = xprod(bx1-ax1, by1-ay1, bdx, bdy) / rxs;
        u = xprod(bx1-ax1, by1-ay1, adx, ady) / rxs;

        xi = ax1 + t*adx;
        yi = ay1 + t*ady;
        if t>=0 && t <= 1 && u >= 0 && u <= 1
            flag = 1;
        else
            flag = 2;
        end
    end
end
function vect = vectorToPoints(vect,hparent)
    % First step: data units to relative inside axis
    if strcmpi(get(hparent, 'Type'), 'axes')
        xspan = diff(get(hparent, 'XLim'));
        yspan = diff(get(hparent, 'YLim'));
        try 
            vect = vect./[xspan yspan];
        catch
            vect = vect./[xspan;yspan];
        end
    end
    
    % Second step: go up to figure
    while(~strcmpi(get(hparent, 'Type'), 'figure'))
        punit = get(hparent, 'Units');
        ppos = get(hparent, 'Position');
        if strcmpi(punit, 'normalized')
            %Multiply to go to parent
            try
                vect = vect.*ppos([3 4]);
            catch
                vect = vect.*ppos([3;4]);
            end
        else
            % Convert units to points and add up
            error('Please use only normalized units');
        end
        hparent = get(hparent, 'Parent');
    end
    
    % Third step: Convert from relative to figure to paper point
    figunit = get(hparent, 'PaperUnits');
    figsize = get(hparent, 'PaperPosition');
    figsize = convertToPoint(figsize([3 4]), figunit);
    try
        vect = vect.*figsize;
    catch
        vect = vect.*figsize';
    end
end


function outdim = convertToPoint(indim, unit)
% Convert to point from input unit
    switch(lower(unit))
        case 'normalized'
            error('Normalized units are not supported');
        case 'inches'
            outdim = 72*indim;
        case 'centimeters'
            outdim = 72*0.393700787*indim;
        case 'points'
            outdim = indim;
        case 'pixels'
            error('Pixel units are not supported');
        otherwise
            error('Unrecognized input unit');
    end
end