[완료] 3차원 평면에서의 좌표 변환?
1. 현재 코딩 상황
개발 환경은 Python 2.7 - IDLE 을 사용하고 있습니다.
외부 라이브러리로 pyproj, numpy를 사용하고 있습니다.
2. 프로그램의 Input
NMEA0183 GPS 데이터 양식을 파일로 받습니다. 그 중 현재는 GPGGA값을 받으며, 이 값에는 "경,위도,고도, 위성 상태, 시간(5Hz)"을 입력받습니다.
3. 프로그램의 Output
하기 명시된 계산결과를 포함해 CSV 파일로 추출합니다.
4. 프로그램이 나타낼 데이터
GPS 데이터는 WGS84 로 구성되어 있어 계산의 편의를 위해 UTM Zone 52 평면좌표계로 변환합니다.
이 좌표계는 각 좌표의 단위가 m로 되어있는 3차원 평면좌표계이며, 거리 계산시에는 오차 보정을 위한 Cale Factor 0.996을 곱해줍니다.
이렇게 '변위'가 있고 '시간' 이 있으므로, 이 값을 통해 다음과 같은 연산을 수행하고자 합니다.
완성된 값 : 각 지점별 속도(Km/h), 가속도(G)
작업중인 값 : 'X/Y/Z축 각각의 가속도', Y축 가속도를 통해 예측 가능한 코너반경. 속도의 제곱 / 선회반경 = 횡가속도 (m/s^2)
5. 현재의 문제
3차원 평면 좌표 X,Y,Z 가 일정한 시간 간격으로 3개 이상 있을때, 이 좌표에서의 "차량 진행방향"을 X축으로 둔 상태에서 "각 축의 가속도"를 알고 싶습니다.
X/Y/Z축 각각의 가속도라고 해도 이 가속도는 '지구에서 봤을때'가 아니라 차량으로 봤을때. 가 필요하고, 이때 X축을 '차량 진행방향'으로 놓을 필요가 있습니다.
GPS 데이터는 수십만 샘플링이 있고 -1과 +1을 모두 볼 수 있으므로 '직전의 진행방향'과 '직후의 진행방향'을 알아내는것 자체는 어려운 일이 아닙니다.
그런데 이걸 벡터상에서 원점 이동 시키는거야 당연히 단순 뺄셈으로 되는데 "직전의 진행방향이 X축이 되도록" 선회시키는것이 어렵네요.
이 '회전'을 하려면 이 값들을 어떻게 처리해서 "회전할 각"을 얻어내고 어떤 연산을 통해 "회전"을 시킬수 있을까요. 도움 부탁드립니다 ㅠㅠ
질문을 제대로
질문을 제대로 이해한건지 모르겠습니다.
t1, t2의 시간에 대한 속도 v1, v2가 있을때, v2와 v1이 이루는 각도가 궁금하신게 맞나요?
알고리즘은 잘 모르겠지만, 수학적으로는, 두 벡터사이의 각도theta는 내적으로 구할수 있습니다.
v1*v2=|v1||v2|cos(theta)
가 성립합니다. 여기서 *는 내적을 나타내는 기호로 썼습니다.
그리고 벡터를 회전시키는건 3차원 회전 행렬을 곱하면 됩니다.
근데 3차원 평면은 어떤 평면인가요?
t1,t2,t3 라는 시간이
t1,t2,t3 라는 시간이 있고, 그에 따른 UTMx,y,z 값이 있습니다.
구체적으로는 아래와 같습니다.
시간은 0.0sec, 0.2sec, 0.4sec입니다. 단위는 m.
274271.9385 4142156.152 13.081
274274.3325 4142154.748 13.051
274276.7641 4142153.324 13.084
각각 X,Y,Z 값이며, 3차원 '구면'이 아닌 XYZ 좌표계가 직교하는 공간(=XYZ 절대좌표계)이라는 점에서 "평면"이라는 단어를 썼습니다.
이 값은 각 시간에 차량이 위치하는 점을 나타내고 있으며, 0.2sec 상황에서 0.0 sec와의 진행방향을 X축 기준으로 잡고 0.4sec와의 가속도X,Y,Z를 구하려고 합니다.
--
from bzImage
It's blue paper
from bzImage
It's blue paper
Z 축을 평면에
Z 축을 평면에 수직이라고 하면,
먼저 진행 방향 벡터가 X-Z 평면 상에 놓이도록 1차 회전하고,
다시 X-Z 평면을 회전하여 X 축과 진행 벡터를 일치시키면 되지 않을까요?
학부때 워낙 놀아서
학부때 워낙 놀아서 그 회전을 어떻게 해야할지 잘 모르겠습니다 ㅠ_ㅠ
--
from bzImage
It's blue paper
from bzImage
It's blue paper
원점을 일치시킨
원점을 일치시킨 상태에서 진행 방향 벡터를 (x1 y1 z1) 이라 하면
각도 arctan(y1/x1) 만큼 Z축을 중심으로 회전 시키면 진행 방향 벡터가 XZ 평면과 일치하게 되는데,
회전된 벡터를 (x2 0 z2) 라 하고
다시 각도 arctan(z2/x2) 만큼 Y 축을 중심으로 회전시킵니다.
평면 상의 회전은 벡터에 회전 행렬을 곱해 주면 되니까 간단합니다.
ps. 임의의 축을 기준으로 회전하는 루틴은 많이 나와 있습니다.
위분 말씀대로 방향 벡터와 x축의 내적으로 각을 구하고, 외적으로 두 벡터에 수직인 벡터를 구해서 이 벡터를 중심으로 구한 각 만큼 회전시키면 될 듯.
이게 맞는건지 영
이게 맞는건지 영 자신이 없네요...
아시는대로 오일러 회전은 '기준이 되는 축' 자체는 회전을 하지 않기 때문에...
XY평면 회전 후 XZ 평면을 회전시키긴 했는데 검산이 안되니...;
--
from bzImage
It's blue paper
from bzImage
It's blue paper
좌표변환
cartesian coordinate -> spherical coordinate
변환 과정을 알아보세요
회전행렬은 한국어로도 구글링하시면 많이 나올거에여
아니 그건 아닌것
아니 그건 아닌것 같아요.
직전에서야 간신히 죽을고생해서 (...) 구형좌표계의 WGS84 값을 직교좌표계의 UTM 값으로 바꿨는데
그걸 다시 구형 좌표계에서 R과 두 각을 구해도 그걸 다시 세 방향의 가속도로 분해하려면 다시 직교좌표계로 돌아와야 하고... 음...;;;;
--
from bzImage
It's blue paper
from bzImage
It's blue paper
직교좌표계에서 x축과 나란한 단위벡터는
구면좌표계에서 (1, 0, pi/2)와 같죠.
만약 직교좌표계에서 vector A를 어떻게 돌려야 x축과 나란해지는지를 알고 싶다면 구면좌표계를 쓰는 방법밖에 없습니다.
데이터로 주어진 구면좌표계 상의 좌표는 위치값이고, 지금 구면 좌표계로 나타내고자 하는 것은 속도값이니 별개죠.
(물론 구면좌표계로 나타낸 벡터 그대로 미분을 하면 cartesian coordinate로의 변환 없이 바로 '회전해야 할 각'을 알 수 있겠지만, 좌표값을 알아내기 위해 3변수 일차연립방정식을 풀어야 할 것 같군요-_- 구면좌표계 상에서는 axis vector의 미분값이 다른 axis vector에 의존하니..)
한번 프로그램을 짜봅시다.
직교좌표계로 바꾼 속도 벡터 V가 있다고 할 때,
V를 (r, phi, theta)로 나타낸다고 해 봅시다. 그러면,
r^2 = (V_x^2 + V_y^2 + V_z^2),
V_x = r sin(theta) cos(phi)
V_y = r sin(theta) sin(phi)
V_z = r cos(theta)
겠죠.
(theta는 z축과 이루는 각, phi는 V의 x-y평면으로의 projection vector와 x축이 이루는 각입니다. 그림이 없으니 이거참-_-;;;)
r = |V|,
sin theta = (V_x^2 + V_y^2) / (r^2), cos theta = V_z / r
sin phi = V_y / (r sin(theta)), cos phi = V_x / (r sin theta) // sin theta는 바로 윗줄에서 구했으므로
arcsin을 이용해 각을 바로 구하지 않고 sin, cos값을 구하는 이유는
1. 어짜피 앞으로 사용할 수식에서 sin, cos값만 이용하고,
2. 역함수는 비용이 크고 precision이 떨어지기 때문입니다.
자, 이제 (r, phi, theta)를 구했으니, 이 값을 x축과 나란하게, 즉 (r, 0, pi/2)로 만들어 봅시다..
과정은 다음과 같습니다 :
1. (r, phi, theta) -> (r, 0, theta)로 만듭니다. z축을 중심으로 -phi만큼 회전시키면 됩니다.
회전 행렬 M_1 =
[ cos (-phi), -sin (-phi) 0
sin (-phi), cos (-phi), 0
0, 0, 1]
이 되겠군요. 회전한 벡터
V' = M_1 V
로 둡시다.
2. (r, 0, theta) -> (r, 0, pi / 2)로 만듭니다. y축을 중심으로 pi/2 - theta만큼 회전시키면 되겠군요.
회전행렬 M_2 =
[ cos (pi/2 - theta), 0, -sin(pi/2 - theta)
0, 1, 0
sin(pi/2 - theta), 0, cos(pi/2 - theta)]
이 되겠군요. 회전한 벡터
V'' = M_2 (M_1 V) = M_2 M_1 V
가 됩니다.
결국, 속도 벡터 V가 있을때, V를 x축과 나란하게 만드는 회전행렬
M = M_2 M_1
이 되는군요.
만약 가속도 벡터 A를 V로 나란하게 만들고자 하면, 회전할 벡터
A' = M A
를 구하면 되겠군요.
오타 있으면 낭패..
그리고
쓰면서 글 쓰신걸 다시 한번 제대로 읽어 보았는데 -_-;
그냥 속도 단위벡터와 가속도 벡터를 외적하면 요구하시는 구심가속도가 나오지 않을까요..
V cross A = |V| |A| sin theta이니, |V| = 1로 맞추면..
A, B, C 세 좌표를
A, B, C 세 좌표를 안다고 하죠. B가 현재 위치이고 A는 바로 직전, C는 바로 직후의 위치 벡터라고 합니다.
그럼, 벡터 B-A는 현재 진행 방향이 됩니다. 물론 vA = (B-A)/dt는 직전의 속도, vB = (C-B)/dt는 직후의 속도가 되겠죠.
그리고 [vB-vA]/dt 는 가속도 벡터가 됩니다. (수치적으로 미분하는 법은 더 다양하게 있지만 그냥 간단히 적었습니다.)
이제, B-A를 X축으로 만들고 싶다는 얘기는, B-A = (a, b, c) = R로 주어진 벡터를 X=(x, 0, 0)으로 만드는 정규직교 변환을 생각하면 됩니다. 그리고 같은 변환을 다른 벡터들(속도벡터, 가속도벡터)에 모두 적용하면 되죠. 그럼 가속도 벡터에 대해서도 원하는 대로 R을 X축으로 하는 경우의 각 성분을 알아낼 수 있게 됩니다.
그럼
(F|G|H).(a,b,c)=(x,0,0)이고
F.R = x
G.R = 0
H.R = 0
F.F = 1
G.G = 1
H.H = 1
F.G = 0
F.H = 0
G.H = 0
이 되는 직교 변환을 아무거나 하나만 찾으면 됩니다.
F, G, H는 각각 열벡터이고, 어쨌든 미지수 9개에 방정식이 9개 있으므로 이 연립 방정식은 잘 풀립니다.
F, G, H를 찾아냈으면, 이제 어떤 위치벡터 R=(a,b,c)에 대해서도 그 벡터를 X축으로 만드는 행렬 (F|G|H)를 a,b,c의 함수로 찾아낼 수 있습니다. 한번만 풀면 되니까 조금 복잡해도 풀어볼만 하겠죠.
--------------------------
피할 수 있을때 즐겨라!
http://snowall.tistory.com
피할 수 있을때 즐겨라! http://melotopia.net/b
이 문제를 제대로(?)
이 문제를 제대로(?) 다루고 싶으시다면 Differential geometry에 나오는 Frame field에 대해서 이해하실 필요가 있습니다.
거기서 다루는 문제가 바로 어떤 질점이 곡선 위에서 멈추지 않고 부드럽게 움직일 때, 그 질점이 경험하는 속도, 가속도, 뒤틀림 등을 어떻게 알아내느냐는 것이거든요.
그거 이해하시면 모든 비밀이 풀릴 겁니다. 근데 좀 난해해요 -_-;
--------------------------
피할 수 있을때 즐겨라!
http://snowall.tistory.com
피할 수 있을때 즐겨라! http://melotopia.net/b
계산 완료...
결과적으로...
아기다리님 말씀대로 직전의 진행방향의 각을 얻어내고, 계산하고자 하는 벡터를 구면좌표계로 변환한 뒤, 구면좌표계에서 각을 빼버리고 구하는 방식으로 원하는 값을 얻었습니다.
... 계산은 했는데
문제는 데이터 소스가 GPS로, HDOP가 1m 이상인데다 실제 값은 그보다 큰 오차를 보여줘, 거의 쓰레기값에 가까운 -_-);;; 결과만 잔뜩이네요.
코너 한가운데에서 반대방향으로 가속도가 걸린다거나.. 음.;
아기다리님께 감사드립니다.
--
from bzImage
It's blue paper
from bzImage
It's blue paper
코너 한가운데에서
코너 한가운데에서 반대방향이라 함은 무엇에 대해서 반대방향인지요?
최초 진행방향의 반대 방향으로 가속도가 존재한다면, 코너를 돌기위한 구심가속도로 보입니다.
'코너를 돌기 위한
'코너를 돌기 위한 구심가속도'와 '반대방향의 가속도'가 프레임 중간중간마다 발견이 된다는 이야기지말입니다.
--
from bzImage
It's blue paper
from bzImage
It's blue paper
이미 해결이 되신 듯
이미 해결이 되신 듯 하지만,
단순히 가속도의 진행방향 성분과 그와 수직인 성분을 구하시는 것이 목적이라면 궂이
복잡하게 좌표계의 회전까지 계산하실 필요가 없습니다.
가령, t1,t2,t3 에서의 속도가 (vx,vy,vz), 가속도가(ax,ay,az) 라면, 가속도의
속도방향 성분은, 속도방향으의 단위벡터 (vx,vy,vz)/sqrt(vx^2 + vy^2 + vz^2)과
가속도 벡터의 내적을 구하면 바로 나옵니다. 즉,
a_vert = (ax*vx + ay*vy + az*vz) / sqrt(vx^2 + vy^2 + vz^2)이 되겠죠.
속도와 수직방향 가속도 성분은 피타고라스 정리를 사용하면 쉽게 구할 수 있구요.
a_para = sqrt(ax^2 + ay^2 + az^2 - a_vert^2).
혹은, 아기다리님께서 적었듯이 속도방향의 단위벡터와 가속도벡터의 외적을 구해도 되구요.
직관적으로 "그렇게
직관적으로 "그렇게 되나...?" 란 생각이 들긴 하는데 일단 월요일날 한번 그 로직을 작성해보고 다시 말씀드려보겠습니다 (_ _)
--
from bzImage
It's blue paper
from bzImage
It's blue paper
상당히 유효한 범위
상당히 유효한 범위 내에서 값을 얻을 수 있네요.
둘의 비교는 따로 해봐야게지만... 감사합니다 (_ _)
P.s. 사실 지금까지 경위도 입력을 잘못받아서 노이즈가 폭풍처럼 쩔어주고 계셨습니다 orz
--
from bzImage
It's blue paper
from bzImage
It's blue paper
댓글 달기