function [points, d2, d2INN, d2OUT] = computePlaneRAY(data)

[ma, id] = max(abs([data.plane.x, data.plane.y, data.plane.z]));

if id == 1 
    for i=1:2
        for j=1:2
            x(i,j) = (data.plane.d-data.plane.y*data.ym(i)-data.plane.z*data.zm(j))/data.plane.x;
        end
    end
    points(1,:) = [x(1,1) data.ym(1) data.zm(1)];
    points(2,:) = [x(2,1) data.ym(2) data.zm(1)];
    points(3,:) = [x(2,2) data.ym(2) data.zm(2)];
    points(4,:) = [x(1,2) data.ym(1) data.zm(2)];

elseif id == 2
    for i=1:2
        for j=1:2
            y(i,j) = (data.plane.d-data.plane.x*data.xm(i)-data.plane.z*data.zm(j))/data.plane.y;
        end
    end
    points(1,:) = [data.xm(1) y(1,1) data.zm(1)];
    points(2,:) = [data.xm(2) y(2,1) data.zm(1)];
    points(3,:) = [data.xm(2) y(2,2) data.zm(2)];
    points(4,:) = [data.xm(1) y(1,2) data.zm(2)];

else
    for i=1:2
        for j=1:2
            z(i,j) = (data.plane.d-data.plane.x*data.xm(i)-data.plane.y*data.ym(j))/data.plane.z;
        end
    end
    points(1,:) = [data.xm(1) data.ym(1) z(1,1)];
    points(2,:) = [data.xm(2) data.ym(1) z(2,1)];
    points(3,:) = [data.xm(2) data.ym(2) z(2,2)];
    points(4,:) = [data.xm(1) data.ym(2) z(1,2)];
end

ne = [data.plane.x; data.plane.y; data.plane.z];
nn = 0;
pts = [];
pti = [];

for i=1:size(data.afs.kanreal,1)
    x = data.afs.pts(data.afs.kanreal(i,1),:);
    y = data.afs.pts(data.afs.kanreal(i,2),:);
    if sign(x*ne-data.plane.d)*sign(y*ne-data.plane.d) == -1
        nn = nn+1;

        xn = x*ne;
        yn = y*ne;

        al = (data.plane.d-yn)/(xn-yn);

        pts(nn,:) = al*x+(1-al)*y;
    elseif x*ne == data.plane.d
            nn = nn+1;
            pts(nn,:) = x;
    elseif y*nn == data.plane.d
            nn = nn+1;
            pts(nn,:) = y;
    end
end
d2 = pts;

ptsOUT = [];
nn = 0;
for i=1:size(data.afs.kanrealout,1)
    x = data.afs.pts(data.afs.kanrealout(i,1),:);
    y = data.afs.pts(data.afs.kanrealout(i,2),:);
    if sign(x*ne-data.plane.d)*sign(y*ne-data.plane.d) == -1
        nn = nn+1;

        xn = x*ne;
        yn = y*ne;

        al = (data.plane.d-yn)/(xn-yn);

        ptsOUT(nn,:) = al*x+(1-al)*y;
    elseif x*ne == data.plane.d
            nn = nn+1;
            ptsOUT(nn,:) = x;
    elseif y*nn == data.plane.d
            nn = nn+1;
            ptsOUT(nn,:) = y;
    end
end
d2OUT = ptsOUT;

ptsINN = [];
nn = 0;
for i=1:size(data.afs.kanrealinn,1)
    x = data.afs.pts(data.afs.kanrealinn(i,1),:);
    y = data.afs.pts(data.afs.kanrealinn(i,2),:);
    if sign(x*ne-data.plane.d)*sign(y*ne-data.plane.d) == -1
        nn = nn+1;

        xn = x*ne;
        yn = y*ne;

        al = (data.plane.d-yn)/(xn-yn);

        ptsINN(nn,:) = al*x+(1-al)*y;
    elseif x*ne == data.plane.d
            nn = nn+1;
            ptsINN(nn,:) = x;
    elseif y*nn == data.plane.d
            nn = nn+1;
            ptsINN(nn,:) = y;
    end
end
d2INN = ptsINN;
