SURF(Speed-Up Robust Features)

컴퓨터비전/영상처리 2015. 5. 19. 18:02

SURF는 2008년 Herbert Bay가 제안한 방법입니다.


[Bay 08] Bay, Herbert, et al. "Speeded-up robust features (SURF)." Computer vision and image understanding 110.3 (2008): 346-359.


2004년 SIFT의 등장 이후 물체의 descriptor에 대한 관심이 높아진 이후로 많은 연구가 있었습니다.

SIFT의 최대 단점인 연산의 복잡성을 해소 하기 위한 방법으로 한때 많은 관심(성능)을 보인 

SURF 기술자를 정리해 봅시다.


위 논문에서는 SIFT 외에도 다른 방법들과 비교를 했지만 SIFT 못지 않게 성능을 유지하며 계산 속도를

향상 시키고자 했습니다.


그래서 속도 향상을 위해 제안한 방법을 요약하면 크게 총 3가지로 정리할 수 있습니다.


1. Integral Image(적분 이미지) 이용

 SURF에서는 SIFT와는 다르게 적분 영상을 사용하여 관심점과 영역을 찾습니다. 관심점 계산을 위해 

고속 헤시안 검출(Fast Hessian Detector)를 사용하는데, 이때 적분 영상을 사용합니다.


적분 영상이란 말 그대로 영상을 적분하는것을 말합니다.

적분이란? 고등학교 학생 때의 기억을 더듬어 보면 어떤 영역의 넓이를 구했던 것을 생각해 볼 수 있습니다.

그렇다면 영상의 적분이란 영상이 갖는 영역의 넓이(밝기의 합)을 구하는 것을 추론할 수 있습니다.


영상의 밝기의 합은 그냥 선형으로 영역을 계산하면 되지 않을까 생각하기 쉽지만 많은 영역을

n^2으로 계산하면 연산속도가 매우 느립니다. 따라서 적분 영상을 처음에 하나 만들어 두고 영상의 부분합을

구하면 빠른 속도로 계산할 수 있습니다.


그렇다면 영상의 Integral Image는 어떻게 계산할까요?

아래 그림을 보면 한눈에 이해 할 수 있습니다.



위 그림의 왼쪽 이미지는 기존 이미지를 의미하고, 오른쪽은 적분 이미지를 의미합니다.

적분 이미지는 (0,0)부터 현재 좌표점까지의 밝기값의 누적으로 계산할 수 있습니다.

만약 (1,1)를 구할때는 (0,0) + (0,1) + (1,0) + (1,1)의 총합을 (1,1)에 저장하면 됩니다.


영상의 적분 이미지를 생성해 두고 영상의 부분 합을 계한하려고 하면 

D = D + A - B - C



3번의 더하기 빼기로 영역의 밝기 부분합을 구할 수 있습니다.



빨강 부분에 관한 합을 적분 부분합 식을 이용하면

 D = 21 + 4 - 11 - 11 = 3 

를 빠른 속도로 계산할 수 있습니다.


2. 축소한 detector와 descriptor를 사용하여 빠른 (차원수 축소)

 관심점을 계산하기 위해 SURF는 헤시안(Hessian) 행렬을 사용합니다.

논문에서는 정확성이 좋고 determinant가 최대값인 위치의 blob 구조를 검출하기 위해서라고 나와 있습니다.


[그림출처 : http://www.quora.com/What-are-the-benefits-of-calculating-the-eigenvector-of-a-Hessian-matrix-in-image-processing]


 만약 determinant(행렬식)가 음수이고 eigenvalue가 서로 다른 부호이면 해당 좌표는 관심점이 아니고

determinant가 양수이고, eigenvalue가 둘 다 음수 혹은 양수이면 해당 좌표는 관심점으로 판단합니다.

SURF의 경우 2차 미분을 Gaussian을 사용하는데, Scale-space를 분석하는데 가우시안이 가장 최적이라고

판단하여 사용함을 논문에 명시하였습니다. 따라서 x,y,xy에 대한 헤시안 행렬을 계산하면 관심점(Interest Point)

로 판단합니다.


또한 가우시안 분포 헤시안 행렬을 계산할때 계산을 좀 더 단순화 시키기 위해 근사화한 Dxx와 Dyy 박스필터를 

사용합니다.



 위 사진과 같이 근사화한 박스 필터를 사용하므로써 적분이미지 계산으로 영역의 넓이를 빠르게 구하고 

한 넓이로 헤시안 행렬을 계산해 관심점을 판단합니다.



논문에서 weight값은 0.9를 사용했습니다. 


 SURF도 스케일에 robust하기 위해 SIFT와 비슷하게 여러 스케일에서 관심점을 추출합니다. 

아래 그림을 보면 SIFT 필터 사이즈를 고정시키고, 이미지 사이즈를 줄여 관심점을 검출하는 반면, 

SURF는 SIFT와는 반대로 이미지 스케일을 고정시키고 필터 사이즈를 키워가며 관심점을 추출합니다.



SIFT와 다른 차이점을 둔 이유는 계산의 효율성 때문일 것 같습니다. 적분 이미지를 사용하기 때문에

위치를 변경하는것 만으로 필터 사이즈를 조절할 수 있기 때문입니다. 또한 down sampling하지 않기 

때문에 aliasing이 없는 장점이 있습니다.(scale invariance 가 제한적일 수 있음 - 단점)


아래 그림을 보면 3계의 옥타브가 있고 점점 스케일이 증가하는것을 확인할 수 있습니다. 

옥타브가 증가하면 필터 사이즈가 2배로 증가하는 것을 볼 수 있습니다.

(논문에서는 하나의 옥타브에 4번 scale up)


 

 선정한 관심점을 비 최대치 억제(non maximum suppression)을 통해 크기 공간에 대한 극대값(검출된 관심점)

을 선정합니다. 또한 이 극대값에 대해 보간법을 적용하여 보간된 관심점을 특징 벡터로 저장하여 사용합니다.

(Sub-Pixel로 보간)



 최종 보간된 관심점을 중심으로 6s(scale)의 원 안에 있는 이웃들에 관해 x, y 방향의 Haar wavelet response를 

계산합니다. 오른쪽 원 안의 파란 점이 response값의 분포입니다. 그리고 빨간 화살표는 부채꼴 모양 안의 영역에 있는

response의 합을 벡터로 표현한 것입니다. 60도 크기의 window를 슬라이딩 하면서 x,y축의 response의 합을 계산해서

벡터를 생성시키고 벡터 길이가 가장 긴 것이 orientation으로 할당됩니다.

이때 Haar wavelet filter를 이용하기 때문에 적분 이미지의 장점을 적용할 수 있습니다.(빠른 계산)



 위 그림을 보면 descriptor를 생성하기 위해 관심점 주위 20scale 크기의 window를 구성해 이미지를 회전에 

invariant 하기 위해 회전시킨 뒤 계산합니다. 이때 앞에서 구한 orientation를 활용하면 됩니다.

또한 descriptor window를 4×4로 나눈뒤, 나눠진 구역은 5×5 haar wavelet로 구성합니다.



 논문에서는 또한 contrast에 invariant하다는 것을 보여주기 위해 단순히 x, y 축에 대한 response값을 

합한 것이 아니라, 각 값에 대한 절대값을 합하여 contrast에 강하다는것을 증명하고 있습니다.


3. Contrast를 이용한 빠른 매칭

 hessian 행렬을 계산한 라플라시안 부호를 비교해 간단히 매칭을 시도합니다.(부호 비교)



ubuntu 14.04 LTS iTorch 설치

프로그래밍 2015. 4. 18. 03:56

iTorch 설치 가이드


다양한 포멧 출력 및 시각 효과(그래프) 등을 지원하는 iTorch 설치 방법을 알아보겠습니다.

iTorch 는 사진, 오디오, 비디오 등 다양한 형식 출력을 지원합니다.


사진


동영상


그래프


itorch를 설치하긴 위해선 연관 라이브러리를 설치해야 합니다.



  https://github.com/facebook/iTorch


에서 설치 방법을 보면 Torch-7를 설치해야 합니다.


Torch설치는 지난 포스팅을 참조하면 됩니다.


그리고 요구조건인 IPython을 설치해야 하는데 설치 방법은 다음과 같습니다.


1. IPython 설치


1.1 pip 명령어로 ipython을 설치

 $ pip install ipython


1.2 ipython-notebook을 설치

 $ pip install ipython-notebook


를 설치하면 설치가 완료됩니다.

 $ ipython --version


를 입력해서 만약 버전이 2.2 보다 낮으면 업데이트를 해야합니다.

업데이트를 하려면 다음과 같이 입력하면 됩니다.


 $ pip install --upgrade ipython[all]


업데이를 완료한 뒤 ipython --version를 입력하면 아래 사진과 같이 버전이 업데이트 된것을 볼 수 있습니다.



2. iTorch 설치


2.1 Dependencies 라이브러리를 다운받습니다.

 $ sudo apt-get install libzmq3-dev


2.2 git허브에서 소스를 다운받습니다.

 $ git clone https://github.com/facebook/iTorch.git


2.3 다운받은 폴더로 들어가 make합니다.

 $ cd iTorch


 $ luarocks make


만약 sudo 관리자로 설치했을 경우(혹은 global로 설치했을 경우)

다음과 같이 입력합니다.


 $ sudo env "PATH=$PATH" luarocks make

 $ sudo chown -R $USER $(dirname $(ipython locate profile))


이제 설치가 완료되었습니다.


2.4 실행

 $ itorch notebook # notebook mode

   or

 $ itorch             # consol mode
   or
 $ itorch qtconsole # Qt mode


2.5 예제

iTorch notebook를 이용해서 Demo를 테스트 해볼수 있습니다.(아래 주소 참조)

http://nbviewer.ipython.org/github/facebook/iTorch/blob/master/iTorch_Demo.ipynb


'프로그래밍' 카테고리의 다른 글

ubuntu 14.04 LTS Torch 설치  (0) 2015.04.17
Lua 프로그래밍  (0) 2015.04.16
안드로이드 가상 머신 지니모션(Genymotion) 설치  (0) 2015.04.09
멀티 스레드  (0) 2014.12.13
파이썬 argparse  (0) 2014.11.08

ubuntu 14.04 LTS Torch 설치

프로그래밍 2015. 4. 17. 15:45

1. 터미널 실행


2. 우분투 root 비밀번호 설정

 $ sudo passwd root


3. 관리자 로그인

 $ su

   암호 : ___


4. torch 및 관련 패키지 다운로드

 $ curl -s https://raw.githubusercontent.com/torch/ezinstall/master/install-all | bash


5. 실행

단순히 실행하려면

 $ luajit -ltorch

을 입력하면 아래와 같이 shell이 실행된다.



5.1 th interpreter 설정

 $ th -lparallel -loptim -lpl -limage

 $ th


'프로그래밍' 카테고리의 다른 글

ubuntu 14.04 LTS iTorch 설치  (2) 2015.04.18
Lua 프로그래밍  (0) 2015.04.16
안드로이드 가상 머신 지니모션(Genymotion) 설치  (0) 2015.04.09
멀티 스레드  (0) 2014.12.13
파이썬 argparse  (0) 2014.11.08

Lua 프로그래밍

프로그래밍 2015. 4. 16. 02:12


1. 변수 선언

1
2
3
4
5
6
7
num = 42 (숫자, 모든 숫자는 더블형으로 저장된다)
 
= 'hello world' (문자열)
 
= "double-quotes are also fine" (문자열)
 
t  =  nil (정의되지 않은 t 선언) 
cs


2. 반복문

1
2
3
4
5
while num < 50 do
 
    num = num + (+++=와 같음)
 
end
cs

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
//---------------------------------------------------------------
// 상승 반복문(1씩 증가)
 
upSum = 0
 
for i = 1100 do
   upSum = upSum + i
end
 
//---------------------------------------------------------------
// 하강 반복문(1씩 감소)
 
downSum = 0
 
for j = 1001-do
   downSum = downSum + j
end
cs

1
2
3
4
repeat
   print('down')
   num = num - 1
until num == 0
cs


3. 조건문

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
if num > 40 then
 
    print('over 40')
 
elseif s ~= 'walternate' then (~= 는 같지 않으면을 의미)
 
    io.write('not over 40\n') (stdout)
 
else
 
    thisIsGlobal = 5
 
    local line = io.read (stdin, line 읽기)
 
    print('Hello, ' .. line)
 
end
 
//--------------------------------------------------------------------------------------//
 
aBoolValue = false
 
if not aBoolValue then 
   print('twas false'
end
 
//--------------------------------------------------------------------------------------//
 
ans = aBoolValue and 'yes' or 'no'  -->  no
 
cs


4. 함수

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
//--------------------------------------------------------------------------------------//
// fib함수는 변수로 숫자를 받아 피보나치 수열 출력
 
function fib(n)
   if n < then return end
   return fib(n - 2+ fib(n - 1)
end
 
//--------------------------------------------------------------------------------------//
// 함수 안에 함수를 생성 가능
 
function adder(x)
   return function(y) return x + y end
end
 
a1 = adder(9)
a2 = adder(36)
print(a1(16)) -> 결과 = 25출력
print(a2(64)) -> 결과 = 100출력
cs

1
2
3
4
5
6
7
8
9
10
x, y, z = 1234
 
// a, b, c를 입력받아 6개 숫자를 리턴하는 함수 선언
function bar(a, b, c)
     print(a, b, c)
     return 4815162342
end
 
x, y = bar('zaphod')  -> prints "zaphod nil nil"
 
cs


5. 테이블

1
2
3
4
// 사전(dictionary) 생성
= {key1 = "value1", key2 = false}
 
print(t.key1) -> 'value1'출력
cs

1
2
3
4
5
= {"value1""value2"1.21"endList"}
 
for i = 1, #v do (#v는 리스트 사이즈)
   print(v[i])
end
cs


6. 기타


 PDF 파일 다운로드

 luarefv51.pdf

   


좀더 자세한 내용은


참조


'프로그래밍' 카테고리의 다른 글

ubuntu 14.04 LTS iTorch 설치  (2) 2015.04.18
ubuntu 14.04 LTS Torch 설치  (0) 2015.04.17
안드로이드 가상 머신 지니모션(Genymotion) 설치  (0) 2015.04.09
멀티 스레드  (0) 2014.12.13
파이썬 argparse  (0) 2014.11.08

HoG Descriptor (Histogram of Oriented Gradients descriptors)

컴퓨터비전/영상처리 2015. 4. 15. 02:04

HOG descriptors


이전 포스팅 에서 영상의 기술자에대해서 정리하였습니다.

이번에는 HoG(Histograms of Oriented Gradients)에 정리해봅니다.


HoG는 대표 기술자로 잘 알려져 있는 SIFT[1], SURF[2], GLOH[3]등 다양한 곳에서 사용하고 

있습니다. 그중 가장 친숙한?(혹은 잘 알려진) SIFT에서 사용되는 HoG에 대해 알아보겠습니다.


SIFT는 1999년 캐나다 대학 UBC(University of British Columbia) 교수인 David Lowe가 발표한 논문에서

소개되었습니다. SIFT에서는 두개의 Keypoint(관심점)과 Descriptor(기술자)를 사용하여 두 사진을 매칭 하였습니다.


특징점 추출은 건너 뛰고 기술자 생성에 관한 내용만 정리해 보겠습니다.


첫번째로 SIFT에서는 Scale에 불변한 Keypoint를 찾은후, keypoint위치를 warp 하여 16x16 사이즈의 윈도우를 

만들었습니다.


다음으로 각 픽셀에 대해 기울기를 계산합니다.(Orientation and magnitude) 

각 픽셀에 대하여 아래 사진과 같이 16개의 4x4픽셀 사각형에 포함된 영역에 대한 기울기를 계산합니다.



각각의 사각형에 대해 8방향의 히스토그램을 계산합니다.


위와 같은 방법으로 히스토그램을 구하면 16개의 히스토그램을 8방향으로 표현한 128차원(16*8)의 

벡터를 얻을 수 있습니다.



요약하자면 아래 사진과 같습니다.


위와 같이 히스로그램의 방향을 통해 기술자로 사용하는 방법이 HoG입니다.

SIFT는 illumination(조명)의 변화에 불변(강건) 하며, 물체의 방향에 상관없는 기술자(회전에 불변)를 

생성할 수 있습니다.

References :

[1] Lowe, David G. “Object recognition from local scale-invariant features.”Computer vision, 1999. The proceedings of the seventh IEEE international conference on. Vol. 2. Ieee, 1999.

[2] Bay, Herbert, Tinne Tuytelaars, and Luc Van Gool. “Surf: Speeded up robust features.” Computer VisionECCV 2006. Springer Berlin Heidelberg, 2006. 404-417.

[3] Mikolajczyk, Krystian, and Cordelia Schmid. “A performance evaluation of local descriptors.” Pattern Analysis and Machine Intelligence, IEEE Transactions on27.10 (2005): 1615-1630.



descriptor(기술자)란 무엇인가?

컴퓨터비전/영상처리 2015. 4. 15. 00:28

사진(영상) 에서의 descriptor(기술자)란?


만약 아래 두 사진을 비교한다고 할때 가장 쉽게 생각할 수 있는 방법이 무엇일까요?

아래 두 사진이 같은 사진인지? 아니면 다른 사진인지를 판단하는 방법은 같은 위치에 있는

두 영역(patch)가 모두 같으면 같은 사진일것이고, 다르면 다른 사진일 것 입니다.


       


위와 같이 가장 간단하면서 쉬운 방법은 영상에서의 다양한 위치에 존재하는 patch를 비교하면 

비교 문제를 해결 할 수 있습니다.


영상처리 및 컴퓨터비전에서도 마찬가지로 다양한 부분을 비교하여 많은 부분이 일치하면

동일한 사진이라고 판단하는 방법이 물체의 위치를 찾는데 사용되고 있습니다.


과정은 총 4가지로 이루어져 있습니다.


 1. 두 사진에서 distinctive한 keypoint를 찾습니다. (코너, 에지 등)

 2. 두 사진의 keypoint의 주변을 각각 비교하여 일치 여부를 판단합니다.

 3. 주변 위치가 일치한 keypoint 위치를 쌍으로 homography(호모그래피)를 계산합니다.

 4. 첫번째 이미지의 위치를 두번째 이미지 위치로 조정합니다.



위 사진에서 기술자 함수란 사진에서 특정 영역에 대한 부분을 서로 비교하기 위해

동일한 방법을 통해 사진에서의 특징을 하나의 비교 대상으로 만드는 것을 말합니다.


기술자 함수는 단순하게 영역을 비교하기 위해 설계되어야 하지만 간단한 문제같지만 매우 복잡한

문제일 수도 있습니다. 사진이 회전되었을 경우 혹은 잡음이 섞여 있을 경우, 조명이나 외곡이 있을 경우에도

강건한(robust)한 특징 추출이 가능해야 합니다.


따라서 회전 및 잡음에도 강건한 기술자를 만들기 위해 비교할 좌표(keypoint) 선정에 있어서

코너를 선택하는 경우가 많으며 코너로 부터 기술자를 생성하여 두 이미지를 비교합니다.


pdf ppt 변환

카테고리 없음 2015. 4. 14. 22:10

PDF를 PPT(파워포인트), 워드(word)로 변환하는 방법을 알아봅시다.



홈페이지 이동 : http://smallpdf.com/kr 에서 다양한 변환을 지원합니다.


지원하는 목록을 살펴보면

1. PDF 압축

2. JPG PDF 변환

3. PDF JPG 변환

4. PDF 워드 변환

5. PDF 엑셀 변환

6. PDF PPT 변환

7. 워드 PDF 변환

8. PDF 합치기

9. PDF 분할        등 다양한 기능을 지원합니다.



그중 가장 많이 사용하는 PDF PPT 변환을 해보겠습니다.



PDF PPT변환을 클릭하면 드레그 앤 드랍으로 파일을 선택할 수도 있고, 파일 선택을 누르면 파일 탐색기를 통한

파일 입력이 가능합니다. 또한 Dropbox 및 구글 드라이브 에서 파일을 불러올 수도 있습니다.



파일 선택을 누르면 파일 탐색기가 열리고 변환할 파일을 하나 선택 한 후 열기 버튼을 누르면 업로드가 시작됩니다.



파일 업로드가 시작되면 업로드 중라는 메세지와 함께 변환이 시작됩니다.



업로드 및 변환이 완료되면 파일 다운로드 버튼을 눌러 다운받을 수 있습니다.



youtube 동영상 다운로드

카테고리 없음 2015. 4. 13. 01:39

Youtube에서 동영상 다운로드 하는 방법


Youtube 1080p이상을 받기 위해서는 아래 프로그램이 필요하다.


4k Video Downloader : http://goo.gl/X2pxh3



1. 다운받아 설치해 준다.








2. 설치를 완료했으면 크랙(crack.exe)을 실행시켜 준다.




3. 프로그램 실행



동영상을 다운로드 받으려면 유투브 동영상 주소가 필요하다.


4. 유투브 동영상 주소 복사하기



4.1 위 사진과 같이 유투브 동영상 시청 페이지에 들어가서



4.2 공유 버튼을 누르면 해당 영상에 대한 주소를 가져올수 있다.

     주소를 드레그 한뒤 복사한다(Ctrl + C)



4.3 프로그램으로 돌아와서 Past Link를 클릭하면 Parsing이 시작된다.



4.4 Parsing이 완료되면 다운로드 옵션을 선택할 수 있는 창이 뜬다.



다운로드 받을 동영상 옵션을 선택한뒤 Download버튼을 누르면 다운로드가 시작된다.



안드로이드 가상 머신 지니모션(Genymotion) 설치

프로그래밍 2015. 4. 9. 23:18

상당히 빠른 속도, 다양한 버전의 안드로이드 가상머신을 지원하는 

안드로이드 가상머신 Genymotion을 설치하는법을 소개합니다.



1. https://www.genymotion.com/ 사이트를 방문하면 아래와 같은 화면을 볼수 있습니다.


2. 화면 중앙에 있는 Get Genymotion을 클릭한뒤, 다운로드 버튼을 눌러줍니다.




3. 저는 윈도우개발 환경이므로 운영체제로 윈도우를 선택하였습니다.

   안드로이드 스튜디오, 이클립스와 연동시켜 안드로이드 앱 개발시 가상머신(virtual machine)으로 사용할 수 있습니다.




4. 다운로드를 눌러 설치파일을 받아서 실행시킨 후 다음다음다다음... 을 클릭하면 설치가 완료됩니다.








5. 가상머신을 돌리기 위해서는 Virtual Box가 필요해 또 설치해 줍니다.







6. 설치 완료!



'프로그래밍' 카테고리의 다른 글

ubuntu 14.04 LTS Torch 설치  (0) 2015.04.17
Lua 프로그래밍  (0) 2015.04.16
멀티 스레드  (0) 2014.12.13
파이썬 argparse  (0) 2014.11.08
[STL] List  (0) 2014.03.06

칼만 필터(Kalman Filter)

컴퓨터비전/영상처리 2015. 4. 6. 01:01

칼만 필터는 부정확한 측정값(관찰, 혹은 예측값)으로부터 오차를 최소화 하는 추정치를 반복적으로 

계산하는 방법이다. 가우시안 잡음(Gaussian noise)인 경우에 최적의 추정량을 구할 수 있다.


Wan, Eric A., and Rudolph Van Der Merwe. "The unscented Kalman filter for nonlinear estimation." 

Adaptive Systems for Signal Processing, Communications, and Control Symposium 2000. AS-SPCC. 

The IEEE 2000. IEEE, 2000.


칼만 필터 프로세스 방정식과 측정 방정식은 다음과 같다.





칼만 필터는 초기화 이후 예측(predict) 위 단계를 계속 반복하면서 추정치를 구한다.

측정 잡음의 공분산 R이 아주 크면, 칼만이득 K가 작아져서 다음 단계의 추정치를 계산할 때 현재의 측정값이 무시된다.


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
#include <cv.h>
#include <highgui.h>
#include <iostream>
 
using namespace std;
using namespace cv;
 
int main()
{
    // Greq Welch and Gary Bishop, "An Introduction to the Kalman Filter", 2006.
    // Estimating a Random Constant
 
    IplImage *dstImage = cvCreateImage(cvSize(512512), IPL_DEPTH_8U, 3);
    cvSet(dstImage, CV_RGB(255255255));
 
    IplImage *pImage = cvCreateImage(cvSize(512512), IPL_DEPTH_8U, 3);
    cvSet(pImage, CV_RGB(255255255));
 
    cvNamedWindow("dstImage");
    cvNamedWindow("pImage");
 
    CvRNG rng_state = cvRNG(0xffffffff);
 
    int t, count = 50;
    double x = -0.37727;
 
    int step = dstImage->width / count;
 
    CvMat *= cvCreateMat(count, 1, CV_32FC1);
    cvRandArr(&rng_state, z, CV_RAND_NORMAL, cvRealScalar(x), cvRealScalar(0.1));
 
    ////////////////////////////////
    //        kalman filter       //
    ////////////////////////////////
    double Q = 1e-5// process variance
    double R = (0.1* (0.1); // estimate of measurement variance
    //double R = 1.0; // estimate of measurement variance
    //double R = 0.0001; // estimate of measurement variance
 
    CvMat *xhat = cvCreateMat(count, 1, CV_32FC1); // posteriori estimate of x
    CvMat *= cvCreateMat(count, 1, CV_32FC1); // posteriori error estimate 
    CvMat *xhatminus = cvCreateMat(count, 1, CV_32FC1); // priori estimate of x
    CvMat *Pminus = cvCreateMat(count, 1, CV_32FC1); // priori error estimate
    CvMat *= cvCreateMat(count, 1, CV_32FC1); // kalman gain
 
    // to access as an array in 1D matrix
    float *xhatA = xhat->data.fl;
    float *PA = P->data.fl;
    float *xhatminusA = xhatminus->data.fl;
    float *PminusA = Pminus->data.fl;
    float *KA = K->data.fl;
    float *zA = z->data.fl;
 
    xhatA[0= 0.0;
    PA[0= 1.0;
 
    for (t = 1; t < count; t++)
    {
        //predict
        xhatminusA[t] = xhatA[t - 1];
        PminusA[t] = PA[t - 1+ Q;
 
        //update
        KA[t] = PminusA[t] / (PminusA[t] + R);
        xhatA[t] = xhatminusA[t] + KA[t] * (zA[t] - xhatminusA[t]);
        PA[t] = (- KA[t]) * PminusA[t];
    }
 
 
    //drawing values
    double minVal, maxVal;
    cvMinMaxLoc(z, &minVal, &maxVal);
    double scale = dstImage->height / (maxVal - minVal);
    printf("observations, z : minVal = %f, maxVal = %f\n", minVal, maxVal);
 
    //drawing the truth value, x = -0.37727
    CvPoint pt1, pt2;
    pt1.x = 0;
    pt1.y = dstImage->height - cvRound(scale * x - scale*minVal);
 
    pt2.x = dstImage->width;
    pt2.y = dstImage->height - cvRound(scale * x - scale*minVal);
    cvLine(dstImage, pt1, pt2, CV_RGB(00255), 2);
 
    //drawing the noisy measurements, observations, z
    for (t = 0; t < count; t++)
    {
        pt1.x = t*step;
        pt1.y = dstImage->height - cvRound(scale * zA[t] - scale * minVal);
        cvCircle(dstImage, pt1, 3, CV_RGB(25500), 2);
    }
 
    //drawing the filter estimate, xhat
    pt1.x = 0;
    pt1.y = dstImage->height - cvRound(scale * xhatA[0- scale * minVal);
 
    for (t = 1; t < count; t++)
    {
        pt2.x = t*step;
        pt2.y = dstImage->height - cvRound(scale * xhatA[t] - scale * minVal);
        cvLine(dstImage, pt1, pt2, CV_RGB(02550), 2);
        pt1 = pt2;
    }
 
    //drawing the error covariance, P
 
    cvMinMaxLoc(P, &minVal, &maxVal);
    scale = pImage->height / (maxVal - minVal);
    printf("error covariance, P : minVal = %f, maxVal = %f\n", minVal, maxVal);
 
    pt1.x = 0;
    pt1.y = pImage->height - cvRound(scale * PA[0- scale * minVal);
 
    for (t = 1; t < count; t++)
    {
        pt2.x = t*step;
        pt2.y = pImage->height - cvRound(scale * PA[t] - scale*minVal);
        //printf("PA[%d] = %f\n", t, PA[t]);
        cvLine(pImage, pt1, pt2, CV_RGB(25500), 2);
        pt1 = pt2;
    }
 
    cvShowImage("dstImage", dstImage);
    cvShowImage("pImage", pImage);
 
    cvSaveImage("dstImage.jpg", dstImage);
    cvSaveImage("pImage.jpg", pImage);
    
    cvWaitKey(0);
 
    cvReleaseMat(&z);
    cvReleaseMat(&xhat);
    cvReleaseMat(&P);
    cvReleaseMat(&xhatminus);
    cvReleaseMat(&Pminus);
    cvReleaseMat(&K);
 
    cvReleaseImage(&dstImage);
    cvReleaseImage(&pImage);
}
cs