trimplot(pt, h, ver)

function trimplot(pt, h, ver)
%TRIMPLOT    Trims 2D Plots
%     TRIMPLOT(PT,HANDLE,VERIFY) Trims the curve segment closest to PT at
%     the nearest intersections and/or ends of the curve on either side of
%     the nearest point on the curve to PT. HANDLE specifies the axes with
%     the multiple 2D curves. VERIFY is a logical that can be set to false
%     if the trim operation need not be validated by the user. HANDLE and
%     VERIFY are both optional and default to the current axes and true
%     respectively. PT is specified as a vector [X,Y].

if nargin == 1
    h = gca;
    ver = true;
elseif nargin == 2
    if islogical(h)
        ver = h;
        h = gca;
    else
        ver = true;
    end
end
trimplotcore(pt,h,ver);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function trimplotcore(pt, h, ver)
pt1=pt(1);
pt2 = pt(2);
pt = [pt1;pt2];
ax = axis();
ch = get(h,'Children');
lns = [];
for ind = 1:length(ch)
    if isequal(get(ch(ind),'Type'), 'line')
        lns = [lns ch(ind)];
    end
end
if length(lns) == 0
    error('No lines to be trimmed')
end

[m,proj,pos] = findclosestline(pt,lns)
trl = lns(m); 
lns(m) = [];
xd = get(trl, 'XData');
xdb = xd;
yd = get(trl, 'YData');
ydb = yd;
if size(pos,2) == 2
    xd = [xd(1:pos(1)), proj(1), xd(pos(2):end)];
    yd = [yd(1:pos(1)), proj(2), yd(pos(2):end)];
    pos = pos + 1;
end    
xdn = xd(1:pos);
ydn = yd(1:pos);
xup = xd(pos:end);
yup = yd(pos:end);
posdn = 0;
posup = 0;
vals = [];
for ind = 1:length(xup)-1
    vals = findintersection([xup(ind), xup(ind+1);yup(ind), yup(ind+1)], lns);
    if length(vals)>0
        posup = ind+1;
        xup = [vals(1),xup(posup:end)];
        yup = [vals(2),yup(posup:end)];
        %plot(xup,yup,'g^');
        break;
    end
end
if posup == 0
    xup = [];
    yup = [];
end
vals = [];
for ind = length(xdn):-1:2
    vals = findintersection([xdn(ind), xdn(ind-1);ydn(ind), ydn(ind-1)], lns);
    if length(vals)>0
        posdn = ind-1;
        xdn = [xdn(1:posdn), vals(1)];
        ydn = [ydn(1:posdn), vals(2)];
        break;
    end
end
if posdn == 0
    xdn = [];
    ydn = [];
end


 
set(trl,'XData', []);
set(trl,'YData', []);
l1 = get(trl);
l2 = get(trl);
l1.XData = xup;
l1.YData = yup;
l2.XData = xdn;
l2.YData = ydn;
intsect = intersect({'DisplayName', 'XDataMode','XDataSource','YDataSource','ZDataSource','BeingDeleted','Type'}, fieldnames(l1));
l1 = rmfield(l1,intsect);
intsect = intersect({'DisplayName', 'XDataMode','XDataSource','YDataSource','ZDataSource','BeingDeleted','Type'}, fieldnames(l2));
l2 = rmfield(l2,intsect);
a = line(l1);
b = line(l2);
axis(ax);
if ver
    usin = input('Would you like to keep the modified plot? (Y/y/N/n)\n','s');
    if ~(isequal(usin,'Y') || isequal(usin, 'y'))
        set(a,'XData',[]);
        set(a,'YData',[]);
        set(b,'XData',[]);
        set(b,'YData',[]);
        set(trl,'XData',xdb);
        set(trl,'YData',ydb);
    end
end
axis(ax);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ret, proj, opos] = findclosestline(pt, lns)
x0 = pt(1,1);
y0 = pt(2,1);
dv = NaN;
opos = 0;
p = [];
pos = 0;
for jnd = 1:length(lns)
    xv= get(lns(jnd),'XData');
    yv= get(lns(jnd),'YData');
    mind = NaN;
    ps = 0;
    for ind = 1:length(xv)-1
        lnpx1 = xv(ind);
        lnpy1 = yv(ind);
        lnpx2 = xv(ind+1);
        lnpy2 = yv(ind+1);
        m = (lnpy2-lnpy1)/(lnpx2-lnpx1);
        c1 = lnpy1 - lnpx1*m;
        c2 = y0 + x0/m;
        projx = (c2 - c1)/(m+(1/m));
        projy = m*projx + c1;
        if (((projx-lnpx1)*(projx-lnpx2)) > 0) || (((projy-lnpy1)*(projy-lnpy2)) > 0)
            dist1 = sqrt((x0-lnpx1).^2 + (y0 - lnpy1).^2);
            dist2 = sqrt((x0-lnpx2).^2 + (y0 - lnpy2).^2);
            [d,loc] = min([dist1,dist2]);
            if loc == 1
                projx = lnpx1;
                projy = lnpy1;
                ps = ind;
            else
                projx = lnpx2;
                projy = lnpy2;
                ps = ind+1;
            end
        else
            d = sqrt((x0-projx).^2 + (y0-projy).^2);
            ps = [ind, ind+1];
        end
        if (d<mind) || isnan(mind)
            mind = d;
            p = [projx;projy];
            pos = ps;
        end
    end
    if (mind < dv) || (isnan(dv))
        dv = mind;
        ret = jnd;
        opos = pos;
        proj = p;
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function ret = findintersection(seg, lns); 
x1 = seg(1,1);
x2 = seg(1,2);
y1 = seg(2,1);
y2 = seg(2,2);
ret = [];
mind = NaN;
for jnd = 1:length(lns)
    xv= get(lns(jnd),'XData');
    yv= get(lns(jnd),'YData');
    int = [];
    for ind = 1:length(xv)-1
        lnpx1 = xv(ind);
        lnpy1 = yv(ind);
        lnpx2 = xv(ind+1);
        lnpy2 = yv(ind+1);
        m1 = (lnpy2-lnpy1)/(lnpx2-lnpx1);
        c1 = lnpy1 - lnpx1*m1;
        m2 = (y2-y1)/(x2-x1);
        c2 = y1 - x1*m2;
        x0 = (c2 - c1)/(m1-m2);
        y0 = m1*x0 + c1;
        if ((((x0-lnpx1)*(x0-lnpx2)) <= 0) && (((x0-x1)*(x0-x2)) <= 0)...
            && (((y0-lnpy1)*(y0-lnpy2)) <= 0) && (((y0-y1)*(y0-y2)) <= 0))
            int = [int, [x0;y0]];
        end
    end
    if length(int) ~= 0
        int
        dist = sqrt((int(1,:)-x1).^2 + (int(2,:)-y1).^2)
        mind
        [d,loc] = min(dist);
        if (d<mind) || isnan(mind)
            mind = d;
            ret = [int(1,loc), int(2,loc)];
        end
    end
end