[Matlab으로 STL파일 다루기] 9. 평면으로 STL 자른 단면 만들기 (2) 두 점과 평면의 교점을 구하는 함수

stl 위에서 triangle을 만들고 있는 인접한 두 점의 모든 조합과, 평면과의 교점을 구해야 합니다. 교점들을 연결하면 단면이 됩니다. 

아래와 같은 단계로, 두 점을 이은 선분과 평면사이의 교점 존재 여부를 판별합니다. 

Step) 두 점 사이에 평면이 있는지 확인

→ 있다면 교점 저장

→ 없다면 교점 없음

두 점 사이에 평면이 있는지 여부를 판별해야 합니다. 두 점을 P1, P2, 평면의 법선벡터를 n, 평면 위의 한 점을 Q라고 놓겠습니다. 

평면 위의 한 점 Q와 각 두점을 연결한 벡터 QP1 과 QP2를 만듭니다. 이 벡터를 법선벡터 n에 투영한 벡터를 벡터 Q1P1 과 Q2P2로 놓겠습니다. 이 두 벡터의 방향이 같다면 평면과의 교점이 없는 것이고, 두 벡터의 방향이 다르다면 평면과의 교점이 있는 것입니다. 

#myfun : projected vector of v1 on v2

function f = pjtvec(v1,v2)

    f = norm(v1)*dot(v1,v2)/(norm(v1)*norm(v2))/norm(v2)*v2;

end

만약 두 점사이에 교점이 있다면 교점을 구해야합니다. 두 점과 평면 사이의 교점을 구하는 함수를 정의해봅시다. 

위그림의 I가 교점인데요. 이 점은 OQ1과 OQ2를 m:n으로 내분한는 점입니다. m:n은 P1Q1의 길이와 P2Q2의 길이의 비와 같습니다. 수식으로 표현하면 아래와 같습니다. 

벡터 합을 이용하여 아래와 같이 변형합니다. 

평면과 두 점을 입력하면 교점을 구해주는 함수를 만들어봅시다. 

#plane expression form : ax+by+cz+d=0 -> [a b c d]

function I=plpointersec(pl,P1,P2)

n=[pl(1) pl(2) pl(3)];

Q=[0 0 -pl(4)/pl(3)];

QP1=P1-Q;

QP2=P2-Q;

Q1P1=pjtvec(QP1,n);

Q2P2=pjtvec(QP2,n);

I=(norm(Q2P2)*(P1-Q1P1)+norm(Q1P1)*(P2-Q2P2))/(norm(Q1P1)+norm(Q2P2));

end

이제 판별과 교점 구하기를 동시에 하는 함수를 만들어봅시다. 교점이 있으면 교점을, 교점이 없다면 NULL을 반환하는 함수입니다.  

#plane expression form : ax+by+cz+d=0 -> [a b c d]

#decide if the intersection exist.

function I=plpointersec2(pl,P1,P2)

n=[pl(1) pl(2) pl(3)];

#find arbitrary point on plane

if pl(3)!=0

Q=[0 0 -pl(4)/pl(3)];

elseif pl(2)!=0

Q=[0 -pl(4)/pl(2) 0];

elseif pl(1)!=0

Q=[-pl(4)/pl(1) 0 0];

end

#are the points on the plane?

if (dot(n,P1)+pl(4))*(dot(n,P2)+pl(4))==0

  if (dot(n,P1)+pl(4))==0

    fprintf("P1 is on the plane\n")

    I=P1;

  elseif (dot(n,P2)+pl(4))==0

    fprintf("P2 is on the plane(I=P2 is returned)\n")

    I=P2;

  end

  

# if not, decide wether the points has intersection. 

else

  QP1=P1-Q;

  QP2=P2-Q;

  Q1P1=pjtvec(QP1,n);

  Q2P2=pjtvec(QP2,n);

  if dot(Q1P1,Q2P2) > 0

    I=NA;

  else 

  I=(norm(Q2P2)*(P1-Q1P1)+norm(Q1P1)*(P2-Q2P2))/(norm(Q1P1)+norm(Q2P2));

end

end

작동을 잘 하는지 확인해봅시다. 

PL=[0 0 1 -2];

p1=[0 0 0];

p2=[0 0 3];

plpointersec2(PL,p1,p2)

z=2인 평면이기 때문에, 교점은 [0,0,2] 가 생겨야 합니다. 

>> plpointersec2(PL,p1,p2)

ans =

   0   0   2

맞게 계산됩니다.