이미지의 affine 변형을 계산. 입력 프레임 img 내에서 p파라메터에 있는 affine계수에 따라 일정 부분을 취해 sz(예를 들면, 32x32)크기의 영상으로 만들어 리턴. p 행렬 내에 복수개의 affine 계수가 넘어가면 복수 개의 sz 크기 영상들이 생성되어 리턴.
function wimg = warpimg(img, p, sz)
% function wimg = warpimg(img, p, sz)
%
% img(h,w) : 입력 프레임
% p(6,n) : Affine 인자 행렬(n은 affine 인자 set의 갯수)
% sz(th,tw) ： 입력 영상 내에서 취할 template 크기
%
% 이미지의 affine변형을 고속으로 처리하기 위한 함수
%% Copyright (C) 2005 Jongwoo Lim and David Ross.
%% All rights reserved.
if (nargin < 3)
sz = size(img);
end
if (size(p,1) == 1) % size(param0,1)=1
p = p(:); % param0는 행벡터인데 열벡터로 만들어 줌
end
w = sz(2); % 32
h = sz(1); % 32
n = size(p,2); % size(p,2)=1
[x,y] = meshgrid([1:w]-w/2, [1:h]-h/2);
% x=[ -15 ... 16;
% -15 ... 16;
% ...
% -15 ... 16], size(x)=32x32
% y=[ -15 .... -15;
% -14 ... -14;
% ...
% 16 ... 16]; size(y)=32x32
pos = reshape(cat(2, ones(h*w,1),x(:),y(:)) ...
* [p(1,:) p(2,:); p(3:4,:) p(5:6,:)], [h,w,n,2]);
% size(ones(h*w,1))=1024, size(x(:))=1024
% CAT Concatenate arrays.
% CAT(DIM,A,B) concatenates the arrays A and B along
% the dimension DIM.
% CAT(2,A,B) is the same as [A,B].
% CAT(1,A,B) is the same as [A;B].
% cat(2, ones(h*w,1),x(:),y(:)) =
% [ 1 -15 -15; 1 -15 -14; 1 -15 -14; ...; 1 16 15; 1 16 16]: 1024x3
% [p(1,:) p(2,:); p(3:4,:) p(5:6,:)]
% 177.0000 147.0000
% 0.1166 -0.6368
% -0.6368 3.4771
% p1=cat(2, ones(h*w,1),x(:),y(:))*[p(1,:) p(2,:); p(3:4,:) p(5:6,:)]
% size(p1)=1024x2
% size(pos)=32x32x1x2
%
% (Example) For one pixel and one set of affine params:
% (1 x y).(tx ty; a11 a21; a12 a22)=(tx+a11.x+a12.y ty+a21.x+a22.y)
% For all pixels in sz (32x32=1024 pixels) and one set of affine parameters:
% [1 x1 y1; 1 x2 y2; ....].[tx ty; a11 a21; a12 a22]
% = [tx+a11.x+a12.y ty+a21.x+a22.y; tx+a11.x+a12.y ty+a21.x+a22.y; ...]
% *size: (1024x3).(3x2) = (1024x2)
% For all pixels & n sets of affine parameters:
% [1 x1 y1; 1 x2 y2; ....].[tx(1) tx(2) ... tx(n) ty(1) ty(2) ... ty(n);
% a11(1) a11(2) ...a11(n) a21(1) a21(2) ...a21(n);
% a12(1) a12(2) ...a12(n) a22(1) a22(2) ...a22(n)];
% * size: (1024x3).(3xnx2)=(1024xnx2) --> (32x32xnx2)
wimg = squeeze(interp2(img, pos(:,:,:,1), pos(:,:,:,2)));
% size(wimg)=32x32
% interp2(z,x1,y1): bilinear interpolation for the values of the
% underlying 2-D function z at points in matrices x1 and y1.
%%pos(:,:,:,1) : x麟깃앤黎
%%pos(:,:,:,2) ： y麟깃앤黎
wimg(find(isnan(wimg))) = 0; %%nan인 값을 찾아 0으로 만들어 줌
%ISNAN True for Not-a-Number.
% ISNAN(X) returns an array that contains 1's where
% the elements of X are NaN's and 0's where they are not.
% For example, ISNAN([pi NaN Inf -Inf]) is [0 1 0 0].
% B = SQUEEZE(A) returns an array B with the same elements as
% A but with all the singleton dimensions removed. A singleton
% is a dimension such that size(A,dim)==1;
2013년 11월 27일 수요일
2013년 11월 22일 금요일
Particle filter에 의한 영역 추적
연속 이미지에서 미리 설정해준 물체 영역(patch)을 particle filters에 의해 추적한다. 이 때, 파티클의 발생 방법에 대해 알아 본다.
먼저 사용자가 영상 내에서 물체 영역을 설정한다. 구체적으로 예를 들어 보면,
-red window 중심 위치, 크기: (177(col_center) 147(row_center)), [115(col_length)x 145(row_length)]
(1) 처음 사용자가 설정해 준 영역은 추적 대상이 된다. 현재 설정 영역의 크기는 너무 크므로 가로,세로 폭을 줄여 32x32크기로 만든다. 이를 templates로 한다.
(2) 목표 물체는 3차원 운동을 할 수 있지만 이미지에 맺히는 2D 운동은 affine motion이라 가정한다. 물체의 2d운동을 표현하기 위해 geometric affine parameter를 설정한다.
% 추적 사각영역의 상태변수:
먼저 사용자가 영상 내에서 물체 영역을 설정한다. 구체적으로 예를 들어 보면,
(figure from [2])
-image size: 352(가로)x288(세로)-red window 중심 위치, 크기: (177(col_center) 147(row_center)), [115(col_length)x 145(row_length)]
(1) 처음 사용자가 설정해 준 영역은 추적 대상이 된다. 현재 설정 영역의 크기는 너무 크므로 가로,세로 폭을 줄여 32x32크기로 만든다. 이를 templates로 한다.
(2) 목표 물체는 3차원 운동을 할 수 있지만 이미지에 맺히는 2D 운동은 affine motion이라 가정한다. 물체의 2d운동을 표현하기 위해 geometric affine parameter를 설정한다.
% 추적 사각영역의 상태변수:
% col-center, row-center, 수직 스케일, 회전각 theta, 가로세로비, 각도 파이
% state=(tx, ty, s1, theta, s2/s1, phi)
geo_param = [177, 147, 115/32, 0, 145/115, 0]; % geometric affine 인자
% 가중치를 가진 particle들을 발생시켜 보자
n = 600; % 샘플의 갯수는 600개로 한다.
N = 1024; % 32x32 템플릿을 1차원 벡터로 표현하면 32x32=1024이다.
params = repmat(geo_param, [1,n]); % geo affine param을 600개 복사. 6x600행렬 생성
randMatrix = randn(6,n); % 6x600개 정규 분포의 행렬 생성,
% 정규 분포이니 평균은 0이다: mean(randMatrix,2)=(0.023 0.042 ...-0.036)
% affsig = [4 4 0.02 0.02 0.005 0.001]: 즉, center 위치 변화는 크다.
% 나머지 인자는 약간 변화한다 가정
params = params + randMatrix.*repmat(affsig(:),[1,n]); % 600개의 parameters생성
% 입력 영상을 600개의 affine 인자에 따라 변형 시킴. 즉, 입력 frame을 변형 시킨다.
% 출력은 32x32 크기의 템플릿 600개
wimgs = warpimg(frame, affparam2mat(params), [32,32]);
% template에서 600개 각각을 빼서 잔차 이미지를 구한다.
% template이미지는 기준영상으로 첫번째영상에서 물체
% 이미지, PCA영상, 모델 영상 등이다.
diff = repmat(template(:),[1,n]) - reshape(wimgs,[N,n]);
% 잔차 값을 이용하여 confidence를 계산
confidence = exp(-sum(diff.^2)./condenssig)'; %condenssig=0.25(잔차 분산 값)
% confidence를 확률 분포(pdf)로 만들기 위해 정규화한다.
confidence = confidence ./ sum(confidence);
% 확률 값이 가장 큰 후보를 선택
[maxprob, maxidx] = max(confidence); % maxprob=0.2233(130번째)
% 이 샘플에 대한 affine para 값을 저장
affine_est = affparam2mat(params(:,maxidx)); %[a11 a12 tx; a21 a22 ty]
% 최적 이미지도 저장한다
timgs = [timgs wimgs(:,:,maxidx)];
참고 문헌
[1] D.A. Ross, J. Lim, et al., Incremental Learning for Robust Visual Tracking
[2] Dong Wang et. al., Online Object Tracking with sparse prototypes
% state=(tx, ty, s1, theta, s2/s1, phi)
geo_param = [177, 147, 115/32, 0, 145/115, 0]; % geometric affine 인자
% 가중치를 가진 particle들을 발생시켜 보자
n = 600; % 샘플의 갯수는 600개로 한다.
N = 1024; % 32x32 템플릿을 1차원 벡터로 표현하면 32x32=1024이다.
params = repmat(geo_param, [1,n]); % geo affine param을 600개 복사. 6x600행렬 생성
randMatrix = randn(6,n); % 6x600개 정규 분포의 행렬 생성,
% 정규 분포이니 평균은 0이다: mean(randMatrix,2)=(0.023 0.042 ...-0.036)
% affsig = [4 4 0.02 0.02 0.005 0.001]: 즉, center 위치 변화는 크다.
% 나머지 인자는 약간 변화한다 가정
params = params + randMatrix.*repmat(affsig(:),[1,n]); % 600개의 parameters생성
% 입력 영상을 600개의 affine 인자에 따라 변형 시킴. 즉, 입력 frame을 변형 시킨다.
% 출력은 32x32 크기의 템플릿 600개
wimgs = warpimg(frame, affparam2mat(params), [32,32]);
% template에서 600개 각각을 빼서 잔차 이미지를 구한다.
% template이미지는 기준영상으로 첫번째영상에서 물체
% 이미지, PCA영상, 모델 영상 등이다.
diff = repmat(template(:),[1,n]) - reshape(wimgs,[N,n]);
% 잔차 값을 이용하여 confidence를 계산
confidence = exp(-sum(diff.^2)./condenssig)'; %condenssig=0.25(잔차 분산 값)
% confidence를 확률 분포(pdf)로 만들기 위해 정규화한다.
confidence = confidence ./ sum(confidence);
% 확률 값이 가장 큰 후보를 선택
[maxprob, maxidx] = max(confidence); % maxprob=0.2233(130번째)
% 이 샘플에 대한 affine para 값을 저장
affine_est = affparam2mat(params(:,maxidx)); %[a11 a12 tx; a21 a22 ty]
% 최적 이미지도 저장한다
timgs = [timgs wimgs(:,:,maxidx)];
참고 문헌
[1] D.A. Ross, J. Lim, et al., Incremental Learning for Robust Visual Tracking
[2] Dong Wang et. al., Online Object Tracking with sparse prototypes
2013년 11월 19일 화요일
Affine transformation
Affine 변환은 위치 이동과 선형변환(linear transformation)의 결합이다.
(x' y' 1)' = [a11 a12 tx; a21 a22 ty; 0 0 1].(x y 1)'
로 표현된다. 행렬로 더 축약해서 표현하면
x' = H.x = [A t; 0' 1].x
이다. 여기서 행렬 A는 2x2 크기의 비 특이(non-singular) 행렬이다. 평면상에서 일어나는 affine 변환은 6자유도를 가진다 (행렬 내에 6개의 미지수가 있다). 만일 3개의 대응점이 있다면 하나의 대응점이 2개의 식을 주므로 affine변환의 6개 인자 결정이 가능하다.
선형 변환 행렬 A는 2개의 기하 변환으로 설명할 수 있다.
(1) 회전
(2) 비 등방성(non-isotropic) 크기변화
아래 그림은 강체 회전과 shearing 변형의 2가지 기하 변환을 보여준다.
기하 변환의 값들은 2x2 크기의 A행렬을 특이값 분해하면 얻을 수 있다:
(예1) 파이 각 적용 여부에 따른 차이 점
>> clear all
>> fi=30*pi/180; % phi을 30도로
>> th=60*pi/180; % theta를 60도로
>>
>> Rf=[cos(fi) -sin(fi); sin(fi) cos(fi)];
>> Rt=[cos(th) -sin(th); sin(th) cos(th)];
>> D=[1.5 0; 0 1.0];
>>
>> A=Rt*Rf'*D*Rf
>> [u,d,v]=svd(A);
>> % d는 D이다. Rt==uv', Rf=v'이다.
>>
>> % 실제 변형을 테스트 해보자.
% (아래 그림에서 원형점(정사각형))
>> pts=[1 1 -1 -1; 1 -1 -1 1];
>> pts =
1 1 -1 -1
1 -1 -1 1
>> D*pts =
1.5000 1.5000 -1.5000 -1.5000
1.0000 -1.0000 -1.0000 1.0000
% shearing 변형 없이 스케일링만 일어남 (아래 그림에서 삼각점)
>> Rf'*D*Rf*pts
1.1585 1.5915 -1.1585 -1.5915
0.9085 -1.3415 -0.9085 1.3415
% shearing 변형이 일어남 (아래 그림에서 사각점)
(예2) tracking에 적용하기 위한 geometric parameter로 변환.
영상에서 물체를 추적할 때 보통 물체 영역을 둘러싸는 작은 사각의 윈도(window patch)를 가정한다.
물체가 크기 변화, 회전, 이동 할 때, 이에 따라 윈도 영역도 크기변화, 회전, 이동해야 한다. 이러한 기하학적 변화를 나타내는 parameter를 affine 변환에서 유도한다.
즉, 일반적으로 알고 있는 affine 변환의 6개의 parameter를 geometric affine 인자로 변환하는 matlab코드는 아래와 같다. Geometric parameter 6개는
(x, y, scale, th, aspect, skew) 이다.
여기서, (x,y)는 윈도의 중심 위치, scale은 수직방향 스케일, th는 회전 각, aspect는 수직-수평 비, skew는 파이 각을 나타낸다.
이러한 변환을 통해 추적되는 물체의 상태 변수(state variable)는 물체를 둘러싸는 윈도의 상태를 나타내는 6개의 geometric parameter가 된다.
function q = affparam2geom(p)
% function q = affparam2geom(p)
%
% 입력 ： p，2x3의 affine 계수；
% 출력 ： q，geometric affine 계수로 각도 스케일비 등을 나타냄；
% The functions affparam2geom and affparam2mat convert a 'geometric'
% affine parameter to/from a matrix form (2x3 matrix).
%
% affparam2geom converts a 2x3 matrix to 6 affine parameters
% (x, y, scale, th, aspect, skew), and affparam2mat does the inverse.
%
% p(6) : [p(1) p(3) p(4); p(2) p(5) p(6)]
%
% x' p(3) p(4) p(1) x
% y' = p(5) p(6) p(2) * y
% 1 0 0 1 1
%
% p(3) p(4)
% A =
% p(5) p(6)
%
% q(6) : [dx dy sc th sr phi]
% dx,dy : left,top의 원점을 기준으로 box의 중심 위치 표현
% sc,sr : ramda1(수직방향으로 스케일), ramda2/ramda1 (길이 비율)
% 수직 방향 스케일은 box크기가 32이므로 32의 배수로 주어짐
% th ： 각도 theta
% phi ： 각도 phi
%
% Reference "Multiple View Geometry in Computer Vision" by Richard
% Hartley and Andrew Zisserman.
% 출력은 [tx, ty, ramda1, theta, ramda2/ramda1, phi]
% MVG 책의 39~40페이지에 설명
% Copyright (C) Jongwoo Lim and David Ross. All rights reserved.
A = [ p(3), p(4); p(5), p(6) ]; % A=[3.594 0; 0 4.531]: y(col)크기가 더 큼
%% A = USV' = (UV')(VSV') = R(th)R(-phi)SR(phi)
[U,S,V] = svd(A);
% (Ex.) 축의 순서 결정. For example,
% S=[4.53 0; 0 3.59]; U=[0 1; 1 0];
% 영상내 설정은 x축(col), y축(row) 순서임. Row가 더 크니
% U가 (0 1)', (1 0)' 벡터의 순서로 나옴(즉, y축이 먼저 나옴)
% 따라서, x축이 우선되게 축 순서 조정 필요
if (det(U) < 0) % 항상 x축이 우선되게 하기 위함
U = U(:,2:-1:1); V = V(:,2:-1:1); % 2열을 1열로
S = S(2:-1:1,2:-1:1); % 2행을 1행으로, 2열을 1열로
% U=[1 0; 0 1], S=[3.59 0; 0 4.53]로 수정 됨
end
%%*******중심 위치*******%%
q(1) = p(1);
q(2) = p(2);
%%*******(col, row)*******%%
% [cos(th) -sin(th); sin(th) cos(th)]=uv'
% =[u11 u12; u21 u22][v11 v21; v12 v22]
% =[u11.v11+u12.v21 ...; u21.v11+u22.v12 ...]-->
q(4) = atan2(U(2,1)*V(1,1)+U(2,2)*V(1,2), U(1,1)*V(1,1)+U(1,2)*V(1,2)); % theta
phi = atan2(V(1,2),V(1,1)); % angle phi
if (phi <= -pi/2) % phi가 -90도 이하면 그냥 -90도로
c = cos(-pi/2); s = sin(-pi/2);
R = [c -s; s c]; V = V * R; S = R'*S*R;
end
if (phi >= pi/2) % phi가 90도 이상이면 그냥 90도로 둠
c = cos(pi/2); s = sin(pi/2);
R = [c -s; s c]; V = V * R; S = R'*S*R;
end
%%*******크기비율*******%%
q(3) = S(1,1); % ramda1
q(5) = S(2,2)/S(1,1); % ramda2/ramda1
%%*******각도파이*******%%
q(6) = atan2(V(1,2),V(1,1)); % phi
q = reshape(q, size(p));
% 출력은 [tx, ty, ramda1, theta, ramda2/ramda1, phi]
(예3) Geometric affine parameters를 2x3크기의 affine 행렬로 변환하는 코드
function q = affparam2mat(p)
% function q = affparam2mat(p)
%
% The functions affparam2geom and affparam2mat convert a 'geometric'
% affine parameter to/from a matrix form (2x3 matrix).
%
% 입력 ： q，기하학적 affine 계수；
% 출력 ： p，2x3행렬에서의 affine 계수；
%
% affparam2geom converts a 2x3 matrix to 6 affine parameters
% (x, y, th, scale, aspect, skew), and affparam2mat does the inverse.
%
% p(6,n) : [dx dy sc th sr phi]'
% q(6,n) : [q(1) q(3) q(4); q(2) q(5) q(6)]
%
% Reference "Multiple View Geometry in Computer Vision" by Richard
% Hartley and Andrew Zisserman.
%
% Copyright (C) Jongwoo Lim and David Ross. All rights reserved.
sz = size(p); % size(p)=1x6
if (length(p(:)) == 6) % p(:)= (177.0000; 147.0000; 3.5938; 0; 1.2609; 0)
p = p(:); % length(p(:))=6
% p = (177.0000; 147.0000; 3.5938; 0; 1.2609; 0)
end
s = p(3,:); th = p(4,:); r = p(5,:); phi = p(6,:);
cth = cos(th); sth = sin(th); cph = cos(phi); sph = sin(phi);
ccc = cth.*cph.*cph; ccs = cth.*cph.*sph; css = cth.*sph.*sph;
scc = sth.*cph.*cph; scs = sth.*cph.*sph; sss = sth.*sph.*sph;
q(1,:) = p(1,:); q(2,:) = p(2,:);
q(3,:) = s.*(ccc +scs +r.*(css -scs)); q(4,:) = s.*(r.*(ccs -scc) -ccs -sss);
q(5,:) = s.*(scc -ccs +r.*(ccs +sss)); q(6,:) = s.*(r.*(ccc +scs) -scs +css);
q = reshape(q, sz);
참고 문헌
[1] Hartley and Zisserman, Multiple View Geometry, Cambridge press
(x' y' 1)' = [a11 a12 tx; a21 a22 ty; 0 0 1].(x y 1)'
로 표현된다. 행렬로 더 축약해서 표현하면
x' = H.x = [A t; 0' 1].x
이다. 여기서 행렬 A는 2x2 크기의 비 특이(non-singular) 행렬이다. 평면상에서 일어나는 affine 변환은 6자유도를 가진다 (행렬 내에 6개의 미지수가 있다). 만일 3개의 대응점이 있다면 하나의 대응점이 2개의 식을 주므로 affine변환의 6개 인자 결정이 가능하다.
선형 변환 행렬 A는 2개의 기하 변환으로 설명할 수 있다.
(1) 회전
(2) 비 등방성(non-isotropic) 크기변화
아래 그림은 강체 회전과 shearing 변형의 2가지 기하 변환을 보여준다.
기하 변환의 값들은 2x2 크기의 A행렬을 특이값 분해하면 얻을 수 있다:
(예1) 파이 각 적용 여부에 따른 차이 점
>> clear all
>> fi=30*pi/180; % phi을 30도로
>> th=60*pi/180; % theta를 60도로
>>
>> Rf=[cos(fi) -sin(fi); sin(fi) cos(fi)];
>> Rt=[cos(th) -sin(th); sin(th) cos(th)];
>> D=[1.5 0; 0 1.0];
>>
>> A=Rt*Rf'*D*Rf
>> [u,d,v]=svd(A);
>> % d는 D이다. Rt==uv', Rf=v'이다.
>>
>> % 실제 변형을 테스트 해보자.
% (아래 그림에서 원형점(정사각형))
>> pts=[1 1 -1 -1; 1 -1 -1 1];
>> pts =
1 1 -1 -1
1 -1 -1 1
>> D*pts =
1.5000 1.5000 -1.5000 -1.5000
1.0000 -1.0000 -1.0000 1.0000
% shearing 변형 없이 스케일링만 일어남 (아래 그림에서 삼각점)
>> Rf'*D*Rf*pts
1.1585 1.5915 -1.1585 -1.5915
0.9085 -1.3415 -0.9085 1.3415
% shearing 변형이 일어남 (아래 그림에서 사각점)
(예2) tracking에 적용하기 위한 geometric parameter로 변환.
영상에서 물체를 추적할 때 보통 물체 영역을 둘러싸는 작은 사각의 윈도(window patch)를 가정한다.
물체가 크기 변화, 회전, 이동 할 때, 이에 따라 윈도 영역도 크기변화, 회전, 이동해야 한다. 이러한 기하학적 변화를 나타내는 parameter를 affine 변환에서 유도한다.
즉, 일반적으로 알고 있는 affine 변환의 6개의 parameter를 geometric affine 인자로 변환하는 matlab코드는 아래와 같다. Geometric parameter 6개는
(x, y, scale, th, aspect, skew) 이다.
여기서, (x,y)는 윈도의 중심 위치, scale은 수직방향 스케일, th는 회전 각, aspect는 수직-수평 비, skew는 파이 각을 나타낸다.
이러한 변환을 통해 추적되는 물체의 상태 변수(state variable)는 물체를 둘러싸는 윈도의 상태를 나타내는 6개의 geometric parameter가 된다.
function q = affparam2geom(p)
% function q = affparam2geom(p)
%
% 입력 ： p，2x3의 affine 계수；
% 출력 ： q，geometric affine 계수로 각도 스케일비 등을 나타냄；
% The functions affparam2geom and affparam2mat convert a 'geometric'
% affine parameter to/from a matrix form (2x3 matrix).
%
% affparam2geom converts a 2x3 matrix to 6 affine parameters
% (x, y, scale, th, aspect, skew), and affparam2mat does the inverse.
%
% p(6) : [p(1) p(3) p(4); p(2) p(5) p(6)]
%
% x' p(3) p(4) p(1) x
% y' = p(5) p(6) p(2) * y
% 1 0 0 1 1
%
% p(3) p(4)
% A =
% p(5) p(6)
%
% q(6) : [dx dy sc th sr phi]
% dx,dy : left,top의 원점을 기준으로 box의 중심 위치 표현
% sc,sr : ramda1(수직방향으로 스케일), ramda2/ramda1 (길이 비율)
% 수직 방향 스케일은 box크기가 32이므로 32의 배수로 주어짐
% th ： 각도 theta
% phi ： 각도 phi
%
% Reference "Multiple View Geometry in Computer Vision" by Richard
% Hartley and Andrew Zisserman.
% 출력은 [tx, ty, ramda1, theta, ramda2/ramda1, phi]
% MVG 책의 39~40페이지에 설명
% Copyright (C) Jongwoo Lim and David Ross. All rights reserved.
A = [ p(3), p(4); p(5), p(6) ]; % A=[3.594 0; 0 4.531]: y(col)크기가 더 큼
%% A = USV' = (UV')(VSV') = R(th)R(-phi)SR(phi)
[U,S,V] = svd(A);
% (Ex.) 축의 순서 결정. For example,
% S=[4.53 0; 0 3.59]; U=[0 1; 1 0];
% 영상내 설정은 x축(col), y축(row) 순서임. Row가 더 크니
% U가 (0 1)', (1 0)' 벡터의 순서로 나옴(즉, y축이 먼저 나옴)
% 따라서, x축이 우선되게 축 순서 조정 필요
if (det(U) < 0) % 항상 x축이 우선되게 하기 위함
U = U(:,2:-1:1); V = V(:,2:-1:1); % 2열을 1열로
S = S(2:-1:1,2:-1:1); % 2행을 1행으로, 2열을 1열로
% U=[1 0; 0 1], S=[3.59 0; 0 4.53]로 수정 됨
end
%%*******중심 위치*******%%
q(1) = p(1);
q(2) = p(2);
%%*******(col, row)*******%%
% [cos(th) -sin(th); sin(th) cos(th)]=uv'
% =[u11 u12; u21 u22][v11 v21; v12 v22]
% =[u11.v11+u12.v21 ...; u21.v11+u22.v12 ...]-->
q(4) = atan2(U(2,1)*V(1,1)+U(2,2)*V(1,2), U(1,1)*V(1,1)+U(1,2)*V(1,2)); % theta
phi = atan2(V(1,2),V(1,1)); % angle phi
if (phi <= -pi/2) % phi가 -90도 이하면 그냥 -90도로
c = cos(-pi/2); s = sin(-pi/2);
R = [c -s; s c]; V = V * R; S = R'*S*R;
end
if (phi >= pi/2) % phi가 90도 이상이면 그냥 90도로 둠
c = cos(pi/2); s = sin(pi/2);
R = [c -s; s c]; V = V * R; S = R'*S*R;
end
%%*******크기비율*******%%
q(3) = S(1,1); % ramda1
q(5) = S(2,2)/S(1,1); % ramda2/ramda1
%%*******각도파이*******%%
q(6) = atan2(V(1,2),V(1,1)); % phi
q = reshape(q, size(p));
% 출력은 [tx, ty, ramda1, theta, ramda2/ramda1, phi]
(예3) Geometric affine parameters를 2x3크기의 affine 행렬로 변환하는 코드
function q = affparam2mat(p)
% function q = affparam2mat(p)
%
% The functions affparam2geom and affparam2mat convert a 'geometric'
% affine parameter to/from a matrix form (2x3 matrix).
%
% 입력 ： q，기하학적 affine 계수；
% 출력 ： p，2x3행렬에서의 affine 계수；
%
% affparam2geom converts a 2x3 matrix to 6 affine parameters
% (x, y, th, scale, aspect, skew), and affparam2mat does the inverse.
%
% p(6,n) : [dx dy sc th sr phi]'
% q(6,n) : [q(1) q(3) q(4); q(2) q(5) q(6)]
%
% Reference "Multiple View Geometry in Computer Vision" by Richard
% Hartley and Andrew Zisserman.
%
% Copyright (C) Jongwoo Lim and David Ross. All rights reserved.
sz = size(p); % size(p)=1x6
if (length(p(:)) == 6) % p(:)= (177.0000; 147.0000; 3.5938; 0; 1.2609; 0)
p = p(:); % length(p(:))=6
% p = (177.0000; 147.0000; 3.5938; 0; 1.2609; 0)
end
s = p(3,:); th = p(4,:); r = p(5,:); phi = p(6,:);
cth = cos(th); sth = sin(th); cph = cos(phi); sph = sin(phi);
ccc = cth.*cph.*cph; ccs = cth.*cph.*sph; css = cth.*sph.*sph;
scc = sth.*cph.*cph; scs = sth.*cph.*sph; sss = sth.*sph.*sph;
q(1,:) = p(1,:); q(2,:) = p(2,:);
q(3,:) = s.*(ccc +scs +r.*(css -scs)); q(4,:) = s.*(r.*(ccs -scc) -ccs -sss);
q(5,:) = s.*(scc -ccs +r.*(ccs +sss)); q(6,:) = s.*(r.*(ccc +scs) -scs +css);
q = reshape(q, sz);
참고 문헌
[1] Hartley and Zisserman, Multiple View Geometry, Cambridge press
2013년 11월 9일 토요일
차원 축소(dimensionality reduction)
크기가 32x32인 얼굴 영상이 50장 있다. 한장의 영상은 픽셀이 1024개가 있는 하나의 벡터이다.
각각의 벡터는 1024차원에서 하나의 점이다. 50개가 있으므로 1024차원의 공간에서 50개의 점을 정의한다.
행렬 c는 [5x50]으로, 50장의 영상을 각기 5개의 계수치로 표현한 것이다.
mean(c,2)를 통해 5개 계수치의 평균을 구해서 m2에 저장한다.
어떤 새로운 이미지 t1를 템플릿과 비교할 때는 c2=u1'*t1으로 5개의 계수를 구하고 norm(m2-c2)로 거리를 계산하면 된다.
따라서 1024개의 벡터가 50개 이다. 즉, 행렬 A의 크기는 1024x50이다.
각각의 벡터는 1024차원에서 하나의 점이다. 50개가 있으므로 1024차원의 공간에서 50개의 점을 정의한다.
점 하나에 1024개의 좌표가 필요하다. 점 표현을 위해 사용되는 주축이 1024개이기 때문이다.
데이터를 압축해서 표현하기 위해 주축의 수를 줄여본다.
50개 점의 분포, 즉 분산이 가장 큰 축을 몇 개 찾으면 된다. 분산이 큰 방향이 데이터의 변량을 잘 표현하는 방향이다.
5개을 찾아 보자.
데이터는 크기가 1024x50인 행렬이다. 먼저 평균을 구하자. 즉, 데이터의 중심점을 지나는 주축을 구할 것이다.
mu = mean(data,2);
각 데이터에서 평균을 제거하자.
data = data-repmat(mu, [1,50]);
data를 특이값 분해한다.
[u,d,v] = svd(data); % data=udv'
d 행렬은 특이값이 들어 있는 대각 행렬이다.
u 행렬의 1~5열을 취한다.
u1=u(:,1:5);
주축 5개를 구했으므로 첫번째 데이터를 새 축에서 표시해 보자.
c1=u1'*data(:,1); % 5x1=(5x1024).(1024x1)
5개의 값이 나오는데 이는 5개의 주축에서 좌표 값이다.
즉, 데이터의 표현이 1024개에서 5개로 줄었다.
이제 5개의 좌표 값 만으로 원 데이터를 재생해보자.
o1=c1(1)*u1(:,1)+c1(2)*u1(:,2)+...+c1(5)*u1(:,5)
이다. matlab 명령으로
o1 = u1*c1; % 1024x1=(1024x5).(5x1). o1=u1*c1=u1*u1'*data(:,1)
이다. 원 데이터와 재생 데이터의 차이는
res1=data(:,1)-o1;
이다.
모든 데이터를 새 축에서 표시해 보자.
c = u1'*data; % (5x50) = (5x1024).(1024x50)
o = u1*c; % (1024x50) = (1024x5).(5x50)
res = data-o; % (1024x50)-(1024x50)
모든 데이터를 새 축에서 표시해 보자.
c = u1'*data; % (5x50) = (5x1024).(1024x50)
o = u1*c; % (1024x50) = (1024x5).(5x50)
res = data-o; % (1024x50)-(1024x50)
행렬 c는 [5x50]으로, 50장의 영상을 각기 5개의 계수치로 표현한 것이다.
mean(c,2)를 통해 5개 계수치의 평균을 구해서 m2에 저장한다.
어떤 새로운 이미지 t1를 템플릿과 비교할 때는 c2=u1'*t1으로 5개의 계수를 구하고 norm(m2-c2)로 거리를 계산하면 된다.
Sequential KL 변환
영상 내에 있는 특정 물체(예를 들면 얼굴)를 추적 한다고 하자. 이 때, 물체를 표현하기 위해 여러 개의 영상 템플릿을 사용하는 PCA(principal component analysis)을 적용할 수 있다.
물체의 외관은 추적 동안 여러가지 원인으로 크게 변한다. PCA 기법을 통해 물체를 추적할 때, 물체를 표현하는 주축(principal axis)을 효과적으로 갱신하는 것이 필요하다.
연속 영상에서 물체 추적동안 새로 얻어지는 프레임들을 고려하여 주축을 갱신할 필요가 있을 때 sequential KL 알고리즘을 사용할 수 있다. (Sequential PCA 또는 Sequential visual tracking으로 부를 수도 있다.)
skl 변환은 sequential kl 변환이다. 연속 신호에서 데이터가 차례로 도착한 것을 처리하거나, 큰 행렬의 SVD를 계산해야 할 때, 큰 행렬을 작은 행렬로 분할하여 여러번 반복 처리할 때 사용한다.
예를 들면, A'=[A E] 행렬을 생각해 보자. A행렬은 10장의 영상 템플릿으로 이루어진 행렬이다. E행렬은 새로 매칭된 5장의 프레임에서 추출된 5장의 영상 템플릿으로 이루어진 행렬이다.
이 때, 이용할 수 있는 템플릿이 모두 15개가 되었으므로, 이것 모두가 포함된 A' 행렬에 대해 새로 주축을 구하는 것이 필요하다. 그런데, A행렬에 대해서는 특이값 분해를 통해 이미 구해 놓은 주축이 있다고 가정하자.
효율적 계산을 위해서는 A'행렬을 다시 특이값 분해를 하는 일 없이, 새롭게 추가된 E행렬에 대해서만 특이값 분해를 하고, 그 결과를 A행렬에 대한 특이값 분해의 결과와 결합한다. 이렇게 결합된 결과가 A'행렬에 대한 특이값 분해의 값과 동일하도록 만드는 일이 SKL이다.
R-SVD[2]는 A의 분해에 기반하여 효과적으로 행렬 A'의 분해 A' = [A|E] = U'D'V''를 계산할 수 있다. [U,D,V]=svd(A)이라고 가정하자. 이 때, A'의 특이값 분해는
과 같이 구해진다. 이 알고리즘을 해석해 보자.
MxL크기의 행렬 B가 있다. 특이값 분해하면
B = UDV'
이다. U는 왼쪽 특이 벡터들로 구성된 MxL크기의 직교(orthonormal) 행렬이다. D는 LxL 크기의 행렬이다. 대각요소가 B행렬의 특이 값이다.
V는 각 열이 오른쪽 특이 벡터인 LxL의 직교 행렬이다.
이 때, 분해해야 할 더 큰 행렬 B*가
B*=(B|E)
로 주어 졌다.
E는 P개의 추가 (이미지)데이터가 있는 MxP 크기의 행렬이다. B*의 특이값 분해를 수행하자.
(Step 1) 행렬 (U|E)에 직교화(orthonormalization) 과정(주3)을 수행한다. 직교 행렬 Ua = (U|Et)를 만든다. E에서 Et를 구하기 위해 QR분해[1]를 이용한다. (참고: 표기에서 Et(E tilt)는 E'(transpose)가 아님. 직교 행렬임)
(Step 2) (L+P)x(L+P)크기의 Va행렬을 만든다.
Ip는 PxP의 단위 행렬이다. Va는 직교 행렬이다.
(Step 3) B*=Ua.Da.Va' 이므로, Da = Ua'.B*.Va 이다. 대입하면
이다.
Et'.B.V = 0(주1)이고, U'.B.V=D이다. 따라서
이다. Da는 D 행렬에 오른쪽 P개의 열이 추가된 새로운 행렬이다.
Da 행렬의 각 요소에 대해 살펴 보자.
.D는 기저 U로 표현 했을 때의 B 행렬에 대한 좌표 값이다.
.U'E는 새 데이터 E가 U기저에 투영된 좌표이다. 즉, U기저에서 E의 표현이다.
.Et'E는 새 데이터가 Et 기저에 투영된 좌표이다. 즉, Et기저에서 E의 표현이다.
즉, 새로 추가된 데이터를 새로 확장된 기저 (U|Et)에 대해 표현한 좌표 값이 (U'E; Et'E)이다. 이를 다시 쓰면 (U'; Et').E이다.
대각요소가 아닌 요소들이 있으므로 Da를 SVD하자.
Da = Ub.Db.Vb'
이고, B*의 SVD는
B* = Ua.Da.Va'= Ua.(Ub.Db.Vb').Va' = (Ua.Ub).Db.(Va.Vb)'
이다.
여기서, Db는 (L+P)x(L+P)의 대각 형렬이다. (Ua.Ub)는 Mx(L+P)의 column-orthonormal 행렬이고, (Va.Vb)는 (L+P)x(L+P)크기의 orthonormal 행렬이다.
상기한 바와 같이 B*의 SVD는 Ua, Da, Va의 계산과 작은 크기인 Da행렬의 SVD만을 요구한다. Ua, Da, Va의 계산은 B행렬의 SVD에서 얻어진 데이터를 사용한다.
(예1) 간단한 행렬에서 테스트
>> clear; clc;
>> A = rand(10,5); % 행렬 A를 10x5크기로 생성(old 데이터로 가정)
>> B = rand(10,4); % 행렬 B를 생성 (new 데이터로 가정)
>> [u0,d0,v0] = svd([A B]); % 두 행렬의 결합 확장 행렬의 특이값 분해
% 기존 데이터와 새 데이터의 모든 데이터가 결합된 행렬의 기저벡터를 얻을 수 있음
% A, B가 합해진 전체 데이터이고, 이에 대한 u0, d0가 필요.
% 그런데, 비디오처리를 할 경우, 계속 새 영상
% 데이터가 입력된다고 보면, [A B C ...]로 행렬 계속 커짐
% 큰 행렬의 svd는 계산량, 메모리 문제가 있음
%*** (1) skl 변환 [2]
% 원래의 고전적 방법
>> [ua,da,va]=svd(A); % A를 특이값 분해 (기존 데이터)
>> [q,r]=qr([ua*da B]); % new 데이터 B의 직교벡터 추출
% q행렬 내에는 직교벡터 ua 5개와 B의 직교벡터 4개가 얻어짐
% A대신 직교행렬과 특이행렬의 곱 ua*da 넣음
>> [u1,d1,v1]=svd(r); % q*r = q*u1*d1*v1'이므로
% *** u0 == q*u1, d0 == d1 ***
% A에서 ua, da만 사용했으므로 q*r는 [A B]가 아님
>> ua=q*u1; % 다음 단계를 위해 ua 갱신
>> da=d1; % 다음 단계를 위해 da 갱신
% 즉 기존 데이터의 u,d만 유지하면 새 데이터가 입력될 때 계산량, 메모리 문제 해결
%*** (2) 보다 개선된 방법
% [ua*da B]를 qr분해하는 것은 비효율적
% (이미 ua는 직교행렬로 분해되어 있기 때문)
% 따라서 추가되는 B만 qr분해 하면 됨
>> o1=B-ua*ua'*B; % B를 ua기저에 사영하고, 기저 ua로 표현되지 않는 B의 요소 추출(주4)
>> [q1,r1]=qr(o1); % 잔차를 가지고 새로운 직교벡터 q1을 얻음
>> rr=[da ua'*B; zeros(10,5) q1*o1];
% ua'*B: B를 ua기저에서 표현, q1*o1: 잔차의 q1기저 표현
>> [u2,d2,v2]=svd(rr); % rr(20x20)
>> uu = [ua q1]*u2; % uu=10x20. 결합 데이터의 기저
>> dd = d2; % u2=20x20. 결합 데이터의 특이값
>> ua = uu(:,1:10); % u0와 같음
>> da = dd(1:10,1:9); % d0와 같음
위 알고리즘은 새로운 데이터가 도착해서 이용가능한 데이터가 늘어났을 때, 기존 행렬과 새로 추가되는 행렬 모두에 대해 zero mean을 가정한다. 즉, 평균 값의 변화를 고려하지 않았다. 따라서 sample mean의 변화를 고려하면서 주축을 갱신하는 것이 필요하다.
다시 말하면, C=[A B]행렬의 특이 값 분해를 할 때, C행렬에 대해 C의 평균을 제거한 후, 특이값 분해를 수행해야 한다. A 및 B행렬의 특이 값 분해도 마찬가지이다. 그런데, A행렬, B행렬의 평균은 모두 C행렬의 평균과 다르며 제 각각이다. 즉, 평균은 항상 달라진다. 이를 고려하여 특이값 분해를 할 필요가 있다.
%*** (3) 평균 변화를 고려한 방법 [3]
% 평균 변화를 고려한 방법. 새로운 데이터가 추가되면 (즉, 기존 A
% 행렬에 B행렬이 추가되어 C행렬이 되었다면
% A평균, B평균, C평균 값이 다 다르다) 평균도 달라지므로 이를
% 고려해야 한다 [3].
>> C = [A B];
>> n = size(A,2); m = size(B,2);
>> ma = mean(A,2); mb=mean(B,2);
>> mc = n*ma/(n+m) + m*mb/(n+m);
>> Ma = repmat(ma,[1,5]);
>> Mb = repmat(mb,[1,4]);
>> Mc = repmat(mc,[1,9]);
>> Am = A-Ma; % A 행렬을 zero-mean으로
>> Bm = B-Mb; % B "
>> Cm = C-Mc; % C "
>>
>> [u0,d0,v0]=svd(Cm); % 목표값(true).
% 즉, Am, Bm에 skl 적용으로 u0, d0를 구해야 함.
참고문헌 [3]은 평균변화를 고려한 Sequential KL변환에 대해 잘 설명하고 있다. 참고문헌에서 scatter matrix Sc가 나오는데 covariance matrix와 다음의 관계가 있다:
Sc = (n+m).cov(C);
cov(C) = Sc/(n+m);
cov(C) = CC'=(UDV')(VDU')=U(D^2)U';
따라서 C행렬의 scatter 행렬 Sc의 svd를 구하는 것은 어떤 행렬 C의 U와 D^2을 구하는 것과 같다.
그런데, 참고문헌 [3]의 appendix에 따르면 평균이 제거된 C행렬의 scatter행렬은
Sc = Sa+Sb+ sqrt(.)(Ib-Ia)
= [(A-Ia) (B-Ib) sqrt(.)(Ib-Ia)] [(A-Ia) (B-Ib) sqrt(.)(Ib-Ia)]'
= Bt.Bt';
임이 증명되어 있다.
즉, C의 제곱(C.C')이 행렬 Bt의 제곱으로 이루어져 있다. 따라서 C의 svd는 Bt행렬의 svd로 구하면 된다. Bt 행렬은 (A-Ia)부분과 [(B-Ib) sqrt(.)(Ib-Ia)]의 두 부분의 결합으로 이루어져 있으며, skl을 적용할 수 있다. 즉, (A-Ia)의 svd 결과와 추가되는 [B-Ib sqrt(.)(.)]행렬의 svd를 결합한다.
>> [ua,da,va]=svd(Am);
>> v1 = sqrt((n*m)/(n+m))*(mb-ma);
>> Bh = [Bm v1];
>> o1 = Bh-ua*ua'*Bh; % Bh를 기저 ua로 표현 후의 잔차
>> [q1,r1]=qr(o1); % 잔차에서 새로운 기저 추출
>> rr = [da ua'*Bh; zeros(10,5) q1*o1];
>> [u2,d2,v2]=svd(rr);
>> dd = d2;
>> uu = [ua q1]*u2;
% C는 열이 9개이나 d2의 의미 있는(0이 아닌) 특이값이 9개가 안되는 경우
% 이 경우는 8개만 나옴
% uu의 1~8열까지만 유의미 함
>> dd = diag(dd);
>> cutoff = sum(dd.^2)*1e-6; % 임계치 정의
>> keep = find(dd.^2>=cutoff);
>> da = dd(:,keep);
>> ua = uu(:,keep);
(주1) Et'BV = Et'(UDV')V = (Et'U)D(V'V) = (0)D=0
(주2) SVD(X,0) 명령은 SVD(X)보다 더 작은 크기의 U, D, V행렬을 만든다. 만일 크기 mxn의 X 행렬에서 m>n이면 U의 첫번째 n열만 계산한다. S의 크기는 nxn이다. 만일 m<=n이면 SVD(X,0)는 SVD(X)와 같다.
(주3) [U|E]행렬을 직교화해야 한다. U는 이미 직교행렬이다. E는 새로 추가된 행렬이고 직교행렬이 아니다. 직교화를 위해 qr분해를 이용하는데, qr분해는 기존 벡터들에 수직한 새로운 벡터를 추가해 나가는 방식으로 직교 벡터를 만들어 간다. svd와 qr분해에서 만들어내는 직교 벡터들 사이의 관계를 통해 이해한다 (보충 설명 필요).
(주4) o1=B-ua*ua'*B: 기저 ua로 표현되지 않는 B의 요소 추출한다는 의미를 살펴보자. svd를 통한 행렬의 주축은 데이터의 산포(scatter)가 큰 몇 개의 축을 선택한다는 의미이다. 이 축들의 아닌 방향으로는 데이터의 산포가 없다는 의미이다. 그런데, 새로운 데이터 B는 기존 주축으로 표현이 되는 데이터도 있고 그러지 않은 데이터도 있을 것이다. 따라서, 기저 ua를 따르지 않는 데이터는 이를 표현하기 위한 새로운 주축이 필요하다. 즉, 잔차 o1에 대한 주축을 새로 선정하면 된다.
[1] QR분해: A행렬을 QR 분해하면 각 열이 서로 수직이고 열벡터의 크기가 1인 직교(orthonormal) Q행렬과 상삼각행렬(upper triangular matrix) R의 곱으로 분해 할 수 있다. 즉, A=Q.R을 얻는다.
[2] A. Levy and M. Lindenbaum, Sequential Karhunen-Loeve Basis Extraction and its application to Images, IEEE Trans. Image Processing 9(8), 2000.
[3] D.A. Ross, J. Lim, et al., Incremental Learning for Robust Visual Tracking
물체의 외관은 추적 동안 여러가지 원인으로 크게 변한다. PCA 기법을 통해 물체를 추적할 때, 물체를 표현하는 주축(principal axis)을 효과적으로 갱신하는 것이 필요하다.
연속 영상에서 물체 추적동안 새로 얻어지는 프레임들을 고려하여 주축을 갱신할 필요가 있을 때 sequential KL 알고리즘을 사용할 수 있다. (Sequential PCA 또는 Sequential visual tracking으로 부를 수도 있다.)
skl 변환은 sequential kl 변환이다. 연속 신호에서 데이터가 차례로 도착한 것을 처리하거나, 큰 행렬의 SVD를 계산해야 할 때, 큰 행렬을 작은 행렬로 분할하여 여러번 반복 처리할 때 사용한다.
예를 들면, A'=[A E] 행렬을 생각해 보자. A행렬은 10장의 영상 템플릿으로 이루어진 행렬이다. E행렬은 새로 매칭된 5장의 프레임에서 추출된 5장의 영상 템플릿으로 이루어진 행렬이다.
이 때, 이용할 수 있는 템플릿이 모두 15개가 되었으므로, 이것 모두가 포함된 A' 행렬에 대해 새로 주축을 구하는 것이 필요하다. 그런데, A행렬에 대해서는 특이값 분해를 통해 이미 구해 놓은 주축이 있다고 가정하자.
효율적 계산을 위해서는 A'행렬을 다시 특이값 분해를 하는 일 없이, 새롭게 추가된 E행렬에 대해서만 특이값 분해를 하고, 그 결과를 A행렬에 대한 특이값 분해의 결과와 결합한다. 이렇게 결합된 결과가 A'행렬에 대한 특이값 분해의 값과 동일하도록 만드는 일이 SKL이다.
R-SVD[2]는 A의 분해에 기반하여 효과적으로 행렬 A'의 분해 A' = [A|E] = U'D'V''를 계산할 수 있다. [U,D,V]=svd(A)이라고 가정하자. 이 때, A'의 특이값 분해는
 |
MxL크기의 행렬 B가 있다. 특이값 분해하면
B = UDV'
이다. U는 왼쪽 특이 벡터들로 구성된 MxL크기의 직교(orthonormal) 행렬이다. D는 LxL 크기의 행렬이다. 대각요소가 B행렬의 특이 값이다.
V는 각 열이 오른쪽 특이 벡터인 LxL의 직교 행렬이다.
이 때, 분해해야 할 더 큰 행렬 B*가
B*=(B|E)
로 주어 졌다.
E는 P개의 추가 (이미지)데이터가 있는 MxP 크기의 행렬이다. B*의 특이값 분해를 수행하자.
(Step 1) 행렬 (U|E)에 직교화(orthonormalization) 과정(주3)을 수행한다. 직교 행렬 Ua = (U|Et)를 만든다. E에서 Et를 구하기 위해 QR분해[1]를 이용한다. (참고: 표기에서 Et(E tilt)는 E'(transpose)가 아님. 직교 행렬임)
(Step 2) (L+P)x(L+P)크기의 Va행렬을 만든다.
Ip는 PxP의 단위 행렬이다. Va는 직교 행렬이다.
(Step 3) B*=Ua.Da.Va' 이므로, Da = Ua'.B*.Va 이다. 대입하면
이다.
Et'.B.V = 0(주1)이고, U'.B.V=D이다. 따라서
이다. Da는 D 행렬에 오른쪽 P개의 열이 추가된 새로운 행렬이다.
Da 행렬의 각 요소에 대해 살펴 보자.
.D는 기저 U로 표현 했을 때의 B 행렬에 대한 좌표 값이다.
.U'E는 새 데이터 E가 U기저에 투영된 좌표이다. 즉, U기저에서 E의 표현이다.
.Et'E는 새 데이터가 Et 기저에 투영된 좌표이다. 즉, Et기저에서 E의 표현이다.
즉, 새로 추가된 데이터를 새로 확장된 기저 (U|Et)에 대해 표현한 좌표 값이 (U'E; Et'E)이다. 이를 다시 쓰면 (U'; Et').E이다.
대각요소가 아닌 요소들이 있으므로 Da를 SVD하자.
Da = Ub.Db.Vb'
이고, B*의 SVD는
B* = Ua.Da.Va'= Ua.(Ub.Db.Vb').Va' = (Ua.Ub).Db.(Va.Vb)'
이다.
여기서, Db는 (L+P)x(L+P)의 대각 형렬이다. (Ua.Ub)는 Mx(L+P)의 column-orthonormal 행렬이고, (Va.Vb)는 (L+P)x(L+P)크기의 orthonormal 행렬이다.
상기한 바와 같이 B*의 SVD는 Ua, Da, Va의 계산과 작은 크기인 Da행렬의 SVD만을 요구한다. Ua, Da, Va의 계산은 B행렬의 SVD에서 얻어진 데이터를 사용한다.
(예1) 간단한 행렬에서 테스트
>> clear; clc;
>> A = rand(10,5); % 행렬 A를 10x5크기로 생성(old 데이터로 가정)
>> B = rand(10,4); % 행렬 B를 생성 (new 데이터로 가정)
>> [u0,d0,v0] = svd([A B]); % 두 행렬의 결합 확장 행렬의 특이값 분해
% 기존 데이터와 새 데이터의 모든 데이터가 결합된 행렬의 기저벡터를 얻을 수 있음
% A, B가 합해진 전체 데이터이고, 이에 대한 u0, d0가 필요.
% 그런데, 비디오처리를 할 경우, 계속 새 영상
% 데이터가 입력된다고 보면, [A B C ...]로 행렬 계속 커짐
% 큰 행렬의 svd는 계산량, 메모리 문제가 있음
%*** (1) skl 변환 [2]
% 원래의 고전적 방법
>> [ua,da,va]=svd(A); % A를 특이값 분해 (기존 데이터)
>> [q,r]=qr([ua*da B]); % new 데이터 B의 직교벡터 추출
% q행렬 내에는 직교벡터 ua 5개와 B의 직교벡터 4개가 얻어짐
% A대신 직교행렬과 특이행렬의 곱 ua*da 넣음
>> [u1,d1,v1]=svd(r); % q*r = q*u1*d1*v1'이므로
% *** u0 == q*u1, d0 == d1 ***
% A에서 ua, da만 사용했으므로 q*r는 [A B]가 아님
>> ua=q*u1; % 다음 단계를 위해 ua 갱신
>> da=d1; % 다음 단계를 위해 da 갱신
% 즉 기존 데이터의 u,d만 유지하면 새 데이터가 입력될 때 계산량, 메모리 문제 해결
%*** (2) 보다 개선된 방법
% [ua*da B]를 qr분해하는 것은 비효율적
% (이미 ua는 직교행렬로 분해되어 있기 때문)
% 따라서 추가되는 B만 qr분해 하면 됨
>> o1=B-ua*ua'*B; % B를 ua기저에 사영하고, 기저 ua로 표현되지 않는 B의 요소 추출(주4)
>> [q1,r1]=qr(o1); % 잔차를 가지고 새로운 직교벡터 q1을 얻음
>> rr=[da ua'*B; zeros(10,5) q1*o1];
% ua'*B: B를 ua기저에서 표현, q1*o1: 잔차의 q1기저 표현
>> [u2,d2,v2]=svd(rr); % rr(20x20)
>> uu = [ua q1]*u2; % uu=10x20. 결합 데이터의 기저
>> dd = d2; % u2=20x20. 결합 데이터의 특이값
>> ua = uu(:,1:10); % u0와 같음
>> da = dd(1:10,1:9); % d0와 같음
위 알고리즘은 새로운 데이터가 도착해서 이용가능한 데이터가 늘어났을 때, 기존 행렬과 새로 추가되는 행렬 모두에 대해 zero mean을 가정한다. 즉, 평균 값의 변화를 고려하지 않았다. 따라서 sample mean의 변화를 고려하면서 주축을 갱신하는 것이 필요하다.
다시 말하면, C=[A B]행렬의 특이 값 분해를 할 때, C행렬에 대해 C의 평균을 제거한 후, 특이값 분해를 수행해야 한다. A 및 B행렬의 특이 값 분해도 마찬가지이다. 그런데, A행렬, B행렬의 평균은 모두 C행렬의 평균과 다르며 제 각각이다. 즉, 평균은 항상 달라진다. 이를 고려하여 특이값 분해를 할 필요가 있다.
%*** (3) 평균 변화를 고려한 방법 [3]
% 평균 변화를 고려한 방법. 새로운 데이터가 추가되면 (즉, 기존 A
% 행렬에 B행렬이 추가되어 C행렬이 되었다면
% A평균, B평균, C평균 값이 다 다르다) 평균도 달라지므로 이를
% 고려해야 한다 [3].
>> C = [A B];
>> n = size(A,2); m = size(B,2);
>> ma = mean(A,2); mb=mean(B,2);
>> mc = n*ma/(n+m) + m*mb/(n+m);
>> Ma = repmat(ma,[1,5]);
>> Mb = repmat(mb,[1,4]);
>> Mc = repmat(mc,[1,9]);
>> Am = A-Ma; % A 행렬을 zero-mean으로
>> Bm = B-Mb; % B "
>> Cm = C-Mc; % C "
>>
>> [u0,d0,v0]=svd(Cm); % 목표값(true).
% 즉, Am, Bm에 skl 적용으로 u0, d0를 구해야 함.
참고문헌 [3]은 평균변화를 고려한 Sequential KL변환에 대해 잘 설명하고 있다. 참고문헌에서 scatter matrix Sc가 나오는데 covariance matrix와 다음의 관계가 있다:
Sc = (n+m).cov(C);
cov(C) = Sc/(n+m);
cov(C) = CC'=(UDV')(VDU')=U(D^2)U';
따라서 C행렬의 scatter 행렬 Sc의 svd를 구하는 것은 어떤 행렬 C의 U와 D^2을 구하는 것과 같다.
그런데, 참고문헌 [3]의 appendix에 따르면 평균이 제거된 C행렬의 scatter행렬은
Sc = Sa+Sb+ sqrt(.)(Ib-Ia)
= [(A-Ia) (B-Ib) sqrt(.)(Ib-Ia)] [(A-Ia) (B-Ib) sqrt(.)(Ib-Ia)]'
= Bt.Bt';
임이 증명되어 있다.
즉, C의 제곱(C.C')이 행렬 Bt의 제곱으로 이루어져 있다. 따라서 C의 svd는 Bt행렬의 svd로 구하면 된다. Bt 행렬은 (A-Ia)부분과 [(B-Ib) sqrt(.)(Ib-Ia)]의 두 부분의 결합으로 이루어져 있으며, skl을 적용할 수 있다. 즉, (A-Ia)의 svd 결과와 추가되는 [B-Ib sqrt(.)(.)]행렬의 svd를 결합한다.
>> [ua,da,va]=svd(Am);
>> v1 = sqrt((n*m)/(n+m))*(mb-ma);
>> Bh = [Bm v1];
>> o1 = Bh-ua*ua'*Bh; % Bh를 기저 ua로 표현 후의 잔차
>> [q1,r1]=qr(o1); % 잔차에서 새로운 기저 추출
>> rr = [da ua'*Bh; zeros(10,5) q1*o1];
>> [u2,d2,v2]=svd(rr);
>> dd = d2;
>> uu = [ua q1]*u2;
% C는 열이 9개이나 d2의 의미 있는(0이 아닌) 특이값이 9개가 안되는 경우
% 이 경우는 8개만 나옴
% uu의 1~8열까지만 유의미 함
>> dd = diag(dd);
>> cutoff = sum(dd.^2)*1e-6; % 임계치 정의
>> keep = find(dd.^2>=cutoff);
>> da = dd(:,keep);
>> ua = uu(:,keep);
(주1) Et'BV = Et'(UDV')V = (Et'U)D(V'V) = (0)D=0
(주2) SVD(X,0) 명령은 SVD(X)보다 더 작은 크기의 U, D, V행렬을 만든다. 만일 크기 mxn의 X 행렬에서 m>n이면 U의 첫번째 n열만 계산한다. S의 크기는 nxn이다. 만일 m<=n이면 SVD(X,0)는 SVD(X)와 같다.
(주3) [U|E]행렬을 직교화해야 한다. U는 이미 직교행렬이다. E는 새로 추가된 행렬이고 직교행렬이 아니다. 직교화를 위해 qr분해를 이용하는데, qr분해는 기존 벡터들에 수직한 새로운 벡터를 추가해 나가는 방식으로 직교 벡터를 만들어 간다. svd와 qr분해에서 만들어내는 직교 벡터들 사이의 관계를 통해 이해한다 (보충 설명 필요).
(주4) o1=B-ua*ua'*B: 기저 ua로 표현되지 않는 B의 요소 추출한다는 의미를 살펴보자. svd를 통한 행렬의 주축은 데이터의 산포(scatter)가 큰 몇 개의 축을 선택한다는 의미이다. 이 축들의 아닌 방향으로는 데이터의 산포가 없다는 의미이다. 그런데, 새로운 데이터 B는 기존 주축으로 표현이 되는 데이터도 있고 그러지 않은 데이터도 있을 것이다. 따라서, 기저 ua를 따르지 않는 데이터는 이를 표현하기 위한 새로운 주축이 필요하다. 즉, 잔차 o1에 대한 주축을 새로 선정하면 된다.
[1] QR분해: A행렬을 QR 분해하면 각 열이 서로 수직이고 열벡터의 크기가 1인 직교(orthonormal) Q행렬과 상삼각행렬(upper triangular matrix) R의 곱으로 분해 할 수 있다. 즉, A=Q.R을 얻는다.
[2] A. Levy and M. Lindenbaum, Sequential Karhunen-Loeve Basis Extraction and its application to Images, IEEE Trans. Image Processing 9(8), 2000.
[3] D.A. Ross, J. Lim, et al., Incremental Learning for Robust Visual Tracking
KL 변환
Karhunen-Loeve (KL) transform으로 불린다. 저 차원의 부분 공간으로 아주 큰 영상이나 벡터의 집합을 근사하는 방법이다.
mxn 크기의 2차원 배열로 주어지는 영상은 M(=mxn) 크기의 벡터 하나로 바꿀 수 있다. 만일 영상이 N개가 있다면, 전체 데이터의 크기는 MxN이다.
MxN 데이터에 대해 기저를 찾아보자.
일반적으로 MxN 벡터 집합에 대한 KL 기저(basis)를 계산하는 것은 O(MN^2)의 연산과 O(MN)의 메모리를 요구한다. 이러한 계산량은 데이터의 실시간 처리를 위해서는 상당히 크다.
KL 변환: N개의 (images) 벡터 (Ai),i=1...N가 주어질 때, K (K<N)개의 직교(orthonomal) 기저벡터의 집합을 찾는 것이다. 단, 기저 Uj에 의해 확장(span)된 공간에 사영(projection)된 원 데이터의 좌표로 재생한 데이터는 원 데이터를 잘 근사해야 한다.
행렬 A는 크기가 MxN (M>>N)이다.
A의 특이값 분해(SVD: Singular Value Decomposition)를 이용하여 기저벡터를 얻는다. 이때, KL 기저는 (특이값이 0이 되지 않는) K개의 왼쪽 특이 벡터(U의 첫번째 K개의 열)이다.
mxn 크기의 2차원 배열로 주어지는 영상은 M(=mxn) 크기의 벡터 하나로 바꿀 수 있다. 만일 영상이 N개가 있다면, 전체 데이터의 크기는 MxN이다.
MxN 데이터에 대해 기저를 찾아보자.
일반적으로 MxN 벡터 집합에 대한 KL 기저(basis)를 계산하는 것은 O(MN^2)의 연산과 O(MN)의 메모리를 요구한다. 이러한 계산량은 데이터의 실시간 처리를 위해서는 상당히 크다.
KL 변환: N개의 (images) 벡터 (Ai),i=1...N가 주어질 때, K (K<N)개의 직교(orthonomal) 기저벡터의 집합을 찾는 것이다. 단, 기저 Uj에 의해 확장(span)된 공간에 사영(projection)된 원 데이터의 좌표로 재생한 데이터는 원 데이터를 잘 근사해야 한다.
행렬 A는 크기가 MxN (M>>N)이다.
A의 특이값 분해(SVD: Singular Value Decomposition)를 이용하여 기저벡터를 얻는다. 이때, KL 기저는 (특이값이 0이 되지 않는) K개의 왼쪽 특이 벡터(U의 첫번째 K개의 열)이다.
R-SVM[1]은 SVD를 수정한 것으로 M>>N일 때 효과적이다:
(1) A = QR: qr분해를 통해 (열)직교(column-orthonormal) Q(MxN) 행렬, 상 삼각 행렬 R(NxN)을 얻는다. 계산량은 4MN^2이다.
(2) R = Ua.Da.Va': 표준 SVD알고리즘으로 R을 특이값 분해 한다. 계산량은 O(N^3)이다.
(3) U = Q.Ua, 계산량은 2MN^2이다.
정리해서 다시 쓰면,
A = QR = Q(Ua.Da.Va') = (Q.Ua).Da.Va' = UDV'
이다. 따라서, U = Q.Ua이다.
계산량을 살펴보면 (1)과 (3)이 6MN^2이고, (2)가 N^3이다. 즉 (1),(3)번이 좌우하고 (2)는 무시할 만하다.
(예1) 간단하게 테스트해 보자.
>> clear; clc;
>> A = rand(10,9); % 행렬 A를 10x9크기로 생성
>> [u0,d0,v0] = svd(A); % 행렬의 특이값 분해
>>
>> [q,r] = qr(A); % A 행렬의 QR 분해
>> [u1,d1,v1] = svd(r); % r의 특이값 분해
>> q*u1 % u0와 같음, 즉, u0.d0.v0' = (q*u1).d1.v1'이다.
>> d1 % d0와 같음
(1) A = QR: qr분해를 통해 (열)직교(column-orthonormal) Q(MxN) 행렬, 상 삼각 행렬 R(NxN)을 얻는다. 계산량은 4MN^2이다.
(2) R = Ua.Da.Va': 표준 SVD알고리즘으로 R을 특이값 분해 한다. 계산량은 O(N^3)이다.
(3) U = Q.Ua, 계산량은 2MN^2이다.
정리해서 다시 쓰면,
A = QR = Q(Ua.Da.Va') = (Q.Ua).Da.Va' = UDV'
이다. 따라서, U = Q.Ua이다.
계산량을 살펴보면 (1)과 (3)이 6MN^2이고, (2)가 N^3이다. 즉 (1),(3)번이 좌우하고 (2)는 무시할 만하다.
(예1) 간단하게 테스트해 보자.
>> clear; clc;
>> A = rand(10,9); % 행렬 A를 10x9크기로 생성
>> [u0,d0,v0] = svd(A); % 행렬의 특이값 분해
>>
>> [q,r] = qr(A); % A 행렬의 QR 분해
>> [u1,d1,v1] = svd(r); % r의 특이값 분해
>> q*u1 % u0와 같음, 즉, u0.d0.v0' = (q*u1).d1.v1'이다.
>> d1 % d0와 같음
[1] G. Golub and C. Van Loan, Matrix Computation. Johns Hopkins Univ. Press, 1992.
2013년 11월 7일 목요일
Camera calibration의 이해
camera parameter들을 더 잘 이해 하기 위해 simple calibration을 수행해 보기로 한다.
먼저 가로와 세로의 길이가 알려져 있는 평평한 물체(예를 들면 책)를 카메라 앞에 두고 사진을 찍는다.
사진을 찍을 때, 책의 모서리가 영상에 수평과 수직으로 나타나게 하고, 기울어 지지 않도록 한다. 책의 중심이 카메라의 광 축을 지나가도록 한다.
또한 카메라와 책 사이의 거리를 미리 측정해 둔다.
(images from hhkim)
예를 들면 아래 이미지에서 책의 가로 길이는 100mm, 높이는 130mm, 카메라까지의 거리는 500mm이다.
dX=130, dY=100, dZ=500
(images from hhkim)
이미지에서 책의 크기를 픽셀 단위로 측정하면 dx=494 pixel, dy=363 pixel이다.
3차원 공간에 있는 한 점이 이미지 면에 맺히는 사영(projection)을 나타내는 식은
x=fx(X/Z),
y=fy(Y/Z)
이다. 여기서 (X, Y, Z)'는 공간 점의 좌표이다. (x, y)'는 이미지 좌표이다. 이 식을 다시 쓰면,
fx = (dx/dX).dZ,
fy = (dy/dY).dZ
이다. 위 이미지의 측정 값들을 대입하면,
fx = (494/130).500 = 1900,
fy = (363/100).500 = 1819
이다.
그런데, 이미지 좌표와 카메라 좌표사이의 관계인 x=fx(X/Z), y=fy(Y/Z)는 동차 좌표(homogeneous coordinate)로 (fx.X, fy.Y, Z)'이고 이 식은 행렬 식으로 다음과 같다.
(fx.X fy.Y Z)'=[fx 0 0 0; 0 fy 0 0; 0 0 1 0].(X Y Z 1)'
만일 이미지 좌표의 원점이 영상 중심이 아니라 (left,top)위치라면 중심까지의 offset 거리를 고려해 주어야 한다.
즉, 공간 좌표 (X Y Z)'는 이미지 좌표 (fX/Z+px, fY/Z+py)로 사영된다. (px, py)'는 이미지 중심 위치이다.
이미지 좌표를 행렬 식으로 나타내면
(fx.X+Z.px, fy.Y+Z.py, Z)'=[fx 0 px 0; 0 fy py 0; 0 0 1 0].(X Y Z 1)'
이다.
이 식을 다시 쓰면
x = K[I | 0]Xcam
이다. 여기서 x는 동차좌표로 표현된 이미지 좌표, K=[fx 0 py; 0 fy py; 0 0 1]인 calibration 행렬이고, I는 3x3의 단위 행렬, 0은 3x1의 0벡터, Xcam은 카메라 축 기준에서 공간 점의 좌표이다.
위의 테스트 이미지에서 해상도는 2592(가로)x1936(높이)이다. 따라서 이미지 센터 위치인 px, py를 이미지 중심이라 하면,
px = 1936/2 = 968,
py = 2592/2 = 1296
이다. 따라서,
K = [1900 0 968; 0 1819 1296; 0 0 1]
이다.
(예 1) 카메라 축의 원점을 기준으로 했을 때, 책 중심 위치의 좌표는 (0 0 500)'이다. 이 점의 이미지 위치를 구하시오.
[s.x, s.y, s]'= K.[I | 0].(0 0 500 1)'
이다.
값을 대입하고 계산하면
(s.x, s.y, s)' = [1900 0 968 0; 0 1819 1296 0; 0 0 1 0].(0 0 500 1)'
= (968x500, 1296x500, 500)'
이다. 동차 좌표이므로 마지막 항으로 나누어 주면
(x, y, 1)'=(968, 1296, 1)'
이다. 즉, 물체의 중심 점의 이미지 좌표는 px, py와 같아진다.
(예 2) 이미지 좌표 (0, 0)은 공간 위치가 어디인지 계산하시오.
(0, 0, 1)' = K. [I | 0].(X, Y, 500, 1)'
이다. K 행렬을 알고 있으므로 양변에 inv(K)를 곱해 정규화 좌표(inv(K).x)를 만든다.
inv(K).(0 0 1)' = [1 0 0 0; 0 1 0 0; 0 0 1 0].(X, Y, 500, 1)' = (X, Y, 500)'
이다.
inv(K).(0 0 1)' = (-0.5095, -0.7125, 1)' = (X/500, Y/500, 1)'
이므로, X=-254.75, Y=-356.25이다.
Z=500이므로,
(X, Y, Z)' = (-256.75, -356.25, 500)'
이다.
물론 이 위치는 카메라 축을 원점으로 했을 때의 좌표이다.
먼저 가로와 세로의 길이가 알려져 있는 평평한 물체(예를 들면 책)를 카메라 앞에 두고 사진을 찍는다.
사진을 찍을 때, 책의 모서리가 영상에 수평과 수직으로 나타나게 하고, 기울어 지지 않도록 한다. 책의 중심이 카메라의 광 축을 지나가도록 한다.
또한 카메라와 책 사이의 거리를 미리 측정해 둔다.
예를 들면 아래 이미지에서 책의 가로 길이는 100mm, 높이는 130mm, 카메라까지의 거리는 500mm이다.
dX=130, dY=100, dZ=500
이미지에서 책의 크기를 픽셀 단위로 측정하면 dx=494 pixel, dy=363 pixel이다.
3차원 공간에 있는 한 점이 이미지 면에 맺히는 사영(projection)을 나타내는 식은
x=fx(X/Z),
y=fy(Y/Z)
이다. 여기서 (X, Y, Z)'는 공간 점의 좌표이다. (x, y)'는 이미지 좌표이다. 이 식을 다시 쓰면,
fx = (dx/dX).dZ,
fy = (dy/dY).dZ
이다. 위 이미지의 측정 값들을 대입하면,
fx = (494/130).500 = 1900,
fy = (363/100).500 = 1819
이다.
그런데, 이미지 좌표와 카메라 좌표사이의 관계인 x=fx(X/Z), y=fy(Y/Z)는 동차 좌표(homogeneous coordinate)로 (fx.X, fy.Y, Z)'이고 이 식은 행렬 식으로 다음과 같다.
(fx.X fy.Y Z)'=[fx 0 0 0; 0 fy 0 0; 0 0 1 0].(X Y Z 1)'
만일 이미지 좌표의 원점이 영상 중심이 아니라 (left,top)위치라면 중심까지의 offset 거리를 고려해 주어야 한다.
즉, 공간 좌표 (X Y Z)'는 이미지 좌표 (fX/Z+px, fY/Z+py)로 사영된다. (px, py)'는 이미지 중심 위치이다.
이미지 좌표를 행렬 식으로 나타내면
(fx.X+Z.px, fy.Y+Z.py, Z)'=[fx 0 px 0; 0 fy py 0; 0 0 1 0].(X Y Z 1)'
이다.
이 식을 다시 쓰면
x = K[I | 0]Xcam
이다. 여기서 x는 동차좌표로 표현된 이미지 좌표, K=[fx 0 py; 0 fy py; 0 0 1]인 calibration 행렬이고, I는 3x3의 단위 행렬, 0은 3x1의 0벡터, Xcam은 카메라 축 기준에서 공간 점의 좌표이다.
위의 테스트 이미지에서 해상도는 2592(가로)x1936(높이)이다. 따라서 이미지 센터 위치인 px, py를 이미지 중심이라 하면,
px = 1936/2 = 968,
py = 2592/2 = 1296
이다. 따라서,
K = [1900 0 968; 0 1819 1296; 0 0 1]
이다.
(예 1) 카메라 축의 원점을 기준으로 했을 때, 책 중심 위치의 좌표는 (0 0 500)'이다. 이 점의 이미지 위치를 구하시오.
[s.x, s.y, s]'= K.[I | 0].(0 0 500 1)'
이다.
값을 대입하고 계산하면
(s.x, s.y, s)' = [1900 0 968 0; 0 1819 1296 0; 0 0 1 0].(0 0 500 1)'
= (968x500, 1296x500, 500)'
이다. 동차 좌표이므로 마지막 항으로 나누어 주면
(x, y, 1)'=(968, 1296, 1)'
이다. 즉, 물체의 중심 점의 이미지 좌표는 px, py와 같아진다.
(예 2) 이미지 좌표 (0, 0)은 공간 위치가 어디인지 계산하시오.
(0, 0, 1)' = K. [I | 0].(X, Y, 500, 1)'
이다. K 행렬을 알고 있으므로 양변에 inv(K)를 곱해 정규화 좌표(inv(K).x)를 만든다.
inv(K).(0 0 1)' = [1 0 0 0; 0 1 0 0; 0 0 1 0].(X, Y, 500, 1)' = (X, Y, 500)'
이다.
inv(K).(0 0 1)' = (-0.5095, -0.7125, 1)' = (X/500, Y/500, 1)'
이므로, X=-254.75, Y=-356.25이다.
Z=500이므로,
(X, Y, Z)' = (-256.75, -356.25, 500)'
이다.
물론 이 위치는 카메라 축을 원점으로 했을 때의 좌표이다.
2013년 11월 3일 일요일
연속 영상에서 분류의 결정
연속 영상에서 물체를 인식하거나 분류를 결정할 경우를 생각해 보자.
분류기가 어떤 결과를 줄 것이다. 결과는 1(물체) 또는 0(배경)이다.
영상은 항상 조명이 변하거나 노이즈가 영향을 끼쳐, 실제 물체가 존재한다고 해도 결과는 일관적(consistent)이지 않다.
예를 들면, 연속된 영상 내의 어떤 위치에 물체가 있다.
연속 10 프레임동안 감지를 수행한다고 가정하자.
결과가 1,1,0,1,1,1,0,0,1,1이다. 즉, 물체, 물체, 배경, 물체,... 이다.
첫 부분이 현재 프레임의 결과이고, 나머지는 이전 프레임의 결과들이다. 이 때, 현재 프레임의 결과 1은 과연 '물체'인지, '배경'인지를 판단해야 한다. AMC(additive Markov chain)은 분류 판단을 위해서 이전 프레임들을 참고하는 방법이다.
AMC는 현재 랜덤 변수의 상태가 연속된 이전 랜덤 변수의 상태에 영향을 받는 것을 표현한다 [1].
즉, 연속해서 여러 프레임 동안 '물체'로 감지되었고, 현재 프레임에서도 '물체'로 감지되었다면 실제 '물체'일 가능성이 높다.
AMC의 식은 다음과 같다:
여기서 z는
로 확률(pdf) 생성을 위한 정규화 값이다. xi는 0과 1의 값을 가진다.
f()는 메모리 함수로 이전 결과에 의해 정의된다.
이전의 5프레임(n=5)을 본다고 하면 f함수는 5개가 정의 된다:
f(x1,x0), f(x2,x1), f(x3, x2), f(x4, x3), f(x5, x4).
(xi,xi-1)에 대한 f함수는 다음과 같다:
두 항 중에서 앞의 항은 현 프레임의 결과를 고려하고, 뒤 항은 이전 프레임과의 상관성을 고려한다.
alpa값은 두 항 중 어느쪽에 가중치를 둘지를 결정한다. alpa의 값은 0~1사이 값이다. 또한 (1/2)^(n-i)는 각 프레임 시간에 대한 가중치인데 현재 프레임(n)에 가까울수록 큰 값으로 기여하게 된다.
(예제) 이전 5프레임을 고려한다. alpa=0.7이고 (x2,x1)=(1,1)일 때 f()값을 계산하시오.
x2 프레임에서 물체로 감지(x2=1)가 되고, 이전 프레임(x1)도 물체로 감지(x1=1)되었다. 가중치는 1/2^(n-2)=1/2^3=1/8이다.
각 항을 대입하면,
f(x2,x1)=0.7( (1/8).1 ) + 0.3( (1/8).(1&1) ) = 0.7*(1/8) + 0.3*(1/8)
이다.
예제를 보면, (xi,xi-1)의 값과 프레임 순서가 주어지면 미리 f()값을 계산할 수 있다. 계산된 값을 LUT(lookup table)에 저장해 놓고 필요 시에 사용할 수 있다.
LUT는 (xi,xi-1, tn)의 함수이다. tn은 프레임 순서를 나타낸다.
즉, (xi=1, xi-1=0)이고, n=5에서 첫번째로 가까운 프레임 이라면
f(1,0,1) = 0.7( (1/2)^1 x 1 ) + 0.3( (1/2)^1 x (1 & 0 )) = 0.7x(1/2)
이다. 이런 식으로 미리 (xi, xi-1, tn)의 모든 경우에 대해 계산해 놓으면 된다.
f(0,0,0) = 0
f(0,1,0) = 0
f(1,0,0) = 0.7
f(1,1,0) = 1
f(0,0,1) = 0
f(0,1,1) = 0.3( (1/2)x 0) = 0
f(1,0,1) = 0.7 x 0.5
f(1,1,1) = 0.7 (1/2) + 0.3(1/2) = 0.5
f(0,0,2) = 0
f(0,1,2) = 0
f(1,0,2) = 0.7 (1/4)
f(1,1,2) = 0.7 (1/4) + 0.3 ( (1/4) ) = 1/4
....
online 연산에서는 새로 f()를 계산하지 말고 LUT에 저장되어 있는 값을 호출해서 사용한다.
최종 분류의 결정은 다음과 같다:
분류기가 어떤 결과를 줄 것이다. 결과는 1(물체) 또는 0(배경)이다.
영상은 항상 조명이 변하거나 노이즈가 영향을 끼쳐, 실제 물체가 존재한다고 해도 결과는 일관적(consistent)이지 않다.
예를 들면, 연속된 영상 내의 어떤 위치에 물체가 있다.
연속 10 프레임동안 감지를 수행한다고 가정하자.
결과가 1,1,0,1,1,1,0,0,1,1이다. 즉, 물체, 물체, 배경, 물체,... 이다.
첫 부분이 현재 프레임의 결과이고, 나머지는 이전 프레임의 결과들이다. 이 때, 현재 프레임의 결과 1은 과연 '물체'인지, '배경'인지를 판단해야 한다. AMC(additive Markov chain)은 분류 판단을 위해서 이전 프레임들을 참고하는 방법이다.
AMC는 현재 랜덤 변수의 상태가 연속된 이전 랜덤 변수의 상태에 영향을 받는 것을 표현한다 [1].
즉, 연속해서 여러 프레임 동안 '물체'로 감지되었고, 현재 프레임에서도 '물체'로 감지되었다면 실제 '물체'일 가능성이 높다.
AMC의 식은 다음과 같다:
여기서 z는
로 확률(pdf) 생성을 위한 정규화 값이다. xi는 0과 1의 값을 가진다.
f()는 메모리 함수로 이전 결과에 의해 정의된다.
이전의 5프레임(n=5)을 본다고 하면 f함수는 5개가 정의 된다:
f(x1,x0), f(x2,x1), f(x3, x2), f(x4, x3), f(x5, x4).
(xi,xi-1)에 대한 f함수는 다음과 같다:
두 항 중에서 앞의 항은 현 프레임의 결과를 고려하고, 뒤 항은 이전 프레임과의 상관성을 고려한다.
alpa값은 두 항 중 어느쪽에 가중치를 둘지를 결정한다. alpa의 값은 0~1사이 값이다. 또한 (1/2)^(n-i)는 각 프레임 시간에 대한 가중치인데 현재 프레임(n)에 가까울수록 큰 값으로 기여하게 된다.
(예제) 이전 5프레임을 고려한다. alpa=0.7이고 (x2,x1)=(1,1)일 때 f()값을 계산하시오.
x2 프레임에서 물체로 감지(x2=1)가 되고, 이전 프레임(x1)도 물체로 감지(x1=1)되었다. 가중치는 1/2^(n-2)=1/2^3=1/8이다.
각 항을 대입하면,
f(x2,x1)=0.7( (1/8).1 ) + 0.3( (1/8).(1&1) ) = 0.7*(1/8) + 0.3*(1/8)
이다.
예제를 보면, (xi,xi-1)의 값과 프레임 순서가 주어지면 미리 f()값을 계산할 수 있다. 계산된 값을 LUT(lookup table)에 저장해 놓고 필요 시에 사용할 수 있다.
LUT는 (xi,xi-1, tn)의 함수이다. tn은 프레임 순서를 나타낸다.
즉, (xi=1, xi-1=0)이고, n=5에서 첫번째로 가까운 프레임 이라면
f(1,0,1) = 0.7( (1/2)^1 x 1 ) + 0.3( (1/2)^1 x (1 & 0 )) = 0.7x(1/2)
이다. 이런 식으로 미리 (xi, xi-1, tn)의 모든 경우에 대해 계산해 놓으면 된다.
f(0,0,0) = 0
f(0,1,0) = 0
f(1,0,0) = 0.7
f(1,1,0) = 1
f(0,0,1) = 0
f(0,1,1) = 0.3( (1/2)x 0) = 0
f(1,0,1) = 0.7 x 0.5
f(1,1,1) = 0.7 (1/2) + 0.3(1/2) = 0.5
f(0,0,2) = 0
f(0,1,2) = 0
f(1,0,2) = 0.7 (1/4)
f(1,1,2) = 0.7 (1/4) + 0.3 ( (1/4) ) = 1/4
....
online 연산에서는 새로 f()를 계산하지 말고 LUT에 저장되어 있는 값을 호출해서 사용한다.
최종 분류의 결정은 다음과 같다:
참고 문헌
[1] Fire detection system using random forest classification for image sequences, Optical Engineering, 2013.
2013년 11월 1일 금요일
두 축 사이의 변환 관계
공간에 두 개의 좌표 축이 있다:
1. 두 축에서 바라본 한 점 X의 좌표 사이의 관계는 다음과 같다.
이 식은 두 축 사이의 회전은 없다고 가정한 경우이다.
Oc에서 X로 가는 벡터는 Oc에서 Ow로 가서 다시 Ow에서 X로 가는 벡터의 합이다.
Ow 기준에서 보면 Oc의 좌표는 tw이다. 따라서 -tw 이동하면 두 축은 일치 된다.
(tc가 아니라 tw로 주어 졌다.)
좌표이동을 선형변환으로 표현하였다. tilt는 3-vector이다.
2. 회전을 표현하는 방법은 2가지 이다. 두 좌표의 자세가 다르다고 가정하자.
먼저 C축을 돌리는 경우이다. 즉, C축을 돌려서 W축과 일치시킨다. 여기서 Rc는 W축 기준으로 C 좌표의 축 방향이다.
다음, W축을 돌리는 경우이다. 즉, W축을 돌려서 C축과 일치시킨다. 여기서 Rw는 C축 기준으로 W좌표의 축 방향이다.
많은 문헌에서 이 표현을 주로 사용한다. 이동과 회전의 선형변환을 하나의 식으로 표현해 보자.
1. 두 축에서 바라본 한 점 X의 좌표 사이의 관계는 다음과 같다.
이 식은 두 축 사이의 회전은 없다고 가정한 경우이다.
Oc에서 X로 가는 벡터는 Oc에서 Ow로 가서 다시 Ow에서 X로 가는 벡터의 합이다.
Ow 기준에서 보면 Oc의 좌표는 tw이다. 따라서 -tw 이동하면 두 축은 일치 된다.
(tc가 아니라 tw로 주어 졌다.)
좌표이동을 선형변환으로 표현하였다. tilt는 3-vector이다.
2. 회전을 표현하는 방법은 2가지 이다. 두 좌표의 자세가 다르다고 가정하자.
먼저 C축을 돌리는 경우이다. 즉, C축을 돌려서 W축과 일치시킨다. 여기서 Rc는 W축 기준으로 C 좌표의 축 방향이다.
다음, W축을 돌리는 경우이다. 즉, W축을 돌려서 C축과 일치시킨다. 여기서 Rw는 C축 기준으로 W좌표의 축 방향이다.
많은 문헌에서 이 표현을 주로 사용한다. 이동과 회전의 선형변환을 하나의 식으로 표현해 보자.
% 두 좌표계 사이의 회전 행렬의 표현
% Coded by funMV. 12/2013
clear; clc;
% 두 좌표계 사이의 회전이 없을 경우
xw = [20,20,0]'; % coord of a point w.r.t world coord system
xt = [50 ,10, 100]'; % translational vector between world -> camera
xc1 = xw - xt;
% 만일 C축이 W축에 대해 R만큼 회전된 경우
a2p=pi/180;
Rc_w = rot(20*a2p, 10*a2p, -30*a2p); % C축 3개의 방향(w축에 상대적으로)
Rw_c = Rc_w'; % W축 3개의 방향(c축에 상대적으로)
xc2 = Rw_c*(xw-xt);
피드 구독하기:
글 (Atom)