function vectors = vectorClip( vectors, boundary, is_closed )

% Intersect vectors with boundary
nvectors = numel(vectors);
new_vectors=[];
for s=1:nvectors
    orig_vec = vectors(s);
    assert(~any(isnan(orig_vec.X)));
    
    %Ignore vectors that are completely outside boundary
    inside = isPointInPolygon([orig_vec.X',orig_vec.Y'],boundary);
    if(~any(inside))
        continue;
    end
    inside = inside';
    
    %Get splits
    if(is_closed)
        [orig_vec.X,orig_vec.Y] = poly2cw(orig_vec.X,orig_vec.Y);
        [orig_vec.X,orig_vec.Y] = polybool('intersection',orig_vec.X,orig_vec.Y,boundary(:,1),boundary(:,2));
        if(isempty(orig_vec.X))
            continue;
        end
        splits = find(isnan(orig_vec.X));
        if(~isnan(orig_vec.X(end)))
            splits = [splits,numel(orig_vec.X)+1];
        end
    else
        splits = find(~inside);
        splits = [splits,numel(orig_vec.X)+1];
    end
    
    %Split vectors
    nsplits = numel(splits);
    splits = [0,splits];
    for k=1:nsplits
        vec = orig_vec;
        
        first = splits(k)+1;
        last = splits(k+1)-1;
        
        if(isempty(first:last))
            continue;
        end
        
        if(~is_closed)
            first = max(first-1,1);
            last = min(last+1,numel(vec.X));
        end
        
        vec.X = vec.X(first:last);
        vec.Y = vec.Y(first:last);
        
        vec.BoundingBox(1) = min(vec.X);
        vec.BoundingBox(2) = min(vec.Y);
        vec.BoundingBox(3) = max(vec.X);
        vec.BoundingBox(4) = max(vec.Y);
        
        new_vectors = [new_vectors;vec];
    end
end
vectors = new_vectors;

end


