http://blog.naver.com/potheymd/50149717902

매트랩 기초, 매트랩, 매틀랩, 매틀랩기초, 매트랩기초, Matlab, matlab, 수치해석, 수치해석기초

매트랩을 이용해서 데이터 파일을 읽고 쓰는 작업을 자주하게 되므로, 간단하면서도 효율적인 파일 조작에 대해 배워야 한다.

Load 함수

가장 간단한 파일 읽기 방법은 load 함수를 이용하는 것이다.

현재 작업 중인 디렉터리에 올려놓은 텍스트 첨부파일을 저장한다 ('inhaUpop.txt')

'inhaUpop.txt' 에는 다음과 같이 네 줄의 데이터가 동일한 형식으로 저장되어있다.

'inhaUpop.txt'는 2003년부터 2006년까지 인하대학교 학부생 수를 나타내고 있다. 첫번째 열은 년도, 두번째 열은 총 학생 수이다.

'inhaUpop.txt' 를 읽기위해 명령 프롬프트에서 load 명령을 실행한다.

명령을 실행한 뒤에 명령창에는 별다른 내용이 표시되지 않는다.

이 상태에서 who 명령을 수행하면 현재까지 생성된 변수 외에 새로이 inhaUpop가 만들어 졌음을 알 수 있다

변수 inhaUpop의 내용은 아래와 같이 확인 가능하다

이처럼 load 명령을 사용해 읽은 데이터 파일은 파일의 확장자를 제외한 이름의 변수에 그 내용이 저장된다.

그러므로 다음과 같이 변수 이름 inhaUpop를 이용해 그래프를 그릴 수 있다.

fscanf 함수

파일을 읽는 방법으로 load보다 다양한 기능과 옵션을 가진 함수는 fscanf이다. fscanf의 사용법은 다음과 같다

• 여기서 fid는 파일 ID의 의미로 읽은 파일에 붙여지는 ID로 해석하면 된다.
• format은 데이터 형식을 설정해주는 부분이다. fprintf처럼 %d와 %f 등을 사용한다.
• size는 데이터의 크기를 배열 크기처럼 [mXn]으로 정해주는 것이다.

fscanf 함수의 리턴값은 두개인데, 대부분 count없이 a만 명시한다.

결과로 저장된 데이터는 행렬 a로 리턴되고, count 는 성공적으로 읽은 데이터 개수를 말한다.

그런데, 작업 대상인 데이터 파일의 fid를 얻기 위해서는 먼저 fopen 을 사용해 파일을 열어야 한다. fopen의 일반적인 사용법은 다음과같다.

• 여기서 Filename은 읽고자 하는 파일 이름을 의미한다
• Permission은 읽는 파일에 대해 허용된 기능을 의미한다.

일반적으로 읽기(read)만 하고 수정 혹은 첨부하지 않을 때는 'r'을 선택한다. 가능한 permission의 옵션을 보려면 help fopen을 사용한다.

다음 명령으로 inhaUpop.txt의 fid를 부여받을 수 있다.

만약 ()안에 명시한 파일이 존재하지 않을 때는 fid가 -1로 반환된다.

앞의 명령으로 부여받은 fid로 inhaUpop.txt를 배열 a로 저장하는 방법은 다음과 같다.

위의 명령에서 size에 해당하는 [2 4]에 주의해야한다.

fscanf명령은 열(column)-행(row) 순서로 데이터 파일을 읽기 때문에 열의 개수를 먼저 명시해야 한다.

inhaUpop.txt는 두개의 열과 네개의 행으로 되어있으므로 [2 4]로 size를 정했다.

만약 읽어들이는 데이터 파일의 열의 개수는 알지만 행의 개수가 많아서 구태여 개수를 세어 값을 정하지 않기 위해서는

행의 수에 inf를 사용한다. 즉 데이터 파일이 2개의 열로 구성되었다면 size 부분에 [2 inf]를 쓴다.

fscanf의 결과로 저장된 a 행렬은 실제 데이터 파일 형태의 전치행렬이다. 따라서 다음과 같이 a 행렬을 전치시켜 사용하는 것이 일반적인 방법이다.

전치된 행렬을 사용해 load 로 읽어 그린 그래프와 동일한 그래프를 그릴 수 있다.

예제의 경우에는 열린 파일 inhaUpop.txt를 닫을필요는 없지만 매트랩 프로그램이 복잡해질수록,

작업이 종료된 파일은 닫아주는 것이 좋다. 열린 파일을 fid를 이용해 닫는 방법은 다음과 같다

만약 여러개의 파일이 열려있는 상태라면 한 번의 명령으로 모든 파일을 닫는 것이 편ㄹ하다.

이 경우는 fclose all; 을 사용한다.

[Matlab / 매트랩] 강의 3번째

2차원 그래프 그리기 (2D graph)

우왕ㅋ 굿ㅋ 번개같은 업뎃입니다!!!!!!.............(은)는 훼이크고

사실을 금욜날 3강까지 업뎃할 예정이었습니다... 흐흑ㅠ 이미 한참 지나버렸네요 흐흑ㅠ

어찌되었건 은근슬쩍 넘어가서 이번부터는 그래프 그리는걸 한번 배워보도록 해요.

그래프는 2차원 그래프 그리기, 3차원 그래프 그리기, 도형 그리기 요렇게 3 파트로 나눠서 해볼까 하는데

오늘은 먼저 2차원 그래프 그리는 것 부터 해볼까 합니다.

요런 2차원 그래프도 표현이 당연히! 가능합니다.

후후 아름답군요. 엘레강스합니다.

위에 나온 그래프 중에서 가장 빈티나는(...) 2차원 그래프부터 그려보는 걸 한번 해볼까 하는데요.

그래프 그리는 순서는 다음과 같습니다.

1) x축 데이터 정하기(벡터로 표현).
2) y값을 x 값을 이용하여 나타내기.
3) plot명령어를 사용하여 그리기.
4) 추가적인 옵션(그래프 색, grid, x,y이름표 붙이기, 그래프 이름, 범례 표시 등등)

주로 제가 사용하는 방법은 위의 형식입니다. 물론 다른 방법도 있습니다만...! 하나하나 차근히 해보도록 해요.

일단 예를 한번 들어보겠습니다.

y=3*x 그래프를 그려보도록 하겠습니다.

1) x 데이터 지정하기(벡터로 표현)

그리고자 하는 x축 범위를 예를 들어 -3~3 까지 나타내고 싶다고 합시다.
이때 주의할 점은 x축 데이터를 벡터로 표현해야 한다는 것입니다.

왜냐고요?

매트랩이 그래프를 그려주는 원리를 일단 설명을 해야 할 것 같네요.

매트랩이 2차원 그래프를 그리기 위해선 데이터가 2 세트가 필요합니다.

x 데이터와 그에 대응하는 y 데이터죠.

예를 들어 x = [ 1 2 3 4 5] 라는 데이터가 있고
y=[2 4 6 8 10] 이라는 데이터가 있다면

x의 첫번째 데이터와 y의 첫번째 데이터를 하나의 세트로 생각해서 (1,2)를 그래프에 나타내줍니다.
마찬가지로
x의 두번째 데이터와 y의 두번째 데이터를 하나의 세트로 생각해서 (2,4)를 그래프에 나타내줍니다.
역시나 마찬가지로 (3,6),(4,8),(5,10)을 그래프에 찍어주는 형태입니다.

즉, x데이터와 y의 데이터를 순서대로 1:1 대응을 시켜서 그래프를 그려주는거죠.

따라서 x 데이터와 y의 데이터는 벡터의 형태 (1,m 크기의 행렬)로 표시되어야 하며, 두개의 크기가 같아야 합니다.

물론 위의 경우는 점으로 표시되지만, plot명령어를 쓰게되면 그 사이를 선으로 연결시켜 줍니다.

직접 매트랩 명령어와 그래프를 보여 드릴게요.
x와 y를 벡터로 나타낸 다음, plot이란 명령어를 썼습니다.

>> x=[1 2 3 4 5];
>> y=[2 4 6 8 10];
>> plot(x,y)

네 요렇게 나와주시는걸 볼 수 있습니다.

기본적으로 plot명령어는 파란색 선으로 그려주기 때문에 위와 같이 나타나지만 다르게도 표현할 수 있습니다.

자, 여튼 데이터를 나타내기 위해서 벡터 형식으로 나타내야 한다는 것은 아시겠죠?

그런데, -3~3 까지 그리기 위해서 벡터형식으로 나타내려면 숫자를 일일이 쳐줘야 한다는 단점이 있습니다.

x=[-3 -2.99 -2.98 ..... 2.99 3]

요런식으로 둬야 하는데 말이죠.

하지만 현실적으로 일일이 칠 수가 없기 때문에!! 바로 입력하는 방법이 있습니다. 바로 요렇게 말이죠.

x=[시작 값 : 데이터 하나당 크기 : 끝값]

예를 들어 볼게요.

>> x=[-3: 1: 3]

x =
-3 -2 -1 0 1 2 3

쉽죠? -3에서 3까지 1의 단위로 벡터를(혹은 행렬을) 만들어 줘라는 의미가 됩니다.

위에서 x데이터와 y의 데이터를 1:1대응을 해서 나타낸다고 했으므로, 데이터가 촘촘할수록 더 정확한 그래프가 그려지겠죠?

y=3*x를 나타내기 위해 여기선 x의 범위를 -3에서 3까지 0.01씩 지정할게요.

>> x=[-3:0.01:3];

2) y값을 x 값을 이용하여 나타내기

y=3*x의 그래프를 그리고 싶습니다. x는 벡터로 지정되었기 때문에, 이제는 y를 벡터로 지정해줄 차례입니다.

단순히
>> y=3*x;

요렇게 쳐도 괜찮습니다만, 요전 앞 강의에서 행렬의 성분마다 계산!!! 을 하기 위해서는 점!! '.' 을 찍어야 한다는 것. 잊지 않으셨죠?

위의 경우는 숫자 곱하기 벡터라는 단순한 상황이기 때문에 위의 명령어도 먹혀 들어가지만,

2차 그래프나 3차그래프, 분수 등의 다른 형태가 나오게 되면 안되는 경우가 태반이기 때문에,
항상 점 '.'을 찍어주시는 것을 습관화합시다.

즉 요렇게 저는 칠께요.

>> y=3.*x;

이렇게 되면, y는 x라는 벡터에 3을 곱한 결과를 가지는 벡터가 되는 것입니다.

3) plot명령어를 사용하여 그리기
앞에서 이미 설명을 해버렸지만, plot(x,y) 이라는 명령어를 쓰면 x 데이터와 y 데이터를 1:1 대응을 시켜서 그래프를 매트랩이 알아서 그려주게 됩니다. 요렇게 입력해주면 되겠죠?

>>plot(x,y)

여기까지 하면 이런 그래프가 나오게 됩니다.
>> x=[-3:0.01:3];
>> y=3.*x;
>> plot(x,y)

[축하드립니다!! 그래프 그리기 기본 스킬을 습득하셨습니다.]

여기까지 하면 그래프 그리기는 사실상 99% 정도는 끝난거지요.

이제 여기서 그래프를 스리슬쩍 꾸미면 2차원 그래프는 완성되게 됩니다.

4) 추가적인 옵션(그래프 색, grid, x,y이름표 붙이기, 그래프 이름, 범례 표시 등등)

살짜쿵 위의 그래프를 꾸며볼까요? 일단 색깔과 형태를 바꿔 볼게요.
명령어는 다음과 같습니다.
plot(x,y,'색깔')
색깔 부분에는 빨간색이면 'r' 초록색은 'g' 노란색을 'y'등으로 나타낼 수 있어요. 한번 빨간색으로 해볼께요.

>> plot(x,y,'r')

쉽죠? 형태도 바꿀 수 있습니다. 점이나 +, 사각형, 원 등으로 표현할 수도 있어요. 색깔 치는 부분에 형태에 관한 명령어나, '색깔,형태'의 형식으로 나타내주면 됩니다. 한번 보여 드릴게요.
(형태에 관한 명령어 점 -> . 사각형 ->s 원->o 등등. help plot을 치면 쉽게 알 수 있습니다.)

>> plot(x,y,'ro')

데이터 개수가 많아서 원이 겹쳐 굵은 선으로 보이는 군요. 데이터 개수가 적으면 원으로 보입니다 하핫

다음은 다른 옵션들을 적용해볼게요. 순서대로 적용하겠습니다.

x축 이름 표시 (xlabel('이름'))
y축 이름 표시 (ylabel('이름'))
그래프 제목 표시 (title('이름'))
경계 조절(그래프 나타낼 범위 조절) (axis([x최소 x최대 y최소 y최대])
그리드 표시 (grid on)

>> xlabel('x axis')
>> ylabel('y axis')
>> title('graph')
>> axis([-2 2 -8 8])
>> grid on

여러가지 옵션들이 적용된 것을 볼 수 있죠?

pdf 파일

http://web.yonsei.ac.kr/hgjung/Ho%20Gi%20Jung%20Homepage/Lectures/2011%20Fall%20PME/2B%20MATLAB%20%EA%B8%B0%EC%B4%88%20%EC%82%AC%EC%9A%A9%EB%B2%95.pdf

3DMP engines

3D그래픽스 물리 수학, 프로그래밍 GPU Shader 게임엔진 알고리즘 디자인패턴 matlab etc.

회전사원수의 곱을 행렬로 만들다보면 X*Y*Z 또는 Z*Y*X 오일러 회전행렬을 곱한 순서가아니라는 것을 알수있다

회전사원수의 곱은 Z*X*Y 행렬 곱으로 나타난다

[왼손기준 = Z*X*Y 행렬 ]

{{cosz,sinz,0},{-sinz,cosz,0},{0,0,1}}*{{1,0,0},{0,cosx,sinx},{0,-sinx,cosx}}*{{cosy,0,-siny},{0,1,0},{siny,0,cosy}}

[오른손 기준 = Z*X*Y 행렬 ]

{{cosy,0,siny},{0,1,0},{-siny,0,cosy}}*{{1,0,0},{0,cosx,-sinx},{0,sinx,cosx}}*{{cosz,-sinz,0},{sinz,cosz,0},{0,0,1}}

[게임수학] 짐벌락과 쿼터니언 - 1 게임수학 / 개인 2012/08/30 13:15 http://blog.naver.com/qkfkf123/60170135935

아래 내용은 위 내용에 조금 추가하면

위의 이미지는 x 축을 회전 시키는 과정이고 아래 그림은 x 축을 90 또는 -90 도로 회전시켰을때 이고

이때 y축과 z 축이 같은 회전 선상에 있다는 것을 알 수 있다

이때 짐벌락이 발생한다

출처 : http://www.galexandria.com/

처음에

저는 쿼터니언이라는 용어가 마치 비밀스러운 어둠의 힘을 가지고 있는 암흑물체에 관한 양자론 용어처럼 느껴집니다. 여러분도 역시 이 어둠의 힘에 매료되었다면 이 글이 도움이 될 것입니다(그렇게 되기를 바랍니다). 이 글은 쿼터니언을 이용해서 회전을 하는 방법을 보여주어 쿼터니언을 좀더 쉽게 이해하는데 도움을 줄 것입니다. 만약 여러분이 글에서 잘못된 내용을 발견하시면 robin@cyberversion.com로 메일을 보내주세요. 또, 이 기사를 여러분의 사이트에 실으려고 하신다면 저에게 메일을 보내주세요. 저는 저의 글이 어디에 퍼져있는지 알고 싶답니다.

왜 쿼터니언을 사용할까?

이 질문에 대답하기 위해서 우선 방향을 나타내는 일반적인 방법에 대해 논의해 보도록 합시다.

오일러 표현

이것이 현재까지 알려진 방법중 방향을 나타내는 가장 간단한 방법입니다. 오일러 표현은 각각의 축마다 축주위의 회전량을 가지고 있습니다. 따라서, 다음과 같은 0도에서 360도(혹은 0~2π)의 범위에서 변하는 3개의 값을 가집니다.

x, y, z <-- 전역 좌표축 주위의 회전 각도

이 값은 각각 롤, 피치, 요(혹은 피치, 롤, 요, 등 등)를 표현합니다. 방향은 3개의 각도로부터 생성되는 3개의 회전 행렬을 지정한 순서대로 곱해서 구할수 있습니다.

Note: 회전은 전역 좌표계의 좌표축을 기준으로 회전합니다. 즉 최초의 회전이 그 후의 2번째 3번째의 회전축에 영향을 주지 않습니다. 이 때문에 짐벌락(gimbal lock) 문제가 발생하게 됩니다. 짐벌락에 대해서는 조금 후에 자세히 설명하겠습니다.

회전축과 각도(Axis Angle)에 의한 표현

이 방법은 짐벌락을 피할수 있으므로 오일러각 방법보다 좀더 낫습니다. 회전축과 각도에 의한 표현은 임의의 회전축을 나타내는 단위 벡터와 단위 벡터 주위의 회전을 나타내는(0~360) 값으로 구성됩니다.

x, y, z <-- 임의의 축을 나타내는 단위 벡터 angle <-- 바로 윗줄에서 정의한 축 주위로 회전 각도

왜 이 방법이 나쁠까요?

짐벌락

오일러 표현에 있어서 회전은 전역 좌표계에서 일어나기 때문에 한 축의 회전이 다른 축의 회전과 겹치는(override) 문제가 발생합니다. 따라서 회전을 할 수 있는 축 하나를 잃게 됩니다. 이것을 짐벌락이라고 합니다.

예를 들어, X축과 평행한 어떤 벡터를 Y축주위로 회전해서 그 벡터가 Z축과 평행하게 되었다고 하면 이 때, Z축주위로 아무리 회전시켜도, 벡터의 방향은 전혀 변하지 않게 됩니다.

나중에 여러분에게 짐벌락의 예와 쿼터니언을 사용해서 짐벌락을 해결하는 방법을 보여드리겠습니다.

보간 문제

회전축 표현이 짐벌락 문제로 부터 자유롭긴 하지만 두개의 회전을 보간하는 경우, 또 다른 문제가 발생합니다. 보간계산을 마친 회전들이 부드럽지 못하여 회전 애니메이션이 삑사리가 날 수도 있습니다. 오일러 표현도 마찬가지로 이 문제를 가지고 있습니다.

자. 쿼터니언으로 빠져 봅시다 !

시작하기에 앞서, 몇가지 가정을 세우겠습니다. 저는 수학적인 이해가 중요시되는 글에서 수학적인 내용들을 대충 생략해버리는 글들에 진절머리가 납니다. 이런 글들은 독자들을 혼란에 빠뜨리기 쉽상입니다.

좌표계 - 이 글은 OpenGL과 같은 오른손 좌표계를 사용합니다. 만약 여러분이 Direct3D 같은 왼손 좌표계를 사용하고 계시다면 행렬들을 전치(transpose)하셔야 합니다. 이미 Direct3D의 샘플들은 쿼터니언 라이브리러를 가지고 있다는 사실에 유념해주시기 바랍니다. 그럼에도 불구하고 저는 여러분이 그 라이브러리를 사용하시기 전에 그 내부가 어떻게 구성되는지를 한번 짚고 넘어갔으면 하는 바램입니다.

회전순서 - 오일러 표현에서 회전 순서는 X->Y->Z의 순입니다. 행렬 형태로 아래와 같이 표현합니다.

RotX * RotY * RotZ <-- 매우 중요

행렬 - 행렬은 OpenGL 처럼 열우선 방식으로 합니다.

예 [ 0 4 8 12 1 5 9 13 2 6 10 14 3 7 11 15 ]

벡터와 점 - 변환을 위해 벡터와 점은 4x1의 행렬로 표현합니다. 다음과 같은 모습입니다.

회전 행렬 * [ vx vy vz 1 ] <-- 4x1 벡터

저는 특히 Direct3D보다 OpenGL을 선호하지는 않습니다. 그저 제가 OpenGL을 먼저 배웠고, 쿼터니언도 OpenGL을 통해 익혔을 뿐입니다.

Note: 만약 X->Y->Z 순서가 아닌 다른 회전 순서를 지정하면 몇몇 쿼터니언의 함수들을 다시 구현해야 합니다. 특히 오일러 표현을 다루는 함수들이 그렇습니다.

쿼터니언이란 무엇인가?

복소수는 i라는 기호를 사용하여 정의하는 허수(가상의 수)입니다. i는 i * i = -1 라는 성질을 가지고 있습니다.

쿼터니언은 복소수의 확장입니다. i만 사용하는 것이 아니라 제곱근이 ―1 이 되는 3개의 허수를 가집니다. 이 3개의 수는 보통 i, j, k로 표기합니다. 다시 말해 이것은 다음과 같은 성질을 가집니다.

j * j = -1 k * k = -1

따라서 쿼터니언을 아래와 같이 표현할 수 있습니다.

q = w + xi + yj + zk

여기서 w는 실수, x, y, z는 복소수입니다.

흔히 사용하는 또 다른 표현은 아래처럼 벡터형태로 표현할수 있습니다.

q = [ w, v ]

여기서 v = (x, y, z)는 "벡터"라고 말하며 w는 "스칼라"입니다.

v를 벡터라고 부르지만 이것은 일반적인 3차원 벡터가 아닌 4 차원 공간상에 벡터를 표현한 것으로 직관적으로 시각화할 수 없습니다.

항등 쿼터니언

벡터와 다르게 2개의 항등 쿼터니언이 있습니다.

곱 항등 쿼터니언은 아래 처럼 표현합니다.

q = [1, (0, 0, 0)]

그래서 이 곱 항등 쿼터니언과 곱해진 어떤 쿼터니언도 변하지 않습니다.

가산 항등 쿼터니언은 (저희는 사용하지 않습니다)은 아래처럼 표현합니다.

q = [0, (0, 0, 0)]

방향 표현으로 쿼터니언 사용하기

우선 저는 쿼터니언이 벡터가 아니라는 사실을 말씀드리고 싶습니다. 따라서 여기서 벡터 수학을 사용하지 말아주십시요.

이 부분은 매우 수학적인 내용을 다룰 것입니다. 참을성을 가지고 제 설명을 읽어주세요.

우선 쿼터니언의 크기를 정의합니다.

|| q || = Norm(q) = sqrt(w2 + x2 + y2 + z2)

단위 쿼터니언은 아래와 같은 속성을 가지고 있습니다.

w2 + x2 + y2 + z2 = 1

그 때문에 쿼터니언의 정규화는 아래와 같이 구합니다

q = q / || q || = q / sqrt(w2 + x2 + y2 + z2)

이 단위 쿼터니언은 특별합니다. 왜냐하면 단위 쿼터니언은 3D 공간에서 방향을 표현할 수 있기 때문입니다. 따라서 앞에서 논의된 두가지 방법 대신에 단위 쿼터니언을 통해 방향을 표현할 수 있습니다. 단위 쿼터니언으로 방향을 표현하려면 단위 쿼터니언을 다른 표현(예: 행렬)으로 변환하거나 반대로 변환하는 방법이 필요합니다. 이것에 대해서는 곧 설명드리겠습니다.

단위 쿼터니언 시각화 하기

단위 쿼터니언은 (x, y, z) 요소를 임의의 축, w요소를 회전 각도로 하는 4차원 공간상의 회전으로서 시각화할 수 있습니다. 모든 단위 쿼터니언들은 4D 공간에서 단위 길이를 가지는 구를 형성하게 됩니다. 다시한번 이게 무슨 소리지 하고 여러분은 직관적으로 이해가 안되실것입니다. 하지만 제가 정말 말하고 싶은 것은 쿼터니언의 스칼라 요소(w)의 부호를 반전시키는 것만으로 180도 회전한 쿼터니언을 얻을 수 있다는 것입니다.

Note: 단위 쿼터니언만을 방향표현에 사용할 수 있습니다. 이 뒤에 논할 모든 내용은 모두 단위 쿼터니언을 사용한다고 가정합니다.

쿼터니언으로부터의 변환

효과적으로 쿼터니언을 사용하려면 결국 쿼터니언을 다른 표현으로 변환해야 할 필요가 있을 것입니다. 키 눌림을 쿼터니언으로 해석할 수는 없지 않습니까? 할수 있으신 분 계신가요? 글쎄여. 아직까지는 없는 듯 하죠?

쿼터니언에서 행렬로의 변환

OpenGL와 DirectX는 행렬로 회전을 표현하기 때문에 쿼터니언->행렬 변환이 아마 가장 중요한 변환 함수일 것입니다. 왜냐하면 동차 행렬이 3D의 기본 표현이기 떄문입니다.

쿼터니언 회전과 동등한 회전 행렬은 다음과 같습니다.

행렬 = [ w2 + x2 - y2 - z2      2xy - 2wz       2xz + 2wy             0

            2xy + 2wz          w2 - x2 + y2 - z2   2yz - 2wx             0

            2xz - 2wy                2yz + 2wx    w2 - x2 - y2 + z2      0

                  0                            0                   0         w2 - x2 - y2 + z2]

w2 + x2 + y2 + z2 = 1이 되는 단위 쿼터니언의 속성을 사용하면 위 식을 다음과 같이 간단하게 만들 수 있습니다.

행렬 = [ 1 - 2y2 - 2z2          2xy - 2wz         2xz + 2wy            0

            2xy + 2wz            1 - 2x2 - 2z2       2yz - 2wx            0

            2xz - 2wy              2yz + 2wx        1 - 2x2 - 2y2         0

                  0                          0                     0                   1 ]

쿼터니언으로부터 회전축과 각도에 의한 표현으로 변환

한 쿼터니언을 삼차원 공간에서의 임의축 주위의 회전축과 각도에 의한 표현으로 변환하는 방법은 다음과 같습니다.

회전축이 (ax, ay, az)이고 각도가 theta (라디안)이면 각도는 angle= 2 * acos(w)가 됩니다. 그 때 ax= x / scale ay= y / scale az= z / scale 이 됩니다. scale은 scale = sqrt (x2 + y2 + z2)입니다..

제가 알고 있는 다른 방법은 scale = sin(acos(w))을 사용하는 것입니다. 저 스스로 수학적으로 동치관계를 증명하려고 하지는 않았지만 두 방법 다 결과는 같을 것이라 생각합니다.

어쨌든, scale이 0이라면 회전이 없다는 뜻입니다. 그리고 여러분이 특별한 조치를 취하지 않는다면 회전축이 무한대가 됩니다. 따라서 scale이 0인 경우에는 언제나 회전각이 0인 임의의 단위벡터를 그 축에 설정하시면 됩니다.

간단한 예

제가 설명하려고 하는 것들에 대해서 뭐가 뭔지 모르겠다 하시는 독자분들을 위해 이제부터 간단한 예를 보여드리겠습니다.

우선 카메라 방향을 오일러 각으로 표현한다고 해봅시다. 그러면 렌더링 루프에서 다음과 같은 식을 이용하여 카메라를 위치시킵니다.

RotateX * RotateY * RotateZ * Translate

이 때 각각의 요소는 4x4 행렬입니다.

이제 단위 쿼터니언을 사용해서 카메라의 방향을 표현하려면 먼저 쿼터니언을 행렬로 변환해야 합니다. 그러면

(쿼터니언으로부터 변환한 회전 행렬) Rotate * Translate

같은 것이 생기게 됩니다

OpenGL에 특화한 예는 다음과 같게 될것입니다.

위의 표현은 모두 같습니다. 제가 말할려고 하는 것은 방향에 쿼터니언을 사용하는 것은 오일러나 회전축과 각도에 의한 표현과 완전히 같고, 앞서 언급한 변환 함수를 통해 상호교환이 가능하다라는 것입니다.

다만, 위의 쿼터니언에 의한 표현에는 오일러 표현법과 마찬가지로 짐벌락의 위험성이 있습니다.

역자주 – "왜 짐벌락 위험성이 있지"라는 답을 글내용상 알수가 없습니다.

물론, 아직은 회전을 쿼터니언으로 만드는 방법을 모르시겠지만 그것은 아래 부분에서 설명하겠습니다.

Note: Direct3D나 OpenGL를 사용하시는 독자분들은 API가 행렬 연결을 처리해주기 때문에 직접 행렬을 취급하는 일은 없겠지만, 이것에 대해서 알아둘 필요가 있습니다.

쿼터니언 곱

단위쿼터니언이 3D공간에서의 한 방향을 표현하기 때문에 2개의 단위 쿼터니언간의 곱은 2개의 단위 쿼터니언의 회전을 결합한 회전을 나타내는 단위 쿼터니언이 됩니다. 놀라운 일이지만, 사실입니다.

다음과 같은 2개의 쿼터니안이 있다고 가정해봅시다.

Q1=(w1, x1, y1, z1); Q2=(w2, x2, y2, z2);

이 2개의 단위 쿼터니언을 결합한 회전은 아래와 같이 구해집니다

Q1 * Q2 =( w1.w2 - v1.v2, w1.v2 + w2.v1 + v1*v2)

이 때 ,

v1 = (x1, y1, z1) v2 = (x2, y2, z2)

이 됩니다.

여기서 .와 *은 내적과 외적을 나타냅니다.

하지만 위 식을 아래와 같이 최적화할 수 있습니다.

w = w1w2 - x1x2 - y1y2 - z1z2 x = w1x2 + x1w2 + y1z2 - z1y2 y = w1y2 + y1w2 + z1x2 - x1z2 z = w1z2 + z1w2 + x1y2 - y1x2

물론, 여타 쿼터니언들과 마찬가지로 결과 단위 쿼터니언을 다른 회전 표현으로 변환하는 것도 가능합니다. 이런 점이 바로 쿼터니언의 진정한 미학입니다. 2개의 단위 쿼터니언은 구 상에 위치하기 때문에 4차원 공간에서 이들을 곱하는 방법은 짐벌락의 문제를 해결합니다.

여기서 곱하는 순서가 중요합니다. 쿼터니언의 곱은 교환칙이 성립되지 않습니다. 즉 이것은 아래의 식을의미합니다.

q1 * q2 ≠ q2 * q1

Note: 2개의 쿼터니언은 동일한 좌표계 축을 참조해야 합니다. 저는 서로 다른 좌표계에서 2개의 쿼터니언을 합성하는 실수를 범한 적이 있고, 이 때 그 결과 쿼터니언이 특정 각도에서만 실패하는 이유를 알아내기 위해 많은 시간을 고민했습니다.

쿼티니언으로의 변환

이제 다른 표현법을 쿼터니언으로 변환하는 방법을 배워봅시다. 저는 샘플 프로그램에 있는 모든 변환들을 사용하지는 않지만 역운동학 등의 보다 진보된 용도로 쿼터니언 방향을 사용하려고 할 때 이것들이 필요할 것입니다.

회전축과 각도(Angle Axis)로부터 쿼터니언으로 변환

3D공간에서의 임의의 회전축을 도는 회전은 아래처럼 쿼터니언으로 변환됩니다.

회전축이 (ax, ay, az) 이고 (반드시 단위 벡터여야 함) 회전 각도를 theta (라디안)이면 w = cos(theta/2) x = ax * sin(theta/2) y = ay * sin(theta/2) z = az * sin(theta/2)이 됩니다.

우선 회전축이 정규화되어 있어야 합니다. 만약 정규화된 회전축 벡터의 길이가 0이라면 회전이 없는 것을 의미하므로 쿼터니언은 단위(identity) 쿼터니언으로 설정되어야 합니다.

오일러로부터 쿼터니언으로 변환

오일러로부터 쿼터니언으로의 변환은 회전순서가 올바르지 않으면 안 되기 때문에 조금더 힘듭니다. 임의축으로부터 3개의 좌표축을 분해하면 오일러 각을 3개의 독립적인 쿼터니언으로 변환해서 이 3개의 쿼터니언을 곱하면 최종 쿼터니언을 얻을 수 있습니다.

따라서 오일러 각 (a, b, c)으로 세개의 독립적인 쿼터니언을 만들수 있으며

Qx = [ cos(a/2), (sin(a/2), 0, 0)] Qy = [ cos(b/2), (0, sin(b/2), 0)] Qz = [ cos(c/2), (0, 0, sin(c/2))]

최종 쿼터니언은 Qx * Qy * Qz로 구할수 있습니다.

역자주 –이것은 오일러 회전 순서가 x->y->z의 순서인 경우에만 해당됩니다. 다른 순서로 된다면 Qx * Qy * Qz이 다르게 표현됩니다

데모 - 짐벌락 피하기

드디어 "어떻게 쿼터니언이 짐벌락을 피할수 있지?"에 대해서 모두가 기다렸던 궁금증에 대답할 때가 왔습니다.

기본적인 아이디어는

1. 회전을 표현하는 데 쿼터니언을 사용한다.
2. 현재 방향에서 새로운 방향으로 변경하기 위해 하나의 임시 쿼터니언을 생성한다.
3. 임시 쿼터니언과 원래의 쿼터니언을 곱한다. 이렇게 하면 양쪽 모두의 회전을 합성한 새로운 방향이 얻어진다.
4. 이 결과 쿼터니언을 행렬로 변환해서 평소처럼 행렬 곱을 사용한다.

우선, 저는 샘플 코드에 대해서는 어떠한 책임도 지지 않겠습니다. 이 코드는 난잡스러고 제대로 구성되어 있지도 않습니다. 이것은 쿼터니언 테스트할 때 제 프로그램에서 사용했던 쿼터니언 코드를 간략히 줄여놓은 버전일 뿐입니다. 따라서 코드에 대해서 크게 신경을 쓰지 않았다는 점만 기억해 주셨으면 좋겠습니다(돈받고 파는 코드가 아니니까요 ^^)

저는 2개의 실행 가능한 샘플을 준비했습니다. 첫번째 프로그램인CameraEuler.exe (http://www.gamedev.net/reference/programming/features/qpowers/CameraEuler.zip)는 오일러 각을 사용해서 카메라를 구현한 예입니다.

여러분이 주의깊게 봐야 할 부분은 main.cpp의 Main_Lopp 함수입니다.

여러분이 while 문에서 눈여겨 봐야할 곳은 다음과 같습니다.

1. X축, Y축, Z축 각각의 회전을, 3개의 각도로 보관 유지 부분
2. 키를 누르면 관련된 회전량을 조정하는 부분
3. while 문에서, 3 개의 오일러 각도를 회전 행렬로 변환해서 그것을 최종 변환 행렬에 곱하는 부분

위/아래 화살표 키는 X축의 회전, 왼쪽/오른쪽 화살표 키는 Y축회전, Insert/PageUp 키는 Z축회전을 담당합니다.

이 프로그램은, 짐벌락을 일으킵니다. 여러분이 짐벌락 현상을 보고 싶으면, 요를 90도로 취하고 X축이나 Z축회전을 시도해 보십시오. 그리고 무슨일이 일어나는지 살펴보세요.

지금부터는 쿼터니언으로 짐벌락 문제를 해결한 프로그램을 살펴봅시다. 이 프로그램은CameraQuat.exe (http://www.gamedev.net/reference/programming/features/qpowers/CameraQuat.zip)이며 앞의 프로그램을 조금 변경한 프로그램입니다.

여러분이 while문에서 눈여겨봐야하는 곳은 다음과 같습니다.

1. 카메라 방향을 쿼터니언으로 표현한 부분
2. 키 입력으로부터 3개의 각도를 얻는 부분. 이 3개의 각도는 on/off 스윗치이며, 누적되지 않는 다는 것에 주의하십시오. 저는 while문에 이것들을 리셋했습니다. 물론 이것이 최선의 방법은 아닙니다. 그냥 신속히 처리하기 위해서 그랬을 뿐입니다.
3. 3개의 각도를 임시 쿼터니언으로 변환하는 부분
4. 임시 쿼터니언과 카메라 쿼터니언을 곱해서 합성한 방향을 얻는 부분. 곱 순서에 주의.
5. 이렇게 해 얻은 카메라 회전을 최종 변환 행렬을 위해 회전축과 각도에 의한 표현으로 변환하는 부분.

키가 눌러졌을 때, 저는 키 입력에 대응하는 고유의 축에서의 작은 회전을 나타내는 임시 쿼터니언을 생성했습니다. 그리고 임시 쿼터니언과 카메라 쿼터니언을 곱했습니다. 이 4차원 공간에서의 회전을 합성한 것이 짐벌락을 피하게 해줍니다. 여러분이 직접 이것을 해보고 자신의 눈으로 확인해주세요.

카메라 쿼터니언은 최종 변환 행렬로 합성할 수 있도록 행렬이나 그에 상응하는 형태로 변환되야 합니다. 여러분이 쿼터니언을 사용할 때, 3차원 공간과 4차원 공간은 섞일 수가 없기 때문에 항상 이 작업을 해주셔야 합니다. OpenGL의 경우에 저는 그냥 쿼터니언으로부터 회전축을 변경하기만 했고 나머지는 API에게 위임했습니다.

제가 두번째 프로그램에서 오일러 각을 회전에 사용하지 않았지만 저는 여러분이 첫번째 프로그램에서의 오일러 회전에 대해서 볼 수 있게 하기 위해서 오일러 회전 부분을 남겨 두었습니다. 오일러 각은 하나 이상의 축으로 회전시키면 올바르게 되지 않을 것입니다. 왜냐하면 쿼터니언이 카메라 쿼터니언으로부터 오일러 각을 얻는 대신 키 입력을 통해서 계산하기 때문입니다. 두번째 프로그램은 프로그램을 시작했을 때 그냥 여러분이 요 값을 90도로 회전하려는 경우에 짐벌락이 더 이상 일어나지 않는다는 것을 보기 위해서 만든 참고자료일뿐입니다.

Note: 저는 여러분에게 저의 수학 라이브러리를 사용하라고 추천하고 싶지 않습니다. 여러분이 쿼터니언을 이해하고 나서 스스로 짜 보십시오. 참고로 저는 이 코드를 전부 버리고 다시 만들 예정입니다. 이 코드들은 저에게 있어서 너무나 지저분하고 난잡한 코드입니다.

제가 보여드리지 않은 것

여러분이 이미 눈치채셔겠지만 저는 쿼터니언을 오일러 각으로 변환하는 방법에 대해서는 설명하지 않았습니다. 그 이유는 아직까지도 완벽하게 동작하는 변환을 제가 알아내지 못했기 때문입니다. 제가 알고 있는 유일한 방법은 쿼터니언으로부터 행렬을 얻고, 그 행렬로부터 오일러 각을 추출하는 것입니다. 그러나, 오일러에서 행렬 변환은 다대일 관계(sin과 cos으로 인해)이기 때문에 저는 atan2를 사용해서 그 역을 구하는 방법을 알지 못합니다. 만약 정확하게 행렬로부터 오일러 각을 추출하는 방법을 알고 계시다면 저에게 알려주세요.

제가 보여드리지 않은 다른 내용은 행렬을 쿼터니언으로 변환하는 방법입니다. 여러분이 오일러 표현과 회전축과 각도에 의한 표현을 쿼터니언으로 변환할 때 행렬을 거치지 않고도 직접변환이 가능하므로 굳이 행렬을 쿼터니언으로 변환할 필요가 없기 때문입니다.

과제 - SLERP

여러분이 쿼터니언을 완전히 이해했다고 생각하신다면 다시 한번 생각해 보십시오. 아직도 쿼터니언에 대해서 더 배워야 할 게 남아 있습니다. 제가 전에 회전축과 각도(Axis Angle)에 의한 표현이 왜 나쁠까라고 말했던 것을 기억하시나요? '보간'이라는 단어가 갑자기 떠오르지 않나요?

저는 쿼터니언을 사용한 보간에 대해서 설명한 시간을 갖지 못했습니다. 이 글은 제가 예상한 것보다 시간이 오래 걸렸습니다. 여기서 SLERP(구면 선형 보간)에 대한 기본적인 아이디어를 드리겠습니다. 이것은 기본적으로 2개의 쿼터니언 방향 사이에 일련의 쿼터니언들을 생성합니다. 일련의 쿼터니언들은 처음과 마지막 쿼터니언사이에서 부드러운 모션(움직임)의 결과를 제공합니다. 오일러와 회전축 표현으로는 일관성있게 이러한 부드러운 움직임을 만들 수 없습니다.

마지막으로

저는 이 기사가 쿼터니언 이론 뒤에 숨겨진 알수 없는 미스테리들을 속시원하게 없애주었길 바랍니다. 다시한번 마지막으로 여러분에게 당부하고 싶은 말은 서로 다른 좌표계에서 2개의 쿼터니언을 곱하지 말아주십시오. 이렇게 하시면 여러분은 고통의 시간을 맛보시게 될것입니다.

그러면 새로 찾아낸 쿼터니언의 능력에 여러분이 기뻐하시길 바라면서 저는 이제 그만 여러분과 헤어질까 합니다. 몸조심하시구여. 그리구 다시 만나뵙기를 바라면서...

ttp://cs.sungshin.ac.kr/~hkim/LECTURE/CG/Shading/RayTracing.htm

http://blog.naver.com/rlawldyd0806/50036277924

행렬식

-행과 열이 같은 정방 행렬에서만 구할 수 있다

여인자를 이용한 행렬식 표현

동일한 열 성분(a_1j, a_2j, •••, a_nj)에 대해서 표현

동일한 행 성분(a_j1, a_j2, •••, a_jn)에 대해서 표현

1X1정방 행렬  A=[a]의 행렬식 |A|=a

 Direct3D 함수

FLOAT D3DXMatrixDeterminant(CONST D3DXMATRIX *pM);

사용자 함수

//행렬의 행렬식을 구한다.

float MYMatrixDeterminant(D3DXMATRIX *pM)
{
 return (pM->_11*pM->_22*pM->_33) + (pM->_12*pM->_23*pM->_31) + (pM->_13*pM->_21*pM->_32)
      - (pM->_11*pM->_23*pM->_32) - (pM->_12*pM->_21*pM->_33) - (pM->_13*pM->_22*pM->_31);
}

http://blog.naver.com/astone94/90101243937

에르미트(Hermite)[1]곡선과 베지에(Bezier)[2]곡선은3차 곡선의 두 예입니다. 이들은 네가지 매개변수에 의하여 정의됩니다. 에르미트 곡선은 두 끝점과 이 점에서의 탄젠트 벡터(tangent vectors)로 정의됩니다. 베지에 곡선은4개의 점으로 정의됩니다. 이 둘은 수학적으로 다르지만 그 둘의 특징과 한계는 매우 유사합니다. 

---

[1]더 많은 정보를 위해서는 http://en.wikipedia.org/wiki/Cubic_Hermite_spline를 참조하세요.

[2]더 많은 정보를 위해서는 http://en.wikipedia.org/wiki/B%C3%A9zier_curve를 참조하세요.

대부분의 경우, 곡선은 여러가지 요소들에 의하여 만들어집니다. 이 요소들은 구간적3차 곡선(piecewise Bezier curve)을 만들기 위해서 필요합니다. 아래는 구간적 베지에 곡선의 예로  10개의 점(storage point)을 이용하여3개의 구간을 가진 곡선(a three-segment curve)를 나타낸 것입니다.이 곡선들이 아래 그림과 같이 연결(join)되어도 이것이 부드럽게 연결되지는 않습니다.

에르미트 곡선이 베지어 곡선과 같은4개의 매개변수를 사용한다고 하여도, 탄젠트라는 추가적인 방향 정보를 가지고 있습니다. 에르미트 곡선들을 연결할 때 이러한 방향정보를 이용하여 각 끝점의 방향을 일치키시면 적은 용량으로도 더욱 부드러운 곡선을 표현할 수 있습니다.

[1]더 많은 정보를 위해서는 http://en.wikipedia.org/wiki/B%C3%A9zier_curve를 참조하세요.

더욱 부드럽고 연속적인 곡선을 표현하기 위해 찾아진 방식 중 가장 강력한 것으로, Non Uniform Rational B-Spline[1](NURBS)가 있습니다. 곡선 조각들은 적은 정보량으로 더욱 부드러운 곡선을 포현하기 위해 더 많은 컨트롤 포인트를 공유하게 됩니다.

NURBS 곡선과 곡면(surface)들은Rhinoceros에서 기하체를 표현하기 위해 사용됩니다. NURBS곡선들의 특징과 요소들은 이 챕터의 뒷 부분을 통해 더욱 자세하게 다뤄질 것입니다.

3차 곡선의 매개변수 방정식은 다음과 같습니다. 일반적인 모델링 작업에서 이러한 공식을 사용할 가능성은 거의 없지만, 이것을 참고적으로 이해하는 것은 도움이 될 것입니다.

3차 곡선의 매개변수 방정식은 다음과 같습니다.

이를 아래와 같이 바꿔 쓸 수 있습니다.

x(t) = axt ³+ bxt ²+ cxt + dx

y(t) = ayt ³  + byt ²+ cyt + dy

z(t) = azt ³+ bzt ²+ czt + dz

Q(t) 방정식은 아래와 같이 써질 수 있습니다.

Q(t) = C. T

T가 다음과 같고.

C가 계수 행렬(the matrix of coefficients)이라고 한다면

 

[1]더 많은 정보를 위해서는 http://en.wikipedia.org/wiki/Non-uniform_rational_B-spline를 참조하세요.

행렬의 곱을 이용한 곡선 방정식은 다음과 같습니다.

 

[출처] 3차 곡선(Cubic polynomial curves )|작성자 U 기찬 BIM

ttp://en.wikipedia.org/wiki/B%C3%A9zier_curve#Cubic_B.C3.A9zier_curves

브레슨햄 동영상 강의

http://forum.falinux.com/zbxe/?document_srl=406146

http://blog.naver.com/nowcome/130019875263

#include <iostream>
using namespace std;

#define max( x, y) (((x) > (y)) ? (x) : (y)) // 큰수 찾는거 당연히 아시겠줘? 

#define abs( x ) (((x) > 0) ? (x) : -(x)) // 절대값. 

#define sign( x ) ((x) > 0 ? 1 : ((x) == 0 ? 0: (-1))) // 양수:1, 0, 음:-1 

#define WIDTH   30
#define HEIGHT  25

char array[HEIGHT][WIDTH]={{'\0',},};

void show(){

 for (int j=0;j<HEIGHT;++j)
 {
  for (int i=0;i<WIDTH;++i)
  {
   cout<<array[j][i]<<" ";
  }
  cout<<"\n";
 }
 

}

void initial(){
 for (int j=0;j<HEIGHT;++j)
 {
  for (int i=0;i<WIDTH;++i)
  {
   array[j][i]=' ';
  }
 }
}

void drawPoint(int x,int y){

 array[y][x]='*';

}

void line(int x1,int y1,int x2,int y2){

 int x,y;
 int dx = x2 - x1;  
 int dy = y2 - y1;  
 int ix = abs(dx);  
 int iy = abs(dy);  
 int inc = max(ix, iy);  
 int plotx = x1;  
 int ploty = y1; 

 x = y = 0;  
 for(int  I = 0; I <= inc ; I++)  
 {  
  x += ix;  
  y += iy;  
  if( x > inc)  
  {  
   x -= inc;  
   plotx += sign(dx);  
  }  
  if( y > inc )  
  {  
   y -= inc;  
   ploty += sign(dy);  
  }  
  drawPoint(plotx, ploty);
  //그리는 함수 ( plotx, ploty); 

 }
}

int main(){

 initial();

 //line( 0,0,10,10 );

 show();

 return 0;
}

http://withwani.springnote.com/pages/159993

브레슨햄 알고리즘

컴퓨터 화면에 선을 빠르게 그리는 방법 중에 하나가 브레슨햄(Bresenhan) 알고리듬이 있습니다. 실수 나누셈을 사용하는 직선 공식보다 정수 계산만을 하는 브레슽핸 알고리듬이 매우 빠르죠.

이번에 임데디드 포럼에 그래픽 관련 문서를 만들면서 브레슨햄 선긋기 알고리듬을 정리해 보았습니다. 포럼에는 동영상까지 올렸으므로 관심있는 분은 방문해 보시기 바랍니다.

브레슨햄(Bresenham)알고리듬

직선의 공식을 사용하여 머릿 속으로 그래프를 그려 보면, 무수히 많은 점을 생각할 수 있읍니다만, 컴퓨터 그래픽에서는 스크린 좌표가 모두 정수이므로, 복잡하고 계산을 느리게 만드는 실수 계산을 정수 계산 만으로도 구현할 수 없을까하는 생각에서 만들어진 것이 브레슨햄 알고리듬입니다.

직선의 공식을 이용하여 계산된 좌료값이 소주점 이하 몇개의 자리로 된 실수값을 어렵게 구했다고 하더라도 컴퓨터화면에 출력할 때에는 소수점이하를 모두 버림한단든지 반올림해서 정수로 만들게 됩니다. 이렇게 버려지는 소수점 이하의 복잡한 계산을 브레슨햄 공식을 이용하여 없애 보도록 하겠습니다.

*** 이하의 설명은 dx 가 dy 보다 큰 경우일 때 입니다. ***

컴퓨터 화면에 직선을 긋기위해 이미 (x_i, y_i) 좌표에 점을 찍었습니다. 브레슨햄의 기본 원리는 선의 공식에 따라 다음 좌표값을 구하는 것이 아니라 현재 점의 위치에 대해서 어느 좌표에 다음 점을 찍을 지를 판단하는 알고리듬입니다.

그림에서 (x_i,y_i)의 다음 점의 위치는 녹색점(x_i+1, y_i+1) 보다는 빨강색점(x_i+1, y_i) 가 실제 점에 가깝습니다. 그렇다면 여기서 녹색점과 빨강색점에 대해 실제 직선의 위치와의 차이 d1, d2를 구해 보겠습니다.

d1, d2를 구하는 이유는, d1 < d2 라면 실제 y 점이 y_i에 가깝다는 얘기가 되고, y_i 에 점을 찍으면 됩니다. 반대로 d1 > d2 라면 y_i+1에 점을 찍으면 되기 때문입니다.

즉, 0 < d1-d2 이면 y_i+1 에, 0 > d1-d2 이면 y_i 에 점을 찍으면되므로 (1) d1 과 d2의 차, 특히 (2) 차이값이 양수인지, 음수인지를 구하면 되겠습니다.

d1 = y -y_i

여기서 y=m(x_i+1) +b 이므로,

d1 = m(x_i+1) +b -y_i

d2 는,

d2 = (y_i+1) -y
   = y_i +1 -( m(x_i+1)+b)
   = y_i -m(x_i +1) -b +1

두 거리의 차이를 구해 보겠습니다.

d1 - d2 = m(x_i+1) +b -y_i -(y_i -m(x_i +1) -b +1)
        = m(x_i+1) +b -y_i -y_i +m(x_i+1) +b -1
        = 2 m(x_i+1) -2y_i +2b -1

m 은 기울기를 dy/dx 이므로,

  d1 - d2 = 2(dy/dx)(x_i+1) -2y_i +2b -1

여기서 다음 점을 찾기 위해서는 두 거리의 차이가 궁굼한 것이 아니라 두 차이에 대한 부호, 즉 양수인지 음수인지만이 필요하므로 정수 나눗셈이 있는 부분을 양편에 dx를 곱해서 없애겠습니다.

dx(d1-d2) = 2dy(x_i+1) -2dxy_i +2dxb -dx
          = 2dyx_i +2dy -2dxy_i +2dxb -dx

다시 이 공식에서 변수 부분인 x_i 또는 y_i 가 들어간 부분과 한번 계산되면 다시 계사할 필요가 없는 부분으로 정리해 보겠습니다. 그리고 부호값 만을 생각할 수 있는 변수 p_i 로 정리해 보겠습니다.

dx(d1 - d2) = p_i = 2dyx_i +2dy -2dxy_i +2dxb -dx
            = 2dyx_i -2dxy_i +2dy +2dxb -dx
            = 2dyx_i -2dxy_i +2dy +dx(2b -1)

상수부분인 2dy +dx(2b-1)를 C 로 대치하면,

p_i = 2dyx_i -2dxy_i +C   <--- 공식 1

가 됩니다.

우리가 이 알고리듬을 이용해서 선을 긋게 되면, 루프를 이용하게 될 것이고, 이전에 계산된 p_i 값을 이번 계산에 대입하고 적용할 수 있습니다.

그러므로 다음 p_i 값인 p_i+1 값을 아래와 같이 구할 수 있습니다.

p_i+1 = 2dyx_i+1 -2dxy_i+1 +C

이전 p 값을 계속 참조하여 계산 하므로 p_i 와 p_i+1 값의 차를 구하겠습니다.

p_i+1 - p_i= 2dyx_i+1 -2dxy_i+1 +C - 2dyx_i +2dxy_i -C
        = 2dy(x_i+1-x_i) - 2dx(y_i+1-y_i)

여기서 x_i+1 와 x_i는 서로 이웃하는 점이므로 x_i+1 - x_i는 1입니다.

그러므로 아래와 같이 정리할 수 있습니다.

p_i+1 - p_i= 2dy - 2dx(y_i+1-y_i)

이제 p_i+1 값만 구한다면,

p_i+1 = p_i + 2dy - 2dx(y_i+1-y_i)

아직까지 복잡하게 보이지만 역시 스크린 좌표는 정수이고 y_i+1와 y_i 의 차이가 0 또는 1이라는 점을 상기한다면,

• 계산된 pi 값이 양수라면 d1이 더 크므로 다음 점은 현재의 yi 보다 1 만큼 크게되므로,
  
  pi+1 = pi + 2dy - 2dx
  
  가 되며,
  
• pi 값이 음수라면 d1이 더 작으므로 yi 값과 같아지므로 yi+1 - yi 는 0이 되므로,
  
  pi+1 = pi + 2dy
  
  값이 됩니다.

여기서 계산된 pi+1은 다음 루프의 pi 값이 되므로 이를 루프에서 계속 계산되고 다음 계산에 대입된다면 처음 p₀ 값을 구하고 점을 찍고 다음 pi 값으로 대입하면서 점 위치를 찾아 점을 찍으면 됩니다.

p₀ 구하기

p₀값을 구하기 위해 공식 1 에 처음 좌표 (x₁, y₁)을 대입하고, 다음 예상 좌표를 (x₁+1, y₁+1/2)로 생각하고 p₀을 구합니다..

p(x₁, y₁) = 2dyx₁ -2dxy₁ +C

p(x₁+1, y₁+1/2)= 2dy(x_i+1) -2dx(y_i+1/2) +C

p₀  = p(x₁+1, y₁+1/2) - p(x₁, y₁)
    = 2dy(x_i+1) -2dx(y_i+1/2) +C - ( 2dyx₁ -2dxy₁ +C)
    = 2dyx₁ +2dy -2dxy₁ -dx -2dyx₁ +2dxy₁ +C -C
    = 2dy(x₁-x₁) +2dy -2dx(y₁-y₁) -dx
    = 2dy -dx

p₀ 을 구하기 위해 x₁ 의 다음 좌표는 당연히 1 을 증가한 x₁+1 이 되겠지만 , y₁ 축의 값은 y₁ 이 될지 y₁+1이 될지 모르므로 확률 상의 경우를 1/2로 보고 계산한 것입니다.

선을 긋기 전에 필요한 상수 값

• 최초의 p₀ 을 구합니다.
• pi 의 부호 값에 따라 아래의 2가지 값이 필요합니다.
  • 2(dy -dx)
  • 2dy

이후 선을 긋는 루프에서

• 계산된 p 값이 양수이면,
  • y 좌표값을 1 증가하고,
  • p 값에 2(dy-dx) 값을 더합니다.
• 음수 이면,
  • p 값에 2dy 만 더합니다.

          #include <stdio.h>
          #include <stdlib.h>
          #include <string.h>        // abs
          
          #include <unistd.h>        // open/close
          #include <fcntl.h>         // O_RDWR
          #include <sys/ioctl.h>     // ioctl
          #include <sys/mman.h>      // mmap PROT_
          #include <linux/fb.h>
          
          int   screen_width;
          int   screen_height;
          unsigned short *fb_mapped;
          
                                   
          void  dot( int x, int y)
          {
             unsigned short *ptr;
          
             ptr   = fb_mapped + screen_width * y + x;
             *ptr  = 0xffff;
          }
          
          void  line( int x1, int y1, int x2, int y2)
          {        
             int      dx, dy;
             int      p_value;
             int      inc_2dy;
             int      inc_2dydx;
             int      inc_value;
             int      ndx;
                             
             dx       = abs( x2 -x1);
             dy       = abs( y2 -y1);
             
             
             if ( dy <= dx)
             {
                inc_2dy     = 2 * dy;
                inc_2dydx   = 2 * ( dy - dx);
                
                if ( x2 < x1)
                {
                   ndx   = x1;
                   x1    = x2;
                   x2    = ndx;
                   
                   ndx   = y1;
                   y1    = y2;
                   y2    = ndx;
                }
                if ( y1 < y2)  inc_value   = 1;
                else           inc_value   = -1;
                   
                dot( x1, y1);
                p_value  = 2 * dy - dx;    
                for (ndx = x1; ndx < x2; ndx++)
                {
                   if ( 0 > p_value)    p_value  += inc_2dy;
                   else
                   {
                      p_value  += inc_2dydx;
                      y1       += inc_value;
                   }
                   dot( ndx, y1);
                }
             }
             else
             {
                inc_2dy     = 2 * dx;
                inc_2dydx   = 2 * ( dx - dy);
                
                if ( y2 < y1)
                {
                   ndx   = y1;
                   y1    = y2;
                   y2    = ndx;
                   
                   ndx   = x1;
                   x1    = x2;
                   x2    = ndx;
                }
                if ( x1 < x2)  inc_value   = 1;
                else           inc_value   = -1;
                   
                dot( x1, y1);
                p_value  = 2 * dx - dy;    
                for (ndx = y1; ndx < y2; ndx++)
                {
                   if ( 0 > p_value)    p_value  += inc_2dy;
                   else
                   {
                      p_value  += inc_2dydx;
                      x1       += inc_value;
                   }
                   dot( x1, ndx);
                }
                
             }
             
          }  
          
          
          int   main( int argc, char **argv)
          {
             int      fb_fd;
             struct   fb_var_screeninfo  fbvar;
             struct   fb_fix_screeninfo  fbfix;
             int      bytes_per_line;
             int      mem_size;
          
             if ( access( "/dev/fb", F_OK))
             {
                printf( "프레임 버퍼 장치가 없습니다.\n");
                return 0;
             }
          
             if ( 0 > ( fb_fd = open( "/dev/fb", O_RDWR)))
             {
                printf( "프레임 버퍼에 접근할 수 없습니다.\n");
                return 0;
          
             }
          
             if ( ioctl( fb_fd, FBIOGET_VSCREENINFO, &fbvar))
             {
                printf( "FBIOGET_VSCREENINFO를 실행하지 못했습니다.\n");
                return 0;
             }
             if ( ioctl( fb_fd, FBIOGET_FSCREENINFO, &fbfix))
             {
                printf( "FBIOGET_FSCREENINFO 실행하지 못했습니다.\n");
                return 0;
             }
          
             screen_width   = fbvar.xres;                    // 스크린의 픽셀 폭
             screen_height  = fbvar.yres;                    // 스크린의 픽셀 높이
             bytes_per_line = fbfix.line_length;             // 한개 라인 당 바이트 개수
          
             mem_size       = bytes_per_line * screen_height;
          
             fb_mapped      = ( unsigned short *)mmap( 0,
                                                    mem_size,
                                                    PROT_READ|PROT_WRITE,
                                                    MAP_SHARED,
                                                    fb_fd,
                                                    0);
          
             
             if ( 5 != argc)
             {
                printf( "사용방법: $]./sample x1 y1 x2 y2 \n");
             }
             else
             {
                line( atoi( argv[1]), atoi( argv[2]), atoi( argv[3]), atoi(argv[4]));
             }
          
             munmap( fb_mapped, mem_size);
             close( fb_fd);
          
             return 0;
          }
          

void bs_BresenhamLine(vector2i s, vector2i e, CImage* img, COLORREF color) {  vector2i vStart = s;  vector2i vEnd = e;  int dx = vEnd.x-vStart.x;  int dy = vEnd.y-vStart.y;  int stepx = 1;  int stepy = 1;  if(dx<0)  {   dx = -dx;   stepx = -1;  }  if(dy<0)  {   dy = -dy;   stepy = -1;  }  int x = vStart.x;  int y = vStart.y;  int eps = 0;  if(dx>dy)  {   do    {    //m_vLinelist[m_nLineSize++] = vector2f((float)x,(float)y);    img->SetPixel(x,y,color);    if(dx==0)     break;    if((eps+=dy) << 1 >= dx)    {     y += stepy;     eps -= dx;    }   } while((x+=stepx)!=vEnd.x);  }  else  {   do   {    //m_vLinelist[m_nLineSize++] = vector2f((float)x,(float)y);    img->SetPixel(x,y,color);    if(dy==0)     break;    if((eps+=dx) << 1 >= dy)    {     x += stepx;     eps -= dy;    }   } while((y+=stepy)!=vEnd.y);  } }     http://blog.naver.com/bbosasi8/110015833947 브레슨햄 알고리즘  프로그래밍이야기  2007/03/26 11:08 http://blog.naver.com/bbosasi8/110015833947

        for( ;x1<=x2;x1++) 
        { 
                errterm += delta_y; 
                 
                if ( errterm > halfx ) 
                { 
                        errterm -= delta_x; 
                        y1++; 
                } 
        } 

#define max( x, y) (((x) > (y)) ? (x) : (y)) // 큰수 찾는거 당연히 아시겠줘?  
 
#define abs( x ) (((x) > 0) ? -(x) : (x)) // 절대값.  
 
#define sign( x ) ((x) > 0 ? 1 : ((x) == 0 ? 0: (-1))) // 양수:1, 0, 음:-1  
 
 
두점 ( x, y ) (x2, y2)  
 
dx = x2 - x1;  
 
dy = y2 - y1;  
 
 
ix = abs(dx);  
 
iy = abs(dy);  
 
inc = max(ix, iy);  
 
 
plotx = x1;  
 
ploty = y1;  
 
x = y = 0;  
 
 
for( I = 0; I <= inc ; I++)  
 
{  
 
x += ix;  
 
y += iy;  
 
 
if( x > inc)  
 
{  
 
x -= inc;  
 
plotx += sign(dx);  
 
}  
 
if( y > inc )  
 
{  
 
y -= inc;  
 
ploty += sign(dy);  
 
}  
 
그리는 함수 ( plotx, ploty);  
 

3DMP engines

3D그래픽스 물리 수학, 프로그래밍 GPU Shader 게임엔진 알고리즘 디자인패턴 matlab etc.. 

결론적으로 이것은   인 결과 즉 사원수3개의 곱으로(또는 p 가 w=0 인 벡터여도 상관없다) 나온 결과들을 

행렬로 옮겨놓은것이다

이하내용은  지금은 알수 없는 어디선가 퍼온자료

쿼터니언은 복소수의 확장으로 i, j, k 3개의 허수와 w, x, y, z 4개의 실수부로 이루어져 있다.

즉 하나의 스칼라와 3개의 벡터 ==> 4개의 구성 요소로 이루어진 복소수를 말한다.

허수 i, j, k는 i * i = -1의 성질을 가진다.

여기서 w는 스칼라 v = (x, y, z)는 벡터이다.

쿼터니언 크기  ||q|| = Norm(q) = sqrt(w² + x^{2 +} y^{2 +} z²)

단위 쿼터니언 성질  w² + x^{2 +} y^{2 +} z² = 1

단위 쿼터니언은 (x, y, z)요소를 임의의 축, w요소를 회전 각도로 하여 4차원 공간에 표현 할 수 있다.

그림으로 간단하게 표현 하면 벡터 축 A를 중심으로 회전하는 θ는

사원수를 이용한 벡터 회전 변환

벡터 p를 사원수 q로 회전할 때 공식과 행렬에 대한 변환은 다음과 같다.

쿼트니언을 행렬 변환하는 함수

행열을 쿼트니언으로 변환하는 함수

3DMP engines

3D그래픽스 물리 수학, 프로그래밍 GPU Shader 게임엔진 알고리즘 디자인패턴 matlab etc..

[평행이동 타원과 선분의 교차판정]

* 설명 : 두 점을 지나는 직선의 y 값을 평행이동 타원의 방정식에 대입하여 x 에 대해 정리하면

           직선과 타원의 교차되는 해를 구할 수 있으며 이 범위가 선분의 범위에 포함된다면 교차된것

           올린 글중