Collaborative filtering

협업 필터링이라고도 부른다.   행렬 분해(matrix decomposition)를 활용한 상품 추천 시스템(recommendation system)에 사용된다.  또는 아주 큰 데이터(big data)에 내재되어 있는 내부 성질을 파악하기 위한 마이닝(mining) 알고리즘으로도 사용될 수 있다.


웹 Market에서 상품을 팔면서 구매자의 만족도를 조사하였다.

6명의 구매자에 대해 5개의 상품을 0~5점으로 나누어 만족도가 높을수록 5점에 가까운 점수를
주게 하였더니 아래와 같은 행렬이 얻어졌다.

>>A=[2 3 5 0 0; 1 2 4 0 0; 4 5 4 0 0; 0 1 0 3 4; 1 0 0 4 5; 0 0 0 5 4];
>> A

A =
     2     3     5     0     0
     1     2     4     0     0
     4     5     4     0     0
     0     1     0     3     4
     1     0     0     4     5
     0     0     0     5     4

행렬에서 행(row)은 구매자이고 열(column)은 개별 상품에 해당한다.

행렬을 보고 대략 추정해보면 1~3번 상품에 대해 호의적인 사용자들과
4~6번의 상품에 호의적인 사용자들의 두 그룹으로 나누어진다는 것을 알 수 있다.


A행렬을 특이값 분해(SVD, singular value decomposition)해보면

>>[u,d,v]=svd(A)

u =
   -0.4719    0.3337   -0.4308    0.0716    0.1067   -0.6810
   -0.3382    0.2434   -0.5548    0.1188   -0.0594    0.7078
   -0.5811    0.3908    0.6720   -0.1899   -0.0353    0.1436
   -0.3004   -0.3804    0.0835    0.4957   -0.7101   -0.0901
   -0.3610   -0.5078    0.1044    0.3629    0.6803    0.0801
   -0.3234   -0.5235   -0.1927   -0.7532   -0.1299   -0.0100

d =
   10.6983         0         0         0         0
         0   10.1159         0         0         0
         0         0    2.4316         0         0
         0         0         0    1.1504         0
         0         0         0         0    0.9895
         0         0         0         0         0

v =
   -0.3709    0.1944    0.5659   -0.1171    0.7005
   -0.4952    0.3026    0.4284   -0.0013   -0.6925
   -0.5643    0.4157   -0.6929    0.0641    0.1564
   -0.3704   -0.5724   -0.1216   -0.7190   -0.0592
   -0.4020   -0.6084    0.0349    0.6821    0.0419

세개의 행렬 u, d, v가 얻어진다.

d행렬 내의 특이값을 보면 2개의 성분이 크고 나머지는 작다.
즉, A행렬의 내재된 성분이 2개가 중요하다는 의미이다.


2개의 성분과 그에 해당하는 u, v행렬의 값만을 이용하여 A행렬을 다시 재생해 보면

>> vt = v';
>> u(:,1:2)*d(1:2,1:2)*vt(1:2,:)

ans =
    2.5285    3.5220    4.2523   -0.0619   -0.0241
    1.8206    2.5374    3.0657   -0.0693   -0.0437
    3.0741    4.2755    5.1518    0.0396    0.0935
    0.4441    0.4272    0.2140    3.3931    3.6334
    0.4338    0.3579    0.0438    4.3709    4.6782
    0.2538    0.1107   -0.2491    4.3124    4.6127

이고 A행렬과 비교해 보면 값이 비슷하다(잘 재생된다).

내재된 성분의 의미를 파악해 보자.

>> u(:,1:2)*d(1:2,1:2)

ans =
   -5.0490    3.3752
   -3.6186    2.4625
   -6.2169    3.9535
   -3.2142   -3.8482
   -3.8621   -5.1373
   -3.4597   -5.2956

행렬의 행은 사용자를 표현하고, 열은 두개의 내재된 성질에 대한 것이다. 즉, 사용자는 2개의 그룹으로 나누어지며
1~3번 사용자는 첫번째 성질에 속하고, 4~6번 사용자는 두번째 성질에 속하는 정도가 크다 (행렬 요소값의 절대치 크기만 보자).

>> vt(1:2,:)

ans =
   -0.3709   -0.4952   -0.5643   -0.3704   -0.4020
    0.1944    0.3026    0.4157   -0.5724   -0.6084

5개의 상품(각 열)에 대해 첫번째 성질(사용자 그룹)은 1~3열의 상품에 대한 선호도가 높고, 2번째 성질은 4~5열의 상품에 대한 선호도가 높다.


이상과 같이 행렬을 분해해 주는 방법인 SVD를 이용하면 데이터에 내재되어 있는 의미 파악이 가능하다. 이를 이용하면 상품 추천시스템을 만들 수 있다.




이제 어떤 사람의 구매 만족도를 조사하였더니 아래와 같은 데이터가 주어졌다.

B=[3 ? 4 ? 2]

이 사용자는 상품 2, 4를 아직 구매하지 않았다.
이 사람에게 상품 2와 4 중에서 어떤 것을 추천할 것인지 계산해 보자.

먼저 평가가 있는 세개의 상품에 대해 평균을 구해 '?'위치에 대치하자.
(3+4+2)/3=3이므로 3을 대입한다.

>>B=[3 3 4 3 2];


이 사용자를 A행렬의 마지막 행에 추가하고 SVD를 다시 수행한다.

>> C=[A; B]

C =
     2     3     5     0     0
     1     2     4     0     0
     4     5     4     0     0
     0     1     0     3     4
     1     0     0     4     5
     0     0     0     5     4
     3     3     4     3     2

>> [u,d,v]=svd(C)

u =
   -0.4049    0.3207    0.4210   -0.1967   -0.0883    0.7109   -0.0642
   -0.2900    0.2341    0.5467   -0.2343    0.0295   -0.6759   -0.2227
   -0.4980    0.3749   -0.6820   -0.0475    0.1547   -0.0733   -0.3388
   -0.2424   -0.3885   -0.0945   -0.6518    0.4371    0.0039    0.4068
   -0.2938   -0.5175   -0.1111   -0.1982   -0.7383   -0.0344   -0.2190
   -0.2644   -0.5321    0.1912    0.4180    0.4667    0.1088   -0.4540
   -0.5376    0.0046    0.0290    0.5151   -0.1157   -0.1394    0.6418

d =
   12.6750         0         0         0         0
         0   10.1164         0         0         0
         0         0    2.4328         0         0
         0         0         0    1.4264         0
         0         0         0         0    1.0250
         0         0         0         0         0
         0         0         0         0         0

v =
   -0.3943    0.1850   -0.5604    0.3711   -0.5988
   -0.4844    0.2896   -0.4361   -0.2823    0.6416
   -0.5781    0.4011    0.6907   -0.0352   -0.1635
   -0.3816   -0.5814    0.1295    0.6218    0.3360
   -0.3606   -0.6188   -0.0456   -0.6282   -0.3004



이제 내재된 2개의 성질을 이용하여 원래 행렬을 재생해 보면

>> vt=v';
>> u(:,1:2)*d(1:2,1:2)*vt(1:2,:)

ans =
    2.6239    3.4256    4.2679    0.0722   -0.1568
    1.8875    2.4663    3.0746    0.0254   -0.1403
    3.1904    4.1556    5.1695    0.2037   -0.0706
    0.4845    0.3498    0.1994    3.4574    3.5400
    0.5003    0.2878    0.0530    4.4650    4.5828
    0.3258    0.0642   -0.2220    4.4083    4.5394
    2.6959    3.3143    3.9577    2.5738    2.4292


마지막 사용자의 2, 4번째 상품에 대한 평가치는 3.3143>2.5738이므로 2번째 상품을 추천한다.



References

[1] 한국정보처리학회, Tutorial(빅데이터 분석 및 소셜 컴퓨팅), 2014.

FastSLAM

EKF SLAM의 문제점

- 2차원 평면상에서 $K$개의 landmark 존재.  LM 파라미터: $(x_i, y_i)$

- robot pose: $(x, y,  \phi)$

- 합하면 $2K+3$개의 인자

- Kalman filter를 통해 이것들을 추정하면 mean값을 $2K+3$개를 유지

  Covariance matrix는 $(2K+3)^2$크기가 요구
- LM의 수가 늘어 감에 따라 유지해야 하는 파라메터의 수가 급격히 증가(계산량 증가)

FastSLAM에서의 개선

- Landmark는 로봇의 pose가 주어진다면 상호 독립이라 가정

- 각 landmark에 대해 $2\times2$ covariance matrix

  전체 인자의 수: $K(2+2\times2)+3$

Robot poses: $s^t=s_1, s_2, ...,s_t$
Robot actions: $u^t=u_1,u_2,...,u_t$
Observations: $z^t=z_1,z_2,...,z_t$
Landmarks: $\boldsymbol{\theta}=\theta_1,...,\theta_k$
Correspondence: $n^t=n_1,n_2,...,n_t$
- $n_t=k$는 $z_t$가 $\theta_k$의 관찰 값임을 의미
- landmark는 한번에 하나씩 관찰된다고 가정

그림에서 shade영역은 $t_1\sim$현재까지의 로봇의 경로이다.
만일 로봇의 실제 경로를 안다면(즉 shade영역) LM $\theta_1$과 $\theta_2$는 서로 독립이다.


SLAM 문제
- action model: $p(s_t|u^t,s^{t-1})$
- sensor model: $p(z_t|s_t,\boldsymbol{\theta},n_t)$

SLAM 문제는 $z^t, u^t, n^t$가 주어질 때 $\boldsymbol{\theta},s^t$를 추정하는 문제이다.  Bayes theory를 이용하면

$p(s^t,\boldsymbol{\theta}|z^t,u^t,n^t)=
{{p(\boldsymbol{\theta},s^t,z^t,u^t,n^t)} \over {p(z^t,u^t,n^t)}}={{p(\boldsymbol{\theta}|s^t,z^t,u^t,n^t)p(s^t,z^t,u^t,n^t)} \over {p(z^t,u^t,n^t)}}$
$=p(\boldsymbol{\theta}|s^t,z^t,u^t,n^t)p(s^t|z^t,u^t,n^t)$

landmark의 독립가정$^{주1}$에 따라
$p(s^t,\boldsymbol{\theta}|z^t,u^t,n^t)
 =p(s^t|z^t,u^t,n^t)\prod_k^{}p(\theta_k|s^t,z^t,u^t,n^t)$

- trajectory 추정 $p(s^t|z^t,u^t,n^t)$
- robot의 자세 추정은 particle filter를 이용하고, 각 particle에 대한 landmark위치와 분산의 갱신의 EKF를 사용
- particle당 $K$개의 Kalman filter가 있으므로 particle의 수가 $M$개라면 갱신과정의 시간 복잡도는 $O(MK)$이다.



(주1) 독립가정에 따라
$p(\boldsymbol{\theta}|s^t,z^t,u^t,n^t)=...$



References
[1] M. Montemerlo, FastSLAM: A Factored Solution to the Simultaneous Localization and Mapping Problem With Unknown Data Association, CMU-CS-03, Ph.D. thesis, 2003.
[2] Benjamin Kuipers, FastSlam, lecture notes, CS 344R/393R: Robotics.

Design Pattern

작성 중 ...


Adaptor 패턴
두 클래스의 인터페이스가 서로 맞지 않을 때 클래스들을 연결하는 역할을 한다.

어떤 클래스(client)가 다른 클래스의 객체를 받아 이 객체의 함수를 호출해서 사용한다고 하자.  그런데 이 client는 이미 완성된 클래스라 수정이 불가하다.

그런데, client가 이미 만들어져 있는 또 다른 클래스(adaptee)의 함수를 사용할 필요가 생겼다.  adaptee 클래스도 이미 완성되어 수정을 할 수 없다.

이때 사용하는 것이 adaptor이다. client가 넘겨받는 객체 클래스(adaptor)에 adaptee 객체를 선언하고 client가 호출하는 adaptor 함수 내에서 adaptee의 함수를 호출하도록 해 준다.

adaptor 클래스는 adaptee클래스를 상속 받아 만들 수도 있다.
결론적으로 adaptee와 client를 수정할 필요 없이 원하는 목적을 달성할 수 있다.

Adaptor 패턴이 잘 설명
a

(예제) testDuck이라는 함수에서 duck객체를 호출한다고 가정하자.

이때 어느 부분도 수정없이 testDuck에서 Turkey라는 객체도 호출(사용)하고 싶다.

(Sol) duck을 상속받은 TurkeyAdaptor를 만든다. 

TurkeyAdaptor 멤버로 타겟인 Turkey 객체를 삽입한다.

TurkeyAdaptor 내의 duck가상함수를 구현할때 Turkey의 함수를 호출하게 한다. 

(사용) testDuck 함수의 인자에 duck 대신 TurkeyAdaptor객체를 넘겨준다. 

Mediator(중계자) 패턴
a, b, c


중개자, 조정자 패턴이라고도 하는 Mediator 패턴이다.

각 클래스끼리의 상호 참조나, 서로가 서로에게 영향을 미치는 클래스들을 잘 분리하고 중개해주는 역할을 한다.
각 클래스끼리의 영향이 적으면 적을 수록, 나중에 모듈을 떼어내기도 수월하고 디버깅도 유리하다.


Mediator 패턴을 쓰면 좋은 경우를 예로 들면,
GUI 화면을 구성할 때, 특정 체크박스의 값에 따라 각 버튼의 Enabled/Disabled 속성 변경, 텍스트 박스의 값에 따라 각 컴포넌트의 속성이 변경되거나 하는 등 다양한 상태에 따라 각 컴포넌트들이 서로가 서로에게 영향을 미치는 경우를 예로 들 수 있다.

이 때, 각 버튼 클릭 이벤트나, 체크 박스 이벤트, 텍스트 박스 이벤트에서 각 컴포넌트들에게 바로 접근해서 직접 다룰 경우 서로간의 종속성이 엄청 높아지게 된다.
중간에 코드에 문제가 생겨서 수정해야 할 때 고쳐야 할 부분이 분산되어 있어서 복잡해지기도 하며, 특정 컴포넌트를 삭제하거나 변경할 때 코드 내의 여기저기서 컴파일 오류가 발생하는 등 다양한 문제들이 발생하게 된다.


Mediator 패턴은 크게 Mediator 인터페이스와 Colleague 인터페이스로 구성된다.
물론, Mediator과 Colleague는 인터페이스(추상,부모 클래스)이기 때문에 이를 상속받고 구현받는 클래스가 추가로 필요하다.


GUI 예에서 각 텍스트 박스나 버튼, 체크 박스 등은 Colleague 인터페이스를 상속받아서 구현한다.
그리고 각각의 객체들의 값이 변경될 때마다 Mediator의 colleagueChanged() 메소드를 호출해서 변경을 알려준다.

그러면 Mediator에서는 각 Colleague들의 속성을 체크해서 각각의 객체들에게 속성 변경 명령을 내려주게 된다.


Mediator 패턴의 요지는 각각의 클래스 A, B, C 들이 A ↔ B, B ↔ C, C ↔ A 이런 식으로 서로 통신하거나 영향을 미칠 경우
그 관계를 모두 끊고, 중간에 Mediator을 둬서 한 번 거쳐서 처리하도록 하는 것이다.

단계는 하나 더 증가하지만, 각 모듈간의 종속성을 끊고 보다 깔끔하고 유지보수가 쉽게 하자는 것이 그 목적이다.


#include <iostream>
#include <list>
#include <string>

using namespace std;

//------------------------------------------------------------------
// Mediator 인터페이스
class Colleague; // 전방 선언
class Mediator
{
public:
virtual void AppendUser(Colleague* colleague) = 0;
virtual void RemoveUser(Colleague* colleague) = 0;
virtual void sendMessage(const char* message, Colleague* sender) = 0;
};

//------------------------------------------------------------------
// Colleague 인터페이스
class Colleague
{
public:
Colleague(Mediator* m, const char* name) : pMediator(m), mName(name) {}

public:
virtual void SendMessages(const char* str) = 0;
virtual void ReceiveMessages(const char* str) = 0;

protected:
Mediator* pMediator;
string mName;
};

//------------------------------------------------------------------
// User 상속 클래스
class User : public Colleague
{
public:
User(Mediator* m, const char* name) : Colleague(m, name) {}

public:

 // 객체들의 값이 변경될 때마다 Mediator의 colleagueChanged()  
 // 메소드를 호출해서 변경을 알려줌
void SendMessages(const char* str) override
{
cout << mName << " send : " << str << endl;
pMediator->sendMessage(str, this);
}

void ReceiveMessages(const char* str) override
{
cout << mName << " recv : " << str << endl;
}
};

//------------------------------------------------------------------
// ChatMediator 상속 클래스
class ChatMediator : public Mediator
{
public:
void AppendUser(Colleague* colleague) override
{
mList.push_back(colleague);
}

void RemoveUser(Colleague* colleague) override
{
mList.remove(colleague);
}

 // Mediator에서는 각 Colleague들의 속성을 체크해서 각각의 객체들에게 
 // 속성 변경 명령을 내려주게 됨
void sendMessage(const char* message, Colleague* sender)
{
for (Colleague* object : mList)
{
if (object != sender)
object->ReceiveMessages(message);
}
}

private:
list<Colleague*> mList;
};

//------------------------------------------------------------------
// Main
int main()
{
ChatMediator mChatMediator;

User mUser1(&mChatMediator, "홍길동");
User mUser2(&mChatMediator, "나이스");
User mUser3(&mChatMediator, "디자인");

mChatMediator.AppendUser(&mUser1);
mChatMediator.AppendUser(&mUser2);
mChatMediator.AppendUser(&mUser3);

mUser1.SendMessages("안녕하세요. 홍길동입니다!");

return 0;
}

실행 결과:
홍길동 send : 안녕하세요. 홍길동입니다!
나이스 recv : 안녕하세요. 홍길동입니다!
디자인 recv : 안녕하세요. 홍길동입니다!

Template(템플릿) 패턴

부모 클래스에서 알고리즘의 골격을 정의한다. 알고리즘의 여러 단계 중 일부는 서브(자식)클래스에서 구현한다. 

알고리즘의 구조는 그대로 유지하면서 서브클래스에서 특정 단계를 재정의할 수 있다.

이 패턴은 알고리즘의 틀을 만들기 위한 것이다. 상위 클래스에는 일련의 단계들로 알고리즘을 정의한 메소드가 있다.  여러 단계 가운데 하나 이상이 추상 메소드로 정의되며, 그 추상 메소드는 하위클래스에서 구현된다.


(예제) 머신비전의 일반적인 과정을 보면 두개의 쓰레드가 필요하다. 

하나는 처리할 데이터인 영상정보를 카메라에서 캡춰하는 것이다.

다른 하나는 작업 목적에 따라 영상을 처리하는 알고리즘 부분이다. 

두 쓰레드는 동시에 실행되어야 하며 서로 공통 부분이 많다. 

상위 클래스는 쓰레드의 동작에 필요한 공통적인 실행절차를 정의하고 두 쓰레드의 상세한 부분은 하위 클래스에서 재 정의 한다.

Command 패턴

작업을 요청한 쪽과 처리하는 쪽을 분리한다.
연속된 명령의 삽입, 삭제, undo, do 등을 처리할수 있다.

단, 호출자에서 하나의 객체로 처리해야 하므로 모든 명령은 하나의 추상(부모)클래스를 이용한다.

작업 처리 쪽:

작업 요청 쪽:

(몇 유사 패턴과의 비교)

스테이트(State) 패턴은 상태 그 자체를 클래스화해서 사용하는 것이고,
스트래터지(Strategy) 패턴은 알고리즘 자체를 클래스화해서 사용하는 것이고,
컴포지트(Composite) 패턴은 각 객체들을 동일시해서 사용하는 것이다.

커맨드 패턴은 명령 그 자체를 클래스화해서 사용하는 것이다. 디자인 패턴들은 대부분 그 개념면에서 비슷하다.


커맨드 패턴에서 명령들을 전부 동일시해서 사용하기 위해서는 역시 하나의 인터페이스를 이용해서 구현을 해야한다.
 




Composite 패턴

추상클래스를 하나 만들고 상속받은 다양한 자식클래스를 만들어 다양한 자식클래스들을 동일 클래스 다루듯 사용하는 패턴이다.

컴포지트 패턴이란 객체들의 관계를 트리 구조로 구성하여 부분-전체 계층을 표현하는 패턴으로, 사용자가 단일 객체와 복합 객체 모두 동일하게 다루도록 한다.

#include <iostream>
#include <vector>
#include <string>
 
using std::cout;
using std::vector;
using std::string;
 
class Component
{
  public:
   virtual void list() const = 0;
   virtual ~Component(){};
};
 
class Leaf : public Component 
{
  public:
   explicit Leaf(int val) : value_(val) 
   {
   }
   void list() const 
   {
      cout << "   " << value_ << "\n";
   }
  private:
   int value_;
};
 
class Composite : public Component 
{
  public:
   explicit Composite(string id) : id_(id) 
   {
   }
   void add(Component *obj)
   {
      table_.push_back(obj);
   }
   void list() const
   {
      cout << id_ << ":" << "\n";
      for (vector<Component*>::const_iterator it = table_.begin(); 
       it != table_.end(); ++it)
      {
         (*it)->list();
      }
   }
  private:
   vector <Component*> table_;
   string id_;
};
 
int main()
{
   Leaf num0(0);
   Leaf num1(1);
   Leaf num2(2);
   Leaf num3(3);
   Leaf num4(4);
   Composite container1("Container 1");
   Composite container2("Container 2");
 
   container1.add(&num0);
   container1.add(&num1);
 
   container2.add(&num2);
   container2.add(&num3);
   container2.add(&num4);
 
   container1.add(&container2);
   container1.list();
   return 0;
}

Leaf이나 Composite모두 Component에서 상속 받아 만든 동일 구상 클래스이다. 
container1은 2개의 Leaf노드를 가진다.
container2는 3개의 Leaf노드를 가진다.

container1은 container2도 가지므로 container1 밑에 2개의 Leaf와 3개의 Leaf를 가진 container2를 가진 Tree구조가 형성된다.

Composite클래스의 list명령으로 Tree를 모두 출력해 볼 수 있다.



팩토리메쏘드 패턴
객체 생성은 추상클래스에서 한다.
어떤 객체를 생성할지는 구상클래스에서 결정한다.

즉, 인터페이스에 객체 생성을 요청하지만, 어떤 클래스의 인스턴스를 만들지는 서브클래스에서 결정한다. 

//------------------------------------------------------------------
// Product 인터페이스 클래스
class Product
{
public:
  virtual void print() = 0;
};

//------------------------------------------------------------------
// Product 상속 클래스
class ConcreteProduct : public Product
{
public:
  void print() override { cout << "ConcreteProduct" << endl; }
};

//------------------------------------------------------------------
// Creator 클래스 (구현 인터페이스 클래스)
class Creator
{
public:
  Product* AnOperation() { return FactoryMethod(); }

protected:
  virtual Product* FactoryMethod() = 0;
};

//------------------------------------------------------------------
// Creator 상속 클래스 (실제 객체 생성 전담)
class ConcreteCreator : public Creator
{
private:
  Product* FactoryMethod() { return new ConcreteProduct; }
};


//------------------------------------------------------------------
// Main
void main()
{
  ConcreteCreator pCreator;

  Product* pProduct = pCreator.AnOperation();
  pProduct->print();

  delete pProduct;
}



References
Factory Method Pattern: 1, 2

2D Disjunctive Normal Form (DNF)

기계학습 알고리즘을 물체 추적에 적용해 본다.

물체 추적 문제는 보통 비디오의 첫번째 프레임에서 추적할 대상 물체를 사용자가 사각형 box(ROI, 즉 Region Of Interest라고도 부른다)를 그려 지정해 준다.  이렇게 물체 영역을 지정하면, box 내부의 영역은 학습을 위한 positive 영역, 물체 주위의 배경 영역은 negative 영역으로 간주된다.  이 정보를 이용하여 물체 추적을 위한 분류기를 학습할 수 있다.

먼저, 물체를 포함하는 영역과 배경(물체 영역을 중시으로 물체 영역의 2배 정도) 영역에서 랜덤하게 복수 개의 positive patch$^{주1}$와 negative patch를 추출한다. 이러한 patch들은 분류기(classifier)를 학습할 수 있는 샘플(training samples)이 된다.


1st frame
추적할 대상을 지정해 주는 프레임이다.
사용자의 ROI 지정 $\rightarrow$ ROI 내의 랜덤한 위치에서 정해진 갯수의 patch를 추출한다.   각각의 patch들이 training sample이 된다.

학습
분류기의 학습과 갱신은 adaptive boosting알고리즘을 이용한다.
(임의로 선택된) 두 특징 함수 $h_i$와 $h_j$(특징 함수는 Harr-like feature, LBP, HOG 등이 될 수 있다)에 patch를 입력하면 특징 값 $v_i, v_j$가 계산된다.  모든 샘플에 대해 얻어진 $v$의 가능한 값의 범위를 여러 bin(일정한 간격을 가진 여러 값의 범위)으로 나눈다.  이때, 어떤 하나의 샘플이 주는 $h_i$와 $h_j$ 값은 2차원 공간 상의 하나의 bin(특징 값이 2개이므로 2차원이다)을 가리킨다. 또 각각의 sample은 부류정보(positive or negative)를 가지고 있으므로 이 부류 정보에 따라 bin에 Positive나 Negative로 보팅한다$^{주2}$.


모든 sample에 대한 voting후 각 bin에서 positive로 보팅한 수가 negative로 보팅한 수보다 $r$만큼 더 많다면 이 bin 영역의 위치는 저장된다. $r$은 미리 지정된 상수이다.

특징은 일반적으로 여러 개이고 특징 들의 쌍(pair)은 아주 많이 구성할 수 있다.
가능한 여러 pair중에서 정, 부의 모든 샘플을 고려했을 때 가장 작은 오류를 주는 특징 pair를 선정한다.  이것이 하나의 weak DNF classifier이다.

이러한 과정을 $M$번 반복한다. 모두 $M$개의 weak 2d DNF classifier를 선정하면 학습이 완료되고, 선정된 weak classifier를 선형 결합하여 strong classifier를 만든다.


2nd frame
물체의 2배 확대 영역 내 모든 patch에 대해 strong classifier를 적용하고 confidence map를 작성한다.  confidence map의 measure는 strong classifier의 값으로 한다.
계산 방법은 각 patch를 2d DNF weak classifier에 입력하고 계산된 값이 학습 시에 저장된 bin에 포함된다면 이 patch는 positive가 된다.

confidence값을 integral image로 합하고 최대의 confidence값을 가진 영역이 추적 위치가 된다.
아래는 c.m(confidence measure)값을 평가하여 현재 프레임에서 물체가 발견되었는지(FOUND) 또는 잃어버렸는지(LOST)에 대한 모드 평가를 보여준다.

c.m < $th_1$ $\rightarrow$ LOST $\rightarrow$ initial state (전 영역 검색)
c.m > $th_2$ & LOST state $\rightarrow$ FOUND state로 초기화
c.m > $th_2$ & FOUND state $\rightarrow$ Update classifier

실제 구현 시의 strong classifier는 2D DNF만 사용하는 것이 아니라 1차원 weak classifier와 2d DNF classifier의 조합으로 구성된다.


(주1) patch는 8x8 pixels과 같은 특징을 계산할 수 있는 작은 크기의 window이다

(주2) Harr-like 특징을 예를 들어 생각해보면, patch 내 어떤 위치에, 어떤 크기로 두 특징 값 $h_i$와 $h_j$가 정의된다고 가정하자. 하나의 샘플이 선정되고 두 특징에 적용되면 2개의 값이 얻어지고 이것은 두 특징이 정의하는 2차원 평면 상의 좌표 값을 가리키게 된다.  이 좌표 평면은 균등한 간격으로 가로 세로로 나누어진 많은 작은 영역(bin)들이 정의되어 있다면 이 값은 이 영역들 중의 하나로 떨어지게 된다.  voting이란 두 특징 값이 가리키는 영역의 counter(초기에는 0) 값을 하나 증가시키는 것이 된다.  counter는 샘플의 부류에 따라 positive, negative의 두 개가 존재한다.



참고 논문
[1] Disjunctive Normal Form of Weak Classifiers for Online Learning based Object Tracking, ZhuTeng et. al., VISAPP'2013.

Markov Chain Monte Carlo

Monte Carlo
컴퓨터의 연산능력을 이용하여 수식만으로 계산하기 어려운 문제가 있을 때, 데이터의 무작위 샘플을 얻고 이를 이용하여 해답을 구하는 방법이다.
예를 들면, 원의 원주율 $\pi$를 구한다고 하자. 샘플링만으로 이 값을 구하기 위하여, 정사각형 내에 반지름이 1인 내접하는 원을 그린 후, 무작위 샘플을 발생시킨다. 그리고 이 샘플이 원 내부에 있는지, 외부인지 판단 만으로 값을 구할 수 있다. 즉, $2^2 - \pi \cdot 1^2 = $ samples between two areas 로 구할 수 있다.
이때 고성능 컴퓨터를 이용하여 샘플을 많이 발생 시킬수록 더 정확한 값을 구할 수 있다.

Markov Chain
연속으로 이어지는 일련의 값이 주어질 때, 다음에 나오는 값은 바로 직전(이전) 값에만 영향을 받는 연속 값의 체인이다.  이 때 값은 상태(state)를 말하고 체인은 상태 값의 연속 시퀀스이다.
값들은 상호 독립적이지 않고 서로 상관성이 존재하는데 연속된 값 사이에서만 그 상관성이 존재하고 멀리 떨어진 값들 사이의 상관성은 무시된다.

MCMC
어떤 수식의 값을 샘플링 만으로 근사하여 구하며, 이때 샘플을 발생시키는 방법을 상호 독립적이지 않고 어떤 관계를 가지도록 발생 시킨다는 의미로 사용된다.  따라서 이름이 MCMC 근사로 붙여졌다.
고차공간에서 확률분포 추정에 사용하며 왜 MCMC인지는 아래 그림 하단부(MH 알고리즘)에서 다시 설명한다.



MCMC(Markov Chain Monte Carlo)는 고차원 공간에서 수행해야 하는 적분이나 최적화 문제를 풀기 위해 사용된다. 주요 적용 분야는 기계학습, 물리, 통계학, 경제학, 결정이론 등이다.
고차공간에서 확률분포 $p(z)$의 추정은 매우 어렵다.  확률분포(probability density function: pdf)는

$p(z)={\tilde{p}(z) \over Z_p}$, where $Z_p = \int_Z {\tilde{p}(z) dz}$

와 같이 계산한다.  여기서 $\tilde{p}(z)$는 정규화가 안된 함수 값이고, $Z_p$는 정규화 상수(normalization constant)이다.  보통 $\tilde{p}(z)$는 쉽게 구할 수 있지만, $Z_p$의 경우는 $Z$ 공간 내의 가능한 모든 영역의 $z$에 대해 적분을 필요로 하며, 적분할 수식의 형태도 주어지지 않는 경우가 많아 어렵다.  이 때 sampling기법을 이용하는 MCMC는 $p(z)$의 추정에 유용한 도구이다.

MCMC가 사용되는 몇 가지 분야를 살펴 보자.

(1) Bayesian inference and learning: 미지의 변수 $x \in X$와 데이터 $y \in Y$가 주어질 때, 다루기 힘든(intractable) 적분의 존재는 Bayesian 통계학에서 자주 나타난다.

$\bullet$ 정규화: Bayes' 이론에서 분포의 정규화 항은 고차원 공간 $X$에 대한 적분을 필요로 한다.

$p(x|y)={p(y|x) \cdot p(x) \over p(y)}={p(y|x) \cdot p(x) \over \int_{X}p(y|x')p(x')dx'}$

$p(y)$를 직접 추정하든지, 주변확률을 이용하여 추정해야 한다.  어느 경우나, 정규화 상수 $p(y)$의 추정은 쉽지 않다.  
분모 $p(y)$ 값 계산을 주변 확률을 이용하려 하려면 적분 인자인 변수 $x'$가 무엇이고 어떤이 있는지 알아야 하는데 불가능하다.

$\bullet$ 주변화(Marginalization): 합이나 적분을 통해 어떤 변수를 제거할 수 있다. 이 때, 공간 $Z$가 크다면 적분이 어려워 진다.

$p(x|y)=\int_{Z} p(x,z|y)dz$

$\bullet$ 기대값(Expectation): 기대값(또는 가중 평균)이 필요할 때 적분을 사용할 수 있다.  공간이 크면 적분이 어렵다.

$E_{p(x|y)}(f(x))=\int_{X}f(x)p(x|y)dx$

  
(2) Statistical mechanics(주1): 상태 $s$와 Hamiltonian $E(s)$를 가진 어떤 시스템의 partition 함수 $Z$를 계산할 때 필요하다.

$Z=\sum_{S} exp[-{E(s) \over kT}]$,

여기서

$k$는 Boltzmann상수, $T$는 시스템의 Temperature를 나타낸다. 상태 $S$의 공간이 아주 크다면 합의 계산은 풀기 어려운 문제가 된다.

(3) 최적화(optimization): 최적화는 무수히 많은 해답 공간 내에서 목적함수를 최소화하는 후보를 하나 골라 내는 것이다. 그런데 이 공간은 연속이고 경계가 없다.  답을 찾기 위해 모든 공간을 뒤지는 것은 어렵다.

(4) 우도(likelihood) 모델의 선정: 두 단계로 구성된다. 먼저, 가능한 각각의 모델에 대해 ML(maximum likelihood) 추정치를 찾는다. 다음, 그 중 하나를 선정하기 위해 비용(penalization) 항(MDL, BIC, or AIC)을 사용한다.  이 방법의 문제점은 모델 집합이 아주 클 경우이다.



Monte Carlo의 예
MC의 아이디어는 어떤 고차 공간에서 정의된 확률 분포 $p(x)$로부터 상호 독립적인 샘플들의 집합 $\{ x^{(i)} \}_{i=1}^N$을 추출하는 것이다. 고차 공간은 시스템의 다양한 구성을 나타내는 공간이거나, 사후확률(posterior)이 정의된 공간, 가능한 해답의 결합 공간이 될 수 있다.
샘플들은 공간의 확률 분포를 근사적으로 표현하기위해 사용된다.

$p_{N}(x)={1 \over N} \sum_{i=1}^N {\delta_{x^{(i)}}(x)}$

여기서 $\delta_{x^{(i)}}(x)$는 Dirac-delta함수로 $x$가 $x^{(i)}$에 위치할 때 1이 되고, 그렇지 않으면 0이 된다.  따라서 $x$가 공간의 확률 분포에 따라 샘플링된다면 $p_{N}(x)$는 공간의 $pdf$에 대한 히스토그램(histogram)이 된다.


(예제) $x \in \{1,2,3 \}, N=100, pdf=Uniform$일 때  $p_{N}(2)$를 구하시오.

$p_{N}(2)={1 \over 100} \sum_{i=1}^N {\delta_{x^{(i)}}(2)}$

이고 이것은  $x^{(i)}=2$인 샘플의 비율이다. 즉, $x=2$인 확률이고, $p_{N}(2) \approx {33 \over 100}$이다.
만일 $pdf=N(2,3^2)$라면 $P_N(2)$는 ${1 \over 3}$보다 휠씬 클 것이다. $x^{(i)}$는 $pdf$에 따라 샘플링되므로 $p_{N}(x)$가 근사적으로 그 $pdf$이다.


Monte Carlo를 이용하면 다루기 힘든(intractable) 함수의 적분도 $I_{N}(f)$로 근사할 수 있다.

$x$의 $pdf$에 따른 $f(x)$의 근사 적분을 이산 합을 통해 구성해 보면

$I_{N}(f)={1 \over N} \sum_{i=1}^N { f(x^{(i)})}$

가 된다.
이 때 $x$의 공간을 이산에서 연속으로 확장하면 $N$은 $\infty$가 되고 이 식은

$I(f)=\int_{X}f(x)p(x)dx$

이다.


Metropolis Hasting 알고리즘
MH는 대표적인 MCMC 알고리즘이다. 샘플링을 임으로 하지 않고 Markov 성질을 가지는 연속 값의 체인을 통해 수행한다.

위 그림을 이해하기 위해 그림의 내부에 표시된 1->2->3 순으로 읽어본다.  MH을 사용하면 최소한의 sampling으로 분포를 근사화 할 수 있다.


그래프와 코드를 보고 간략히 설명한다. 그림에서 2개의 봉우리를 가진 어떤 분포가 있다. 이 분포를 데이터 샘플링을 통해 근사하고자 한다.

분포 $q$는 proposal 확률이라 부르며, 현재 위치 $x$에서 다음 위치 $x^*$을 샘플링하는 것에 관여한다. 연속 샘플링인 마코프 체인의 이동은 acceptance 확률 $A$에 따라 $x$에서 $x^*$로 움직이는데, reject되면 $x$에 그냥 머문다.
$q$가 Gaussian proposal이라면 $q(x^*|x^{(i)})=N(x^{(i)},\sigma)$이므로 이동은 대칭형 random walk이고 $q(x^*|x^{(i)})=q(x^{(i)}|x^*)$이다.

알고리즘을 보면 먼저, $u$를 하나 샘플링하는데 이 값은 0~1사이의 uniform분포에서 얻는다. 다음 $x^*$를 $q$분포에서 하나 샘플링한다.
분포 $q$는 여러가지 형태가 될 수 있으나 여기서는 Gaussian과 같은 대칭형 분포로 $q(x^*|x^{(i)})=q(x^{(i)}|x^*)$ 관계가 성립하므로 수식의 분자, 분모에 있는 $q$는 소거되어 샘플링의 기준치는 $A=min( 1, {p({x^*}) \over {p(x^{(i)})}} )$가 된다.
$x^*$가 $q(x^*|x^{(i)})$에 따라 샘플링된다는 말은 이전 샘플링 값인 $x^{(i)}$근처에서 샘플링된다고 볼 수 있다.  Gaussian 분포의 분산 크기에 따라 샘플링 점과 현재 점 사이의 거리가 결정된다. 타겟 분포 $p$는 모르지만 위치 $x$에서 확률 값 $p(x)$는 측정 가능하다고 가정한다.

어쨌던 현재의 $u$는 이렇게 계산된 $A$보다 작을 수도 있고 클 수도 있다.
만일 $x^*$에서의 확률 값 $p(x^*)$가 $p(x^{(i)})$보다 크다면 $A$는 1이되고 $u<A$이므로 $x^*$는 무조건 accept된다.  즉, 이전 샘플 $x^{(i)}$근처에서 샘플링된 $x^*$의 확률값이 $x^{(i)}$보다 크다면 이 샘플은 분포가 커지는 방향에 있고 accept되어야 좋다. 이는 확률값이 높은 곳에서 더 많은 샘플이 얻어질 수 있다는 것을 표현하게 된다.
 
만일 $x^*$에서의 확률 값 $p(x^*)$가 $p(x^{(i)})$보다 작다면 $A$는 1보다 작아진다. 이 때는 $u$의 값에 따라 $u<A$가 되거나 $u>A$가 될 수도 있다. $u$의 값에 따라 $x^*$가 받아질 수도 있고 아닐 수도 있다.

이러한 방식으로 아주 많은 샘플을 취하면 이 샘플들의 빈도는 확률값의 크기에 따라 얻어지게 되고 높은 확률 값에서 많은 샘플이 얻어지므로 확률분포를 잘 근사하게 된다.

여기서 Monte Carlo는 샘플링으로 어떤 값을 근사한다는 의미로 사용되었고, Markov Chain은 샘플링 시, 현재 샘플링 값은 이전 샘플링 값과 독립적이지 않고 어떤 관계를 가지면서 얻어진다는 의미로 사용된다. 따라서 MCMC 근사가 된다.

위 그림은 분포 $q$의 분산 크기가 미치는 영향을 보여 준다. 분산 값에 따라 MCMC근사의 정밀도가 달라진다.





(주1) Statistical mechanics is a branch of mathematical physics that studies, using probability theory, the average behavior of a mechanical system where the state of the system is uncertain. 


References
[1] http://arongdari.tistory.com/m/post/view/id/62
[2] C. Andrieu, An Introduction to MCMC for Machine Learning, Machine Learning, 2003.

Digital Image Correlation

인장이나 압축력을 받는 시편의 변형 전후의 영상을 DIC로 비교한다. 변형전 영상을 $f$, 변형 후의 영상을 $g$라 했을 때, 영상 계측을 통해 변형을 측정한다. 픽셀 좌표의 변형 전후 관계를 표현하기 위해 좌표의 선형 변환을 나타내는 6개의 변수와 함께 픽셀에서의 영상 밝기 값이 전체적으로 밝아지거나 어두워지는 global offset을 나타내기 위해서 dc값 $w$를 사용한다.   어떤 화소 위치 $x_1$가 선형 변환 $T$에 의해 새로운 위치 $x_2$로 바뀐다면 그 변환 관계는 다음과 같다.

$x_2=T(x_1)$ 이고 $T$가 Affine 변환이라면 가장 일반적으로 사용하는 변형 식은

$\begin{pmatrix} x_2 \\ y_2 \end{pmatrix}=\begin{pmatrix}a & b \\ c & d \end{pmatrix} \begin{pmatrix}x_1 \\ y_1 \end{pmatrix}+\begin{pmatrix} e \\ f \end{pmatrix}$이고 변환에 관계하는 파라메터의 수는 6개이다.

(a)

(b)

Figure: Local part of image for deformation calculation. (a) Dotted box shows local part of image for deformation test. We assume that the small window undergoes linear deformation that is described by six parameters; (b) From left side, four small images show deformation surfaces at same location on 1st, 100th, 200th, 300th sequential frame of specimen under tensile strength, respectively. Image views show that rightward deformation is strong.

기계요소의 표면 변형 계산은 먼저 변형 전 표면에 ROI를 설정하고 이 영역 내의 화소들이 변형 후 위치가 어떻게 변환 되었는지를 따지는 것이다.
이때 변형 전후의 두 window patch 영상을 비교하기 위한 기준 치가 필요하며 식 (1)을 사용한다. global offset 변수 $w$는 scalar 값이다.

$C=\sum_{G_s \in S} {[f(G_s)-\tilde{g}(G_s, P)]^2 \over \sum_{G_s\in S}{f(G_s)^2} }$                  (1)

where

$\tilde{g}(G_s, P)=g(T(G_s))-w$

이다. $P$를 파라메터 벡터, $G_s$를 변형전 영상의 픽셀 좌표라 하자. $S$는 window patch(ROI) 내의 픽셀 집합이다. 위 그림에서 노란 점선 박스로 표시되어 있다.
$T(G_s)$에 의한 화소 위치 $G_s$의 선형 변환은 변형 후 sub-pixel이 되므로 보간이 필요하다.  즉 , 변형 후 영상 $g$를 사용하여 $T(G_s)$ 위치에서의 화소 값을 구해야한다.

Bilinear 보간은 직선과 직선이 만나는 위치에서 불연속이 있으므로 $P$를 구하기 위한 최적화식의 수렴 특성에 문제점이 있어 bicubic 보간을 주로 사용한다.

$P$내의 최적화 변수는 7개(6개+$w$)가 되며, 최적화 방법으로는 Newton법을 사용한다.

Newton법은 Taylor series 전개를 통해 유도되며 2차 편미분이 사용된다. 1차 편미분만 사용하는 gradient descent법보다 효율적이다[2]. 그러나 2차 편미분의 연산은 복잡하고 계산량이 과도한 문제가 있다.  Newton 식을 얻기 위해 식 (1)을 $P_0$ 근방에서 Taylor 전개하자.

$C(P)=C(P_0)+\triangledown C(P_0)^T \delta P+\left. 1 \over 2 \right. \delta P^T\triangledown\triangledown C(P_0)\delta P$.                             (2)

where

$\delta P = P-P_0$

이다. 여기서 $P$는 $P_0$ 근방에서 $C(P)$의 최적값(최소값)을 주는 파라메터 벡터이다.  식 (2)의 양변을 $\delta P$에 대해 미분하면

$\triangledown \triangledown C(P_0) \delta P=-\triangledown C(P_0)$

이고

$\delta P=-{\triangledown \triangledown C(P_0)}^{-1} \cdot \triangledown C(P_0)$

가 되고,

$P=P_0-{\triangledown \triangledown C(P_0)}^{-1} \cdot \triangledown C(P_0)$                      (3)

가 된다.  식 (3)이 최적화를 위한 Newton 반복식이다.


위에서 $\delta P$에 대한 식 (2)의 미분의 의미를 생각해 보기 위해 식(2)를 다시 쓰면,

$C(P_0+\delta P)-C(P_0) = f(\delta P)=\triangledown C(P_0)^T \delta P+...$

이므로 $f'(\delta P)$는

$\lim_{\delta P \to 0}{C(P+\delta P)-C(P) \over \delta P} = C'(P)=\triangledown C(P) +\triangledown \triangledown C(P) \delta P$

이다.  따라서 $C'(P)=0$을 구하면 $C(P)$가 최소가 되는 $P$를 구하는 것이 된다. 즉, 두 영상의 유사도가 높도록 $P$를 구해 내는 것이 된다.



지금 풀어야 하는 식 (3)의 문제는 내부에 2차 편미분이 들어 있는 것이다.  이차 미분식을 다시 쓰면,

$\triangledown \triangledown C(P)=\left( {\partial^2 C \over \partial P_i \partial P_j} \right )_{i=1..7, j=1..7}$

이 값을 구해보기 위해 식 (1)을 $P_i$에 대해 2차 미분을 계산하면

${\partial^2 C \over \partial P_i \partial P_j}$
$=-{2 \over \sum_{G_i \in S}{f^2(G_s)}}$  $\sum_{G_s \in S}{[f(G_s)-\tilde{g}(G_s,P)]}$
$\cdot {\partial^2 \tilde{g}(G_s,P) \over \partial P_i \partial P_j}$
$+{2 \over \sum_{G_i \in S}{f^2(G_s)}} {\partial \tilde{g}(G_s,P) \over \partial P_i}{\partial \tilde{g}(G_s,P) \over \partial P_j}$

이다. 만일 $P$가 exact solution에 충분히 가깝다면

$\tilde{g}(G_s, P) \approx f(G_s)$

이므로(즉, 변형 전후의 영상은 근본적으로 같은 부분을 찍은 것이므로 서로 같다.)

${\partial^2 C \over \partial P_i \partial P_j}$
$ \approx {2 \over \sum_{G_i \in S}{f^2(G_s)}} {\partial \tilde{g}(G_s,P) \over \partial P_i}{\partial \tilde{g}(G_s,P) \over \partial P_j}$

이고 Hessian 행렬은 계산하기 쉬운 두 1차 미분의 곱으로 근사된다. 

Bicubic 스플라인 보간을 사용하므로

$\tilde{g}(G_s, P)=g(\tilde{x}, \tilde{y})-w=\alpha_{mn} \tilde{x}^n \tilde{y}^m-w$            (4)

where $m,n=0,1,2,3$

이다. DIC에서 사용하는 선형 변환 관계식(Affine 변형식)의 일반적 형태는

$\tilde{x}=x+P_1+P_3(x-x_0)+P_5(y-y_0)$
$\tilde{y}=x+P_2+P_4(x-x_0)+P_6(y-y_0)$
 
이다. 여기서 $P=(u_0,v_0,u_x,u_y,v_x,v_y,w)^T$로 표기하여 사용한다.
Newton식의 계산에는 식 (4)와 이 식에 대한 7개 변수의 도함수가 필요하다.  1차 도함수들의 예는

${\partial \tilde{g}(G_s, P) \over \partial P_1}={\partial \tilde{g}(G_s, P) \over \partial \tilde{x}} {\partial \tilde{x} \over \partial P_1}
+{\partial \tilde{g}(G_s, P) \over \partial \tilde{y}} {\partial \tilde{y} \over \partial P_1}={\partial \tilde{g}(G_s, P) \over \partial \tilde{x}}$

${\partial \tilde{g}(G_s, P) \over \partial P_2}={\partial \tilde{g}(G_s, P) \over \partial \tilde{y}}$

${\partial \tilde{g}(G_s, P) \over \partial P_7}=1$

${\partial \tilde{g}(G_s, P) \over \partial P_3}=(x-x_0){\partial \tilde{g}(G_s, P) \over \partial \tilde{x}}$

${\partial \tilde{g}(G_s, P) \over \partial P_4}=(y-y_0){\partial \tilde{g}(G_s, P) \over \partial \tilde{y}}$

${\partial \tilde{g}(G_s, P) \over \partial P_6}=(x-x_0){\partial \tilde{g}(G_s, P) \over \partial \tilde{y}}$

${\partial \tilde{g}(G_s, P) \over \partial P_5}=(y-y_0){\partial \tilde{g}(G_s, P) \over \partial \tilde{x}}$

이다.

Newton법의 유도가 해의 초기 추측치 $P_0$ 근방에서 Taylor전개로 유도 되었으므로 초기치 $P_0$를 구해야 한다.  일단 변위의 도함수들인 $u_x, u_y, v_x, v_y$와 global offset변수인 $w$를 0으로 놓고 탐색 영역 내의 모든 가능한 픽셀 위치에 대해 $u_0, v_0$를 탐색으로 찾아 낸다.  이 때, 가장 작은 계수치 $C$를 주는 위치에서 파라메터를

$(u_0,v_0,0,0,0,0,0)$

로 두고 Newton식을 이용한 최적화를 수행한다.


(추가 작성 예정) DIC의 C++ 구현


References
[1] 위키백과: TeX 문법, Wikipedia: TeX grammar (English)
[2] G. Vendroux and W. G. Knauss, Submicron deformation field measurements: Part 2. Improved digital image correlation, Experimental Mechanics, 38(2), 1998.
[3] Wikipedia Newton's method in optimization page,
      http://en.wikipedia.org/wiki/Newton's_method_in_optimization

PCA LI-tracker의 함수

작성 중 ...


PCA L1 tracker의 메인 함수인 demo.m에서 호출하는 두 주요 함수는

estwarp_condens_PCAL1.m
sklm.m

이다.



function param = estwarp_condens_PCAL1(frm, tmpl, param, opt)
%% function param = estwarp_condens_PCAL1(frm, tmpl, param, opt)
%%      frm:                 Frame image;
%%      tmpl:                PCA model;
%%        -tmpl.mean:           PCA mean vector
%%        -tmpl.basis:          PCA basis vectors
%%        -tmpl.eigval:         The eigenvalues corresponding to basis vectors
%%        -tmpl.numsample:      The number of samples
%%      param:
%%        -param.est:           The estimation of the affine state of the tracked target 
%%        -param.wimg:          The collected sample for update
%%      opt:                 
%%        -opt.numsample:       The number of sampled candidates
%%        -opt.condenssig:      The variance of the Guassian likelihood function
%%        -opt.ff:              Forgotten factor;
%%        -opt.bacthsize:       The number of collected samples for update
%%        -opt.affsig:          The variance of affine parameters
%%        -opt.tmplsize:        The size of warpped image patch
%%        -opt.maxbasis:        The maximum number of basis vectors
%%        -opt.srParam:         The parameters of L1 minimization
%%        -opt.threshold:       The threshold for model update
%%DUT-IIAU-DongWang-2012-05-10
%%Dong Wang, Huchuan Lu, Minghsuan Yang, Online Object Tracking with Sparse
%%Prototypes, IEEE Transaction On Image Processing
%%http://ice.dlut.edu.cn/lu/index.html
%%wangdong.ice@gmail.com
%%

%%************************1.Candidate Sampling************************%%
%%Sampling Number
n = opt.numsample; % 600
%%Data Dimension
sz = size(tmpl.mean); % 32x32
N = sz(1)*sz(2); % 1024

%%Affine Parameter Sampling
% 600개의 affine 파라메터를 랜덤하게 생성
param.param = repmat(affparam2geom(param.est(:)), [1,n]);
% param.est를 600개 복사, size(param.param)=6x600 
% 여기서 주의할 점은 param.est를 그대로 사용안하고 
% affparam2geom으로 변환 후에 600개를 발생한다.
% affparam2geom의 출력은 [tx,ty,ramda1(수직길이),theta, ramda2/ramda1(길이비),phi]이다.
% 따라서 param.est는 [tx,ty,a11,a12,a21,a22]이다. 

randMatrix = randn(6,n); % 6x600, mean(randMatrix,2)=(0.023 0.042 ...-0.036)
param.param = param.param + randMatrix.*repmat(opt.affsig(:),[1,n]);

%%Extract or Warp Samples which are related to above affine parameters
wimgs = warpimg(frm, affparam2mat(param.param), sz);
% size(wimgs)=32x32x600, size(wimgs(:,:,600))=32x32->600개 중 하나의 영상


%%************************1.Candidate Sampling************************%%

%%*******************2.Calucate Likelihood Probablity*******************%%
%%Remove the average vector or Centralizing
diff = repmat(tmpl.mean(:),[1,n]) - reshape(wimgs,[N,n]);
% size(tmpl.mean(:))=1024x1, repmat(..,[1,600])은 이것을 600개의 열로 복사
% 따라서 repmat(...)의 크기는 1024x600, 같은 이미지를 600개의 열에 복사 
% size(reshape(wimgs,[N,n]))=1024x600
% size(diff)=1024x600, 변형에 따라 약간씩 평균과 다른 이미지 600개 

%
if  (size(tmpl.basis,2) > 0)
 
    %%(1)PCA_L1
    if  (size(tmpl.basis,2) == opt.maxbasis) % maxbasis=16
        %(1.1)Calucate representation coefficients of all candidates
        alpha = zeros(N+opt.maxbasis,n); % N(1024), opt.maxbasis(16), n(600)
        % alpha=1040x600
        for num = 1:n
            alpha(:,num) = pca_L1(diff(:,num), tmpl.basis, opt.srParam);
            % size(diff(:,1))=1024x1, tmpl.basis=1024x16,
            % opt.srParam(lambda(0.05), L0(0.05), maxLoopNum(20), tol(10e-3))
        end
        coeff = alpha(1:size(tmpl.basis,2),:);    %%The coefficients of PCA basis vectors
        % size(coeff)=16x600
        err = alpha(size(tmpl.basis,2)+1:end,:);  %%The coefficients of trivial templates
        % size(err) = 1024x600
        %%(1.2)Calucate observation likelihood via Eq.12
        diff = diff-tmpl.basis*coeff-err;         %%Reconstruction, diff=1024x600
        diff = diff.*(abs(err)<opt.srParam.L0);   % abs(err)<opt.srParam.L0 = 1024x600                          
        param.conf = exp(-(sum(diff.^2)
           +opt.srParam.L0*sum(abs(err)>=opt.srParam.L0))./opt.condenssig)';
        % sum(diff.^2)=1x600, sum은 각 열의 합. 
        % 두 항이 다 작아야 되는데, 앞의 항은 err가 작은 항들의 합이므로 가림이 없는 부분
        % 뒤 항은 err가 큰 항들의 합이니 가림 부분
        % 즉, 비 가림부분은 잔차가 작아야 하고, 가림 부분은 err의 합이 작아야(sparse) 해야 함.
    %%(2)Traditional PCA
    % diff는 이미 평균과 임의로 발생된 600개 이미지와의 차이 이미지이다.
    % 그런데 basis를 이용해서 diff를 변환한다. 
    %    
    % basis가 16개가 안되면 현재의 basis로 각 600개의 차 이미지를 표현하고
    % 600개 각각에서 pca로 표현된 이미지를 뺀다. 
    else % basis가 16개가 안될 때 수행
        coef = tmpl.basis'*diff;
        % coff=""5x600"", basis=1024x5, diff=1024x600
        % coff(:,1)->5차원의 축에서 1번째 영상의 좌표 (값 5개로 1번째 영상 표현) 
        % coff(:,2)->5차원의 축에서 2번째 "          (값 5개로 2번째 영상 표현)
        % size(tmpl.basis)=1024x5 (basis축에서 좌표 값으로 diff를 projection)
        % basis는 하나하나가 좌표축의 방향이다. 
        diff = diff - tmpl.basis*coef;
        % tmpl.basis*coef(:,1): 5개의 계수로 1번째 영상을 복원 
        % tmpl.basis*coef(:,2): 5개의 계수로 2번째 영상을 복원
        param.conf = exp(-sum(diff.^2)./opt.condenssig)';
    end
 
else % wimgs가 5개 이상이 되면 basis추출 연산이 수행되므로 wimgs가 5개가 되기 전 단계에서 수행
    %%(3)Square Error
    param.conf = exp(-sum(diff.^2)./opt.condenssig)'; %condenssig=0.25
    % size(diff.^2)=1024x600, size(sum(diff.^2))=1x600 
    % param.conf=600개, 비슷한 것이 큰 값임. 
end
%%*******************2.Calucate Likelihood Probablity*******************%%

%%*****3.Obtain the optimal candidate by MAP (maximum a posteriori)*****%%
param.conf = param.conf ./ sum(param.conf); % pdf를 위한 정규화 
[maxprob,maxidx] = max(param.conf); % maxprob=0.2233(130번째)
param.est = affparam2mat(param.param(:,maxidx));
% param.est에 저장은 [tx,ty,a11,a12,a21,a22]형태로 한다. 

%%*****3.Obtain the optimal candidate by MAP (maximum a posteriori)*****%%

%%************4.Collect samples for model update(Section III.C)***********%%
wimg = wimgs(:,:,maxidx);
if  (size(tmpl.basis,2) == opt.maxbasis)
    err = abs(alpha(size(tmpl.basis,2)+1:end,maxidx)); % err=1024x1
    %Compute the ratio
    errRatio = sum(err>opt.srParam.L0)/length(err); % 7/1024
    %Full update: occlusion이나 background clutter부분이 아주 작아 전체 영역을 갱신
    if  (errRatio < opt.threshold.low) % 0.1
        param.wimg = wimg; % 600개중의 최적의 이미지를 현 단계의 추적 이미지로 할당
        return;
    end
    %No update:
    %Occlusion이나 background clutter부분이 너무 많아 갱신을 하지 않음
    if  (errRatio > opt.threshold.high) % 0.6
        param.wimg = [];
        return;
    end
    %Partial update: 패치 내의 일부분의 영역만 occlusion이나 clutter임
    %이 경우는 가림 영역과 가려지지 않은 영역으로 분리하여
    %가려진 영역은 기존 mean영상의 부분에서 가져오고,
    %가려지지 않은 부분은 현재 프레임의 최적 데이터에서 가져와 합함.
    if  ((errRatio>opt.threshold.low) && (errRatio<opt.threshold.high))
        param.wimg = (1-(err>opt.srParam.L0)).*wimg(:) +
                                      (err>opt.srParam.L0).*tmpl.mean(:);
        % 식의 앞 부분 (1-(err>opt.srParam.L0)): occlusion이 있는 부분은 모두 0으로 
        % 만들고, 0이 아닌 부분만 wimg(:)의 데이터를 취함  
        % 뒷 부분 (err>opt.srParam.L0): 가림이 있는 부분만 취함. 이 부분은 모두 1이므로
        % tmpl.mean에서 데이터를 가져옴
        % 가림이 없는 부분은 현 프레임의 데이터, 있는 부분은 누적 평균 데이터에서 취함
        %  
        % 이 함수의 호출후 메인(demo.m)으로 리턴되면 param.wimg는 누적 됨.
        % 5 프레임까지는 누적 프레임의 평균을 구해 tmpl.mean에 저장(sklm.m에서 mu).
        % 누적 프레임이 5개가 되면 리셋됨. 
        param.wimg = reshape(param.wimg,size(wimg));
        return;
    end
else
    param.wimg = wimgs(:,:,maxidx);
end
%%************4.Collect samples for model update(Section III.C)***********%%