function [ip,dist2ip,normal] = BR_int_box(faces,normals,pos,p_new)
%

LINE_EPS = 0.001;
dir = p_new - pos;
dir = dir/norm(dir);
pos2 = pos + dir;
line = [pos;pos2];

[ip1,st1] = cv_int_line_plane(line,faces(1:3,:),LINE_EPS);
[ip2,st2] = cv_int_line_plane(line,faces(4:6,:),LINE_EPS);
[ip3,st3] = cv_int_line_plane(line,faces(7:9,:),LINE_EPS);
[ip4,st4] = cv_int_line_plane(line,faces(10:12,:),LINE_EPS);
[ip5,st5] = cv_int_line_plane(line,faces(13:15,:),LINE_EPS);
[ip6,st6] = cv_int_line_plane(line,faces(16:18,:),LINE_EPS);
ips(1,:) = ip1;
ips(2,:) = ip2;
ips(3,:) = ip3;
ips(4,:) = ip4;
ips(5,:) = ip5;
ips(6,:) = ip6;
sts(1) = st1;
sts(2) = st2;
sts(3) = st3;
sts(4) = st4;
sts(5) = st5;
sts(6) = st6;

dists = Inf*ones(6,1);
dirs = zeros(6,3);

for p = 1:6
    if sts(p)==1  % Line not allowed in plane
        dir_p = ips(p,:) - pos;
        dist_p = norm(dir_p);
        dir_p = dir_p/norm(dir_p);
        dir_dot = dot(dir,dir_p);
        if dir_dot>0.9
            dists(p) = dist_p;
            dirs(p,:) = dir_p;
        end
    end
end

[min_dist,min_index] = min(dists);
ip = ips(min_index,:);
dist2ip = min_dist;
normal = normals(min_index,:);
