블로그 이미지
Leeway is... the freedom that someone has to take the action they want to or to change their plans.
maetel

Notice

Recent Post

Recent Comment

Recent Trackback

Archive

calendar

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
  • total
  • today
  • yesterday

Category

2010. 6. 22. 18:24 Computer Vision
Material zur Vorlesung "Technische Optimierung'':    
EVOLUTION AND OPTIMUM SEEKING
von     Hans-Paul Schwefel


posted by maetel
2010. 5. 18. 00:26 Computer Vision
ref.
2010/02/10 - [Visual Information Processing Lab] - R. Y. Tsai "A Versatile Camera Calibration Technique for High Accuracy 3-D Maching Vision Metrology Using Off-the-shelf TV Cameras and Lenses"


(1) 고정되어 있는 것으로 가정한 카메라의 내부 파라미터 값들을 구하고 (2) 실시간으로 들어오는 이미지 프레임마다 카메라의 회전과 이동을 계산하기 위하여 Tsai 알고리즘을 쓰기로 하고, C 또는 C++로 구현된 소스코드 또는 라이브러리를 찾아서 붙여 보기로 한다.


Try #1.
처음에는 CMU의 Reg Willson가 C로 짠 Tsai Camera Calibration 코드 에서 필요한 부분을 include하여 쓰려고 했는데, C++ 문법에 맞지 않는 구식 C 문법으로 코딩된 부분이 많아서 고치는 데 애를 먹었다. (Xcode의 C++ 프로젝트에서 .c 파일을 include하면 compile은 되지만, linking error가 난다. 때문에 .c를 .cpp로 바꾸어야 함.)  그런데 결정적으로, "cal_main.cpp" 파일에 정의된, 캘리브레이션의 최종 결과값을 주는 함수들이 호출하는 optimization을 실행하는 함수 lmdif_()가  Fortan 파일 "lmdif.f"에 정의되어 있고, Fortran을 C로 변환해 주는 "f2c.h"에 의해 이것을 "lmdif.c"로 하여 가지고 있다는 문제가 있었다. lmdif.c를 lmdif.cpp 형태로 만들기 위해서는 Fortran 언어와 Fortran을 C++로 변환하는 방법을 알아야 하므로, 결국 포기했다.



Try #2.
Michigan State University Charles B. Owen Display-Relative Calibration (DRC)을 구현한 DRC 프로그램( DRC.zip )에서 카메라 캘리브레이션에 Tsai의 알고리즘 libtsai.zip을 쓰고 있다. 이 라이브러리는 위의 C 코드를 C++로 수정하면서 "CTsai"라는 클래스를 사용하고 여러 함수들을 수정/보완/결합한 것인데, Visual Studio 용 프로젝트 프로그램을 만들면서 Windows 환경에 기반하여 MFC를 활용하였다. 그래서 이것을 나의 Mac OS X 기반 Xcode 프로젝트에서 그대로 가져다 쓸 수는 없다. 용법은 다음과 같다.

DRC/DisplayRelativeCalibration.cpp:
bool CDisplayRelativeCalibration::ComputeCameraCalibration(void)
{
    CTsai tsai;

    tsai.Width(m_camerawid);
    tsai.Height(m_camerahit);

    for(std::list<Corr>::const_iterator i=m_cameracorr.begin();  i!=m_cameracorr.end();  i++)
    {
        tsai.Point(i->x, i->y, i->z, i->u, i->v);
    }

    if(tsai.PointCount() < 8)
        return Error("Didn't get enough points");

    if(!tsai.Compute())
        return Error("Camera calibration failed");

    for(int n=0;  n<tsai.PointCount();  n++)
    {
        double ux, uy;
        tsai.WorldToImage (tsai.PointX(n), tsai.PointY(n), tsai.PointZ(n), ux, uy);

        m_cameraproj.push_back(CGrPoint(ux, uy, 0));
    }

   
    m_cameraf = tsai.F();
    m_cameracx = tsai.Cx();
    m_cameracy = tsai.Cy();
    m_camerakappa1 = tsai.Kappa1();
    m_camerasx = tsai.Sx();
    memcpy(m_cameramatrix, tsai.CameraMatrix(), sizeof(double) * 16);

    return true;
}




문제점#1.

class CTsai 안의 member functions 중에  ncc_compute_exact_f_and_Tz( )와 ncc_compute_exact_f_and_Tz_error( )가 있는데,

libtsai.h:21
class CTsai
{

    bool ncc_compute_exact_f_and_Tz();
    bool ncc_compute_exact_f_and_Tz_error (int m_ptr, int n_ptr, const double *params, double *err);

};

전자인 ncc_compute_exact_f_and_Tz()가 정의된 부분을 보면, 

Tsai_ncc.cpp:274
bool CTsai::ncc_compute_exact_f_and_Tz()
{
    CLmdif<CTsai> lmdif;

    lmdif.Lmdif (this, ncc_compute_exact_f_and_Tz_error,
            m_point_count, NPARAMS, x,
            NULL, NULL, NULL, NULL);
}

클래스 형태의 템플릿( CLmdif )으로 선언된 "lmdif"의 member function "Lmdif"를 호출할 때, 

min/Lmdif.h:48
template<class T> class CLmdif : private CLmdif_
{

int Lmdif(T *p_user, bool (T::*p_func)(int m, int n, const double *parms, double *err),
        int m, int n, double *x, double *fvec, double *diag, int *ipvt, double *qtf)

};

후자인 같은 member function, ncc_compute_exact_f_and_Tz_error()를 인자로 넣고 있고 (위 부분 코드들 중 오렌지 색 부분), 컴파일 하면 이 부분을 <unknown type>으로 인식하지 못 하겠다는 에러 메시지를 보낸다. 그리고 다음과 같은 형태를 추천한다고 한다.
 
note: candidates are: int CLmdif<T>::Lmdif(T*, bool (T::*)(int, int, const double*, double*), int, int, double*, double*, double*, int*, double*) [with T = CTsai]

function pointer의 형태가 틀린 모양인데, 오렌지색 부분을 그냥 함수가 아닌 어떤 class의 non-static member function을 가리키는 pointer로  &CTsai::ncc_compute_exact_f_and_Tz_error 이렇게 바꾸어 주면, 에러 메시지가 다음과 같이 바뀐다.

error: no matching function for call to 'CLmdif<CTsai>::Lmdif(CTsai* const, bool (*)(int, int, const double*, double*), int&, const int&, double [3], NULL, NULL, NULL, NULL)'

연두색 부분 대신 CTsai::ncc_compute_exact_f_and_Tz_error 이렇게 바꾸어 주면, 에러 메시지가 다음과 같다.

error: no matching function for call to 'CLmdif<CTsai>::Lmdif(CTsai* const, bool (&)(int, int, const double*, double*), int&, const int&, double [3], NULL, NULL, NULL, NULL)'

해결:
편법으로, class CLmdif를 클래스 형 템플릿이 아닌 그냥 클래스로 바꾸어서 선언하고 연두색 부분처럼 호출하면 에러는 안 나기에 일단 이렇게 넘어가기로 한다.


문제점#2.
코드에서 Windows OS 기반 MFC를 사용하고 있어 Mac OS X에서 에러가 난다.

해결:
MFC를 사용하는 "StdAfx.h"는 모두 주석 처리한다.


문제점#3.
Lmdif.h



... 기타 등등의 문제점들을 해결하고, 캘리브레이션을 수행한 결과가 맞는지 확인하자.

source code:
           if ( CRimage.size() > 0 ) // if there is a valid point with its cross ratio
            {  
                correspondPoints(indexI, indexW, p, CRimage, linesYorder.size(), linesXorder.size(), world, CRworld, dxList.size(), dyList.size(), iplMatch, scale );
            }  
            cvShowImage( "match", iplMatch );
            cvSaveImage( "match.bmp", iplMatch );
           
            cout << "# of pairs = " << indexI.size() << " = " << indexW.size() << endl;
           
            // # 6. camera calibration
           
            int numPair = indexI.size();
           
            tsai.Clear();
           
            for( int n = 0;  n < numPair;  n++ )
            {
                tsai.Point(world[indexW[n]].x, world[indexW[n]].y, world[indexW[n]].z, p[indexI[n]].x, p[indexI[n]].y);
               
                cout << "pair #" << n << ": " << p[indexI[n]].x << "  " <<  p[indexI[n]].y << "  : "
                    << world[indexW[n]].x << "  " << world[indexW[n]].y << "  " << world[indexW[n]].z << endl;
            }
           
            if( numPair < 8 )
                cout << "Didn't get enough points" << endl;
           
            if(!tsai.Compute())
                cout << "Camera calibration failed" << endl;
           
            cout << endl << "camera parameter" << endl
            << "focus = " << tsai.F() << endl
            << "principal axis (x,y) = " <<  tsai.Cx() << ", " <<  tsai.Cy() << endl
            << "kappa1 (lens distortion) = " <<  tsai.Kappa1() << endl
            << "skew_x = " << tsai.Sx() << endl;
          
            // reproject world points on to the image frame to check the result of computing camera parameters
            for(int n=0;  n<tsai.PointCount();  n++)
            {
                double ux, uy;
                tsai.WorldToImage (tsai.PointX(n), tsai.PointY(n), tsai.PointZ(n), ux, uy);
                CvPoint reproj = cvPoint( cvRound(ux), cvRound(uy) );
                cvCircle( iplInput, reproj, 3, CV_RGB(200,100,200), 2 );
            }
           
// draw a cube on the image coordinate computed by camera parameters according to the world coordinate
            drawcube( tsai, iplInput, patSize );
            cvShowImage( "input", iplInput );  



아래 사진은 구해진 카메라 내부/외부 파라미터들을 가지고 (1) 실제 패턴의 점에 대응하는 이미지 프레임 (image coordinate) 상의 점을 찾아 (reprojection) 보라색 원으로 그리고, (2) 실제 패턴이 있는 좌표 (world coordinate)를 기준으로 한 graphic coordinate에 직육면체 cube를 노란색 선으로 그린 결과이다.

이미지 프레임과 실제 패턴 상의 점을 1 대 1로 비교하여 연결한 16쌍의 대응점

구한 카메라 파라미터를 가지고 실제 패턴 위의 점들을 이미지 프레임에 reproject한 결과 (보라색 점)와 실제 패턴의 좌표를 기준으로 한 그래픽이 이미지 프레임 상에 어떻게 나타나는지 그린 결과 (노란색 상자)

 

위 왼쪽 사진에서 보여지는 16쌍의 대응점들의 좌표값을 "이미지 좌표(x,y) : 패턴 좌표 (x,y,z)"로 출력한 결과:
# of pairs = 16 = 16
pair #0: 7.81919  36.7864  : 119.45  82.8966  0
pair #1: 15.1452  71.2526  : 119.45  108.484  0
pair #2: 26.1296  122.93  : 119.45  147.129  0
pair #3: 36.6362  172.36  : 119.45  182.066  0
pair #4: 77.3832  20.4703  : 159.45  82.8966  0
pair #5: 85.4293  53.7288  : 159.45  108.484  0
pair #6: 97.8451  105.05  : 159.45  147.129  0
pair #7: 109.473  153.115  : 159.45  182.066  0
pair #8: 96.6046  15.962  : 171.309  82.8966  0
pair #9: 105.046  48.8378  : 171.309  108.484  0
pair #10: 118.177  99.9803  : 171.309  147.129  0
pair #11: 130.4  147.586  : 171.309  182.066  0
pair #12: 145.469  4.50092  : 199.965  82.8966  0
pair #13: 154.186  36.5857  : 199.965  108.484  0
pair #14: 168.033  87.5497  : 199.965  147.129  0
pair #15: 180.732  134.288  : 199.965  182.066  0


그런데 위 오른쪽 사진에서 보여지는 결과는 이전 프레임에서 20쌍의 대응점으로부터 구한 카메라 파라미터 값을 가지고 계산한 결과이다.
# of found lines = 8 vertical, 7 horizontal
vertical lines:
horizontal lines:
p.size = 56
CRimage.size = 56

# of pairs = 20 = 20
pair #0: -42.2331  53.2782  : 102.07  108.484  0
pair #1: -22.6307  104.882  : 102.07  147.129  0
pair #2: -4.14939  153.534  : 102.07  182.066  0
pair #3: 1.81771  169.243  : 102.07  193.937  0
pair #4: -10.9062  41.1273  : 119.45  108.484  0
pair #5: 8.69616  92.7309  : 119.45  147.129  0
pair #6: 27.0108  140.945  : 119.45  182.066  0
pair #7: 32.9779  156.653  : 119.45  193.937  0
pair #8: 57.4164  14.6267  : 159.45  108.484  0
pair #9: 77.7374  65.9516  : 159.45  147.129  0
pair #10: 96.3391  112.934  : 159.45  182.066  0
pair #11: 102.524  128.555  : 159.45  193.937  0
pair #12: 76.5236  7.21549  : 171.309  108.484  0
pair #13: 97.5633  58.2616  : 171.309  147.129  0
pair #14: 116.706  104.705  : 171.309  182.066  0
pair #15: 123.108  120.238  : 171.309  193.937  0
pair #16: 125.015  -11.5931  : 199.965  108.484  0
pair #17: 146.055  39.453  : 199.965  147.129  0
pair #18: 164.921  85.2254  : 199.965  182.066  0
pair #19: 171.323  100.758  : 199.965  193.937  0

camera parameter
focus = 3724.66
principal axis (x,y) = 168.216, 66.5731
kappa1 (lens distortion) = -6.19473e-07
skew_x = 1



대응점 연결에 오차가 없으면, 즉, 패턴 인식이 잘 되면, Tsai 알고리즘에 의한 카메라 파라미터 구하기가 제대로 되고 있음을 확인할 수 있다. 하지만, 현재 full optimization (모든 파라미터들에 대해 최적화 과정을 수행하는 것)으로 동작하게 되어 있고, 프레임마다 모든 파라미터들을 새로 구하고 있기 때문에, 속도가 매우 느리다. 시험 삼아 reprojection과 간단한 graphic을 그리는 과정은 속도에 큰 영향이 없지만, 그전에 카메라 캘리브레이션을 하는 데 필요한 계산 시간이 길다. 입력 프레임이 들어오는 시간보다 훨씬 많은 시간이 걸려 실시간 구현이 되지 못 하고 있다.

따라서, (1) 내부 파라미터는 첫 프레임에서 한 번만 계산하고 (2) 이후 매 프레임마다 외부 파라미터 (카메라의 회전과 이동)만을 따로 계산하는 것으로 코드를 수정해야 한다.




Try#3.
OpenCV 함수 이용

1) 내부 파라미터 계산
cvCalib
rateCamera2


2) lens distortion(kappa1, kappa2)을 가지고 rectification
cvInitUndistortRectifyMap

3) line detection

4) 패턴 인식 (대응점 찾기)

5) 외부 파라미터 계산 (4의 결과 & lens distortion = 0 입력)
cvFindExtrinsicCameraParams2

6) reprojection
2)에서 얻은 rectificated image에 할 것


posted by maetel
2010. 4. 27. 22:07

보호되어 있는 글입니다.
내용을 보시려면 비밀번호를 입력하세요.

2010. 4. 22. 20:14

보호되어 있는 글입니다.
내용을 보시려면 비밀번호를 입력하세요.

2010. 2. 10. 15:47 Computer Vision
Seong-Woo Park, Yongduek Seo, Ki-Sang Hong: Real-Time Camera Calibration for Virtual Studio. Real-Time Imaging 6(6): 433-448 (2000)
doi:10.1006/rtim.1999.0199

Seong-Woo Park, Yongduek Seo and Ki-Sang Hong1

Dept. of E.E. POSTECH, San 31, Hyojadong, Namku, Pohang, Kyungbuk, 790-784, Korea


Abstract

In this paper, we present an overall algorithm for real-time camera parameter extraction, which is one of the key elements in implementing virtual studio, and we also present a new method for calculating the lens distortion parameter in real time. In a virtual studio, the motion of a virtual camera generating a graphic studio must follow the motion of the real camera in order to generate a realistic video product. This requires the calculation of camera parameters in real-time by analyzing the positions of feature points in the input video. Towards this goal, we first design a special calibration pattern utilizing the concept of cross-ratio, which makes it easy to extract and identify feature points, so that we can calculate the camera parameters from the visible portion of the pattern in real-time. It is important to consider the lens distortion when zoom lenses are used because it causes nonnegligible errors in the computation of the camera parameters. However, the Tsai algorithm, adopted for camera calibration, calculates the lens distortion through nonlinear optimization in triple parameter space, which is inappropriate for our real-time system. Thus, we propose a new linear method by calculating the lens distortion parameter independently, which can be computed fast enough for our real-time application. We implement the whole algorithm using a Pentium PC and Matrox Genesis boards with five processing nodes in order to obtain the processing rate of 30 frames per second, which is the minimum requirement for TV broadcasting. Experimental results show this system can be used practically for realizing a virtual studio.


전자공학회논문지 제36권 S편 제7호, 1999. 7 
가상스튜디오 구현을 위한 실시간 카메라 추적 ( Real-Time Camera Tracking for Virtual Studio )   
박성우 · 서용덕 · 홍기상 저 pp. 90~103 (14 pages)
http://uci.or.kr/G300-j12265837.v36n07p90

서지링크     한국과학기술정보연구원
가상스튜디오의 구현을 위해서 카메라의 움직임을 실시간으로 알아내는 것이 필수적이다. 기존의 가상스튜디어 구현에 사용되는 기계적인 방법을 이용한 카메라의 움직임 추적하는 방법에서 나타나는 단점들을 해결하기 위해 본 논문에서는 카메라로부터 얻어진 영상을 이용해 컴퓨터비전 기술을 응용하여 실시간으로 카메라변수들을 알아내기 위한 전체적인 알고리듬을 제안하고 실제 구현을 위한 시스템의 구성 방법에 대해 다룬다. 본 연구에서는 실시간 카메라변수 추출을 위해 영상에서 특징점을 자동으로 추출하고 인식하기 위한 방법과, 카메라 캘리브레이션 과정에서 렌즈의 왜곡특성 계산에 따른 계산량 문제를 해결하기 위한 방법을 제안한다.



Practical ways to calculate camera lens distortion for real-time camera calibration
Pattern Recognition, Volume 34, Issue 6, June 2001, Pages 1199-1206
Seong-Woo Park, Ki-Sang Hong




generating virtual studio




Matrox Genesis boards
http://www.matrox.com/imaging/en/support/legacy/

http://en.wikipedia.org/wiki/Virtual_studio
http://en.wikipedia.org/wiki/Chroma_key

camera tracking system : electromechanical / optical
pattern recognition
2D-3D pattern matches
planar pattern


feature extraction -> image-model matching & identification -> camera calibration
: to design the pattern by applying the concept of cross-ratio and to identify the pattern automatically


영상에서 찾아진 특징점을 자동으로 인식하기 위해서는 공간 상의 점들과 영상에 나타난 그것들의 대응점에 대해서 같은 값을 갖는 성질이 필요한데 이것을 기하적 불변량 (Geometric Invariant)이라고 한다. 본 연구에서는 여러 불변량 가운데 cross-ratio를 이용하여 패턴을 제작하고, 영상에서 불변량의 성질을 이용하여 패턴을 자동으로 찾고 인식할 수 있게 하는 방법을 제안한다.


Tsai's algorithm
R. Y. Tsai, A Versatile Camera Calibration Technique for High Accuracy 3-D Maching Vision Metrology Using Off-the-shelf TV Cameras and Lenses. IEEE Journal of Robotics & Automation 3 (1987), pp. 323–344.

direct image mosaic method
Sawhney, H. S. and Kumar, R. 1999. True Multi-Image Alignment and Its Application to Mosaicing and Lens Distortion Correction. IEEE Trans. Pattern Anal. Mach. Intell. 21, 3 (Mar. 1999), 235-243. DOI= http://dx.doi.org/10.1109/34.754589

Lens distortion
Richard Szeliski, Computer Vision: Algorithms and Applications: 2.1.6 Lens distortions & 6.3.5 Radial distortion

radial alignment constraint
"If we presume that the lens has only radial distortion, the direction of a distorted point is the same as the direction of an undistorted point."

cross-ratio  http://en.wikipedia.org/wiki/Cross_ratio
: planar projective geometric invariance
 - "pencil of lines"
http://mathworld.wolfram.com/CrossRatio.html
http://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/MOHR_TRIGGS/node25.html
http://www.cut-the-knot.org/pythagoras/Cross-Ratio.shtml
http://web.science.mq.edu.au/~chris/geometry/


pattern identification

 카메라의 움직임을 알아내기 위해서는 공간상에 인식이 가능한 물체가 있어야 한다. 즉, 어느 위치에서 보더라도 영상에 나타난 특징점을 찾을 수 있고, 공간상의 어느 점에 대응되는 점인지를 알 수 있어야 한다.

패턴이 인식 가능하기 위해서는 카메라가 어느 위치, 어느 자세로 보던지 항상 같은 값을 갖는 기하적 불변량 (Geometric Invariant)이 필요하다.

Coelho, C., Heller, A., Mundy, J. L., Forsyth, D. A., and Zisserman, A.1992. An experimental evaluation of projective invariants. In Geometric invariance in Computer Vision, J. L. Mundy and A. Zisserman, Eds. Mit Press Series Of Artificial Intelligence Series. MIT Press, Cambridge, MA, 87-104.


> initial identification process
extracting the pattern in an image: chromakeying -> gradient filtering: a first-order derivative of Gaussian (DoG) -> line fitting: deriving a distorted line (that is actually a curve) equation -> feature point tracking (using intersection filter)


R1x = 0



http://en.wikipedia.org/wiki/Difference_of_Gaussians



real-time camera parameter extraction

이상적인 렌즈의 optical axis가 영상면에 수직이고 변하지 않는다고 할 때, 영상 중심은 카메라의 줌 동작 동안 고정된 값으로 계산된다. (그러나 실제 렌즈의 불완전한 특성 때문에 카메라의 줌 동작 동안 영상 중심 역시 변하게 되는데, 이 변화량은 적용 범위 이내에서 2픽셀 이하이다. 따라서 본 연구에서는 이러한 변화를 무시하고 이상적인 렌즈를 가정하여 줌동작에 의한 영상 중심을 구하게 된다.)

For zoom lenses, the image centers vary as the camera zooms because the zooming operation is executed by a composite combination of several lenses. However, when we examined the location of the image centers, its standard deviation was about 2 pixels; thus we ignored the effect of the image center change.


calculating lens distortion coefficient

Zoom lenses are zoomed by a complicated combination of several lenses so that the effective focal length and distortion coefficient vary during zooming operations.

When using the coplanar pattern with small depth variation, it turns out that focal length and z-translation cannot be separated exactly and reliably even with small noise.

카메라 변수 추출에 있어서 공간상의 특징점들이 모두 하나의 평면상에 존재할 때는 초점거리와 z 방향으로의 이동이 상호 연관 (coupling)되어 계산값의 안정성이 결여되기 쉽다.


collinearity

Collinearity represents a property when the line in the world coordinate is also shown as a line in the image. This property is not preserved when the lens has a distortion.


Once the lens distortion is calculated, we can execute camera calibration using linear methods.


filtering

가상 스튜디오 구현에 있어서는 시간 지연이 항상 같은 값을 가지게 하는 것이 필수적이므로, 실제 적용에서는 예측 (prediction)이 들어가는 필터링 방법(예를 들면, Kalman filter)은 사용할 수가 없었다.

averaging filter 평균 필터








Orad  http://www.orad.co.il

Evans & Sutherland http://www.es.com









posted by maetel
2009. 7. 3. 10:44 Computer Vision

'Computer Vision' 카테고리의 다른 글

photometric stereo 09-07-03  (0) 2009.07.06
GO MC PS references  (0) 2009.07.03
Linear Least Squares  (0) 2009.07.03
Reinhard Diestel <Graph Theory>  (0) 2009.06.16
photometric stereo 2009-06-10  (0) 2009.06.10
posted by maetel
2009. 4. 16. 15:26 Footmarks
Richard Hartley
http://users.rsise.anu.edu.au/~hartley/

http://users.rsise.anu.edu.au/~jaehak/

http://users.rsise.anu.edu.au/~hongdong/


<Camera Motion Estimation for Multi-Camera Systems>

Urban mapping and construction of 3D urban models has become a popular area of research recently.  The object is to construct models from images taken by a vehicle moving in the urban environment. Typically such a vehicle will take images on both sides simultaneously, with cameras pointing in opposite directions.  Part of the problem is to estimate the motion of the vehicle from a sequence of images taken by the non-overlapping cameras.  I will discuss methods of doing this task, leading to three different algorithms that we have developed for estimating the motion of a rig with non-overlapping cameras.


2009-04-16 나무 @GA705

'Footmarks' 카테고리의 다른 글

RoSEC 2010 winter school  (0) 2010.01.14
CUDA  (0) 2009.04.24
음악의 경계를 넘다: (3) 음악과 소리, 그 경계의 애매함  (0) 2008.05.11
Stephen Lin "Seeing the World as It Really is"  (0) 2008.05.07
라울 자무디오 강연  (0) 2008.03.29
posted by maetel
2009. 3. 31. 20:22 Computer Vision

A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking

Arulampalam, M.S.   Maskell, S.   Gordon, N.   Clapp, T.  
Defence Sci. & Technol. Organ., Adelaide, SA , Australia;


This paper appears in: Signal Processing, IEEE Transactions on [see also Acoustics, Speech, and Signal Processing, IEEE Transactions on]
Publication Date: Feb. 2002
Volume: 50 , Issue: 2
On page(s): 174 - 188
ISSN: 1053-587X
CODEN: ITPRED
INSPEC Accession Number:7173038
Digital Object Identifier: 10.1109/78.978374
Current Version Published: 2002-08-07

posted by maetel

sedumi

Kolman & Beck, Ch.1 Eg.1
: Activity Analysis or Product Mix

>> c = [-120; -100; 0; 0];
>> A = [ 2, 2, 1, 0; 5, 3, 0, 1];
>> b = [8, 15];
>> x= sedumi(A,b,c)
SeDuMi 1.1R3 by AdvOL, 2006 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
eqs m = 2, order n = 5, dim = 5, blocks = 1
nnz(A) = 6 + 0, nnz(ADA) = 4, nnz(L) = 3
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            4.18E-002 0.000
  1 : -5.02E+002 1.25E-002 0.000 0.2989 0.9000 0.9000   1.77  1  1  3.8E+000
  2 : -4.29E+002 3.65E-003 0.000 0.2927 0.9000 0.9000   2.52  1  1  6.4E-001
  3 : -4.29E+002 2.84E-004 0.000 0.0778 0.9900 0.9900   1.19  1  1  4.5E-002
  4 : -4.30E+002 9.99E-007 0.000 0.0035 0.9990 0.9990   1.03  1  1 
iter seconds digits       c*x               b*y
  4      0.2   Inf -4.3000000000e+002 -4.3000000000e+002
|Ax-b| =  2.3e-015, [Ay-c]_+ =  0.0E+000, |x|= 2.9e+000, |y|= 3.6e+001

Detailed timing (sec)
   Pre          IPM          Post
3.125E-002    1.563E-001    0.000E+000   
Max-norms: ||b||=15, ||c|| = 120,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 2.11427.
x =
   (1,1)       1.5000
   (2,1)       2.5000
>>








Kolman & Beck, Ch.1 Eg.3
: The Transportation Problem

>> c = [5; 7; 9; 6; 7; 10; 0; 0; 0; 0; 0];
>> A = [1,1,1,0,0,0,1,0,0,0,0; 0,0,0,1,1,1,0,1,0,0,0; 1,0,0,1,0,0,0,0,-1,0,0; 0,1,0,0,1,0,0,0,0,-1,0; 0,0,1,0,0,1,0,0,0,0,-1];
>> b = [120; 140; 100; 60; 80];
>> x = sedumi(A,b,c)
SeDuMi 1.1R3 by AdvOL, 2006 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
eqs m = 5, order n = 12, dim = 12, blocks = 1
nnz(A) = 17 + 0, nnz(ADA) = 17, nnz(L) = 12
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            3.24E+003 0.000
  1 :  1.16E+003 9.69E+002 0.000 0.2991 0.9000 0.9000   2.37  1  1  1.6E+000
  2 :  1.63E+003 2.05E+002 0.000 0.2114 0.9000 0.9000   1.34  1  1  3.0E-001
  3 :  1.69E+003 3.84E+001 0.000 0.1875 0.9000 0.9000   1.20  1  1  5.1E-002
  4 :  1.70E+003 1.08E+000 0.000 0.0282 0.9900 0.9900   1.03  1  1 
iter seconds digits       c*x               b*y
  4      0.5  15.9  1.7000000000e+003  1.7000000000e+003
|Ax-b| =  5.6e-014, [Ay-c]_+ =  7.6E-016, |x|= 1.1e+002, |y|= 1.4e+001

Detailed timing (sec)
   Pre          IPM          Post
3.125E-001    4.531E-001    1.250E-001   
Max-norms: ||b||=140, ||c|| = 10,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 1.
x =
   (1,1)      68.6083
   (3,1)      51.3917
   (4,1)      31.3917
   (5,1)      60.0000
   (6,1)      28.6083
   (8,1)      20.0000







posted by maetel
SeDuMi - Home
optimization package over symmetric cones

SeDuMi = Self-Dual-Minimization/ Optimization over self-dual homogeneous cones
: MatLab Toolbox developed by Jos F. Sturm


1) SeDuMi를 다운로드 받기 위해 홈페이지에 회원 등록
2) 메뉴의 Download에서 적당한 버전 다운로드
3) Matlab을 실행시키고, 메뉴에서 File > Set Path > Add Folder 클릭, SeDuMi가 들어 있는 폴더 지정 후 저장
4) command window에서 'help sedumi' 입력하여 toolbox 추가 확인


>> help sedumi
                                           [x,y,info] = sedumi(A,b,c,K,pars)
   
SEDUMI  Self-Dual-Minimization/ Optimization over self-dual homogeneous
          cones.
 
  >  X = SEDUMI(A,b,c) yields an optimal solution to the linear program
       MINIMIZE c'*x SUCH THAT A*x = b, x >= 0
       x is a vector of decision variables.
       If size(A,2)==length(b), then it solves the linear program
       MINIMIZE c'*x SUCH THAT A'*x = b, x >= 0
 
  >  [X,Y,INFO] = SEDUMI(A,b,c) also yields a vector of dual multipliers Y,
       and a structure INFO, with the fields INFO.pinf, INFO.dinf and
       INFO.numerr.
 
     (1) INFO.pinf=INFO.dinf=0: x is an optimal solution (as above)
       and y certifies optimality, viz.\ b'*y = c'*x and c - A'*y >= 0.
       Stated otherwise, y is an optimal solution to
       MAXIMIZE b'*y SUCH THAT c-A'*y >= 0.
       If size(A,2)==length(b), then y solves the linear program
       MAXIMIZE b'*y SUCH THAT c-A*y >= 0.
 
     (2) INFO.pinf=1: there cannot be x>=0 with A*x=b, and this is certified
       by y, viz. b'*y > 0 and A'*y <= 0. Thus y is a Farkas solution.
 
     (3) INFO.dinf=1: there cannot be y such that c-A'*y >= 0, and this is
       certified by x, viz. c'*x <0, A*x = 0, x >= 0. Thus x is a Farkas
       solution.
 
     (I)   INFO.numerr = 0: desired accuracy achieved (see PARS.eps).
     (II)  INFO.numerr = 1: numerical problems warning. Results are accurate
           merely to the level of PARS.bigeps.
     (III) INFO.numerr = 2: complete failure due to numerical problems.
 
     INFO.feasratio is the final value of the feasibility indicator. This
     indicator converges to 1 for problems with a complementary solution, and
     to -1 for strongly infeasible problems. If feasratio in somewhere in
     between, the problem may be nasty (e.g. the optimum is not attained),
     if the problem is NOT purely linear (see below). Otherwise, the reason
     must lie in numerical problems: try to rescale the problem.
 
  >  [X,Y,INFO] = SEDUMI(A,b,0) or SEDUMI(A,b) solves the feasibility problem
     FIND x>=0 such that A*x = b
 
  >  [X,Y,INFO] = SEDUMI(A,0,c) or SEDUMI(A,c) solves the feasibility problem
     FIND y such that A'*y <= c
 
  >  [X,Y,INFO] = SEDUMI(A,b,c,K) instead of the constraint "x>=0", this
       restricts x to a self-dual homogeneous cone that you describe in the
       structure K. Up to 5 fields can be used, called K.f, K.l, K.q, K.r and
       K.s, for Free, Linear, Quadratic, Rotated quadratic and Semi-definite.
       In addition, there are fields K.xcomplex, K.scomplex and K.ycomplex
       for complex-variables.
 
     (1) K.f is the number of FREE, i.e. UNRESTRICTED primal components.
       The dual components are restricted to be zero. E.g. if
       K.f = 2 then x(1:2) is unrestricted, and z(1:2)=0.
       These are ALWAYS the first components in x.
 
     (2) K.l is the number of NONNEGATIVE components. E.g. if K.f=2, K.l=8
       then x(3:10) >=0.
 
     (3) K.q lists the dimensions of LORENTZ (quadratic, second-order cone)
       constraints. E.g. if K.l=10 and K.q = [3 7] then
           x(11) >= norm(x(12:13)),
           x(14) >= norm(x(15:20)).
       These components ALWAYS immediately follow the K.l nonnegative ones.
       If the entries in A and/or c are COMPLEX, then the x-components in
       "norm(x(#,#))" take complex-values, whenever that is beneficial.
        Use K.ycomplex to impose constraints on the imaginary part of A*x.
 
     (4) K.r lists the dimensions of Rotated LORENTZ
       constraints. E.g. if K.l=10, K.q = [3 7] and K.r = [4 6], then
           2*x(21)x(22) >= norm(x(23:24))^2,
           2*x(25)x(26) >= norm(x(27:30))^2.
       These components ALWAYS immediately follow the K.q ones.
       Just as for the K.q-variables, the variables in "norm(x(#,#))" are
       allowed to be complex, if you provide complex data. Use K.ycomplex
       to impose constraints on the imaginary part of A*x.
 
     (5) K.s lists the dimensions of POSITIVE SEMI-DEFINITE (PSD) constraints
       E.g. if K.l=10, K.q = [3 7] and K.s = [4 3], then
           mat( x(21:36),4 ) is PSD,
           mat( x(37:45),3 ) is PSD.
       These components are ALWAYS the last entries in x.
 
     (a) K.xcomplex lists the components in f,l,q,r blocks that are allowed
      to have nonzero imaginary part in the primal. For f,l blocks, these
     (b) K.scomplex lists the PSD blocks that are Hermitian rather than
       real symmetric.
     (c) Use K.ycomplex to impose constraints on the imaginary part of A*x.
 
     The dual multipliers y have analogous meaning as in the "x>=0" case,
     except that instead of "c-A'*y>=0" resp. "-A'*y>=0", one should read that
     c-A'*y resp. -A'*y are in the cone that is described by K.l, K.q and K.s.
     In the above example, if z = c-A'*y and mat(z(21:36),4) is not symmetric/
     Hermitian, then positive semi-definiteness reflects the symmetric/
     Hermitian parts, i.e. Z + Z' is PSD.
 
     If the model contains COMPLEX data, then you may provide a list
     K.ycomplex, with the following meaning:
       y(i) is complex if ismember(i,K.ycomplex)
       y(i) is real otherwise
     The equality constraints in the primal are then as follows:
           A(i,:)*x = b(i)      if imag(b(i)) ~= 0 or ismember(i,K.ycomplex)
           real(A(i,:)*x) = b(i)  otherwise.
     Thus, equality constraints on both the real and imaginary part
     of A(i,:)*x should be listed in the field K.ycomplex.
 
     You may use EIGK(x,K) and EIGK(c-A'*y,K) to check that x and c-A'*y
     are in the cone K.
 
  >  [X,Y,INFO] = SEDUMI(A,b,c,K,pars) allows you to override the default
       parameter settings, using fields in the structure `pars'.
 
     (1) pars.fid     By default, fid=1. If fid=0, then SeDuMi runs quietly,
       i.e. no screen output. In general, output is written to a device or
       file whose handle is fid. Use fopen to assign a fid to a file.
 
     (2) pars.alg     By default, alg=2. If alg=0, then a first-order wide
       region algorithm is used, not recommended. If alg=1, then SeDuMi uses
       the centering-predictor-corrector algorithm with v-linearization.
       If alg=2, then xz-linearization is used in the corrector, similar
       to Mehrotra's algorithm. The wide-region centering-predictor-corrector
       algorithm was proposed in Chapter 7 of
         J.F. Sturm, Primal-Dual Interior Point Approach to Semidefinite Pro-
         gramming, TIR 156, Thesis Publishers Amsterdam, 1997.
 
     (3) pars.theta, pars.beta   By default, theta=0.25 and beta=0.5. These
       are the wide region and neighborhood parameters. Valid choices are
       0 < theta <= 1 and 0 < beta < 1. Setting theta=1 restricts the iterates
       to follow the central path in an N_2(beta)-neighborhood.
 
     (4) pars.stepdif, pars.w. By default, stepdif = 2 and w=[1 1].
        This implements an adaptive heuristic to control ste-differentiation.
        You can enable primal/dual step length differentiation by setting stepdif=1 or 0.
       If so, it weights the rel. primal, dual and gap residuals as
       w(1):w(2):1 in order to find the optimal step differentiation.
 
     (5) pars.eps     The desired accuracy. Setting pars.eps=0 lets SeDuMi run
       as long as it can make progress. By default eps=1e-8.
 
     (6) pars.bigeps  In case the desired accuracy pars.eps cannot be achieved,
      the solution is tagged as info.numerr=1 if it is accurate to pars.bigeps,
      otherwise it yields info.numerr=2.
 
     (7) pars.maxiter Maximum number of iterations, before termination.
 
     (8) pars.stopat  SeDuMi enters debugging mode at the iterations specified in this vector.
 
     (9) pars.cg      Various parameters for controling the Preconditioned conjugate
      gradient method (CG), which is only used if results from Cholesky are inaccurate.
     (a) cg.maxiter   Maximum number of CG-iterates (per solve). Theoretically needed
           is |add|+2*|skip|, the number of added and skipped pivots in Cholesky.
           (Default 49.)
     (b) cg.restol    Terminates if residual is a "cg.restol" fraction of duality gap.
           Should be smaller than 1 in order to make progress (default 5E-3).
     (c) cg.refine    Number of refinement loops that are allowed. The maximum number
           of actual CG-steps will thus be 1+(1+cg.refine)*cg.maxiter. (default 1)
     (d) cg.stagtol  Terminates if relative function progress less than stagtol (5E-14).
     (e) cg.qprec    Stores cg-iterates in quadruple precision if qprec=1 (default 0).
 
     (10) pars.chol   Various parameters for controling the Cholesky solve.
      Subfields of the structure pars.chol are:
     (a) chol.canceltol: Rel. tolerance for detecting cancelation during Cholesky (1E-12)
     (b) chol.maxu:   Adds to diagonal if max(abs(L(:,j))) > chol.maxu otherwise (5E5).
     (c) chol.abstol: Skips pivots falling below abstol (1e-20).
     (d) chol.maxuden: pivots in dense-column factorization so that these factors
       satisfy max(abs(Lk)) <= maxuden (default 5E2).
 
     (11) pars.vplot  If this field is 1, then SeDuMi produces a fancy
       v-plot, for research purposes. Default: vplot = 0.
 
     (12) pars.errors  If this field is 1 then SeDuMi outputs some error
     measures as defined in the Seventh DIMACS Challenge. For more details
     see the User Guide.
 
  Bug reports can be submitted at http://sedumi.mcmaster.ca.
 
  See also mat, vec, cellK, eyeK, eigK





posted by maetel

http://en.wikipedia.org/wiki/Linear_programming


1.1 The Linear Programming Problem



general  linear programming problem
- objective function
- constraints

standard form
canonical form

Every linear programming problem that has unconstrained variables can be solved by solving a corresponding linear programming problem in which all the variables are constrained to be nonnegative.

Every linear programming problem can be formulated as a corresponding standard linear programming problem or as a corresponding canonical linear programming problem.   (53p)


Minmization Problem as a Maximization Problem
: To minimize the objective function we could maximize its negative instead and then change the sign of the answer.

Reversing an Inequality

Changing an Equality to Inequality

Unconstrained Variables

Scaling
to make all coefficients in a linear programming problem approximately the same size



1.2 Matrix Notation


feasible solution
: a vector satisfying the constraints of a linear programming problem

optimal solution
: a feasible solution maximizing or minimizing the objective function of a linear programming



1.4




ref.
http://en.wikipedia.org/wiki/Convex_polyhedron
http://mathworld.wolfram.com/ConvexPolyhedron.html
http://en.wikipedia.org/wiki/Extreme_point
http://en.wikipedia.org/wiki/Krein%E2%80%93Milman_theorem
http://library.wolfram.com/infocenter/MathSource/440/
http://www.ifor.math.ethz.ch/~fukuda/fukuda.html



slack & surplus variable in the siimplex algorithm
http://www-fp.mcs.anl.gov/otc/Guide/CaseStudies/simplex/standard.html

posted by maetel
By Bernard Kolman, Robert Edward Beck
Contributor Robert Edward Beck

Edition: 2, illustrated
Published by Academic Press, 1995
ISBN 012417910X, 9780124179103
449 pages



Prologue
: Introduction to Operations Research




Definitions

Operations research (OR) is a scientific method for providing a quantitative basis for decision making (that can be used in almost any field of endeavor).

http://en.wikipedia.org/wiki/Operations_research

The techniques of OR give a logical and systematic way of formulating a problem so that the tools of mathematics can be applied to find a solution (when the established way of doing things can be improved.

A central problem in OR is the optimal allocation of scarce resources (including raw materials, labor, capital, energy, and processing time).

eg. T. C. Koopmans and L. V. Kantorovich, the 1975 Nobel Prize awardees



Development

WWII, Great Britain
the United States Air Force
1947, George B. Dantzig, the simplex algorithm (for solving linear programming problems)
programmable digital computer (for solving large-scale linear programming problems)
1952, the National Bureau of Standards SEAC machine



Phases

> Steps
1. Problem definition and formulation
to clarify objectives in undertaking the study
to identify the decision alternatives
to consider the limitations, restrictions, and requirement of the various alternatives
2. Model construction
to develop the appropriate mathematical description of the problem
: limitations, restrictions, and requirements -> constraints
: the goal of the study -> quantified as to be maximized or minimized
: decision alternatives -> variables
3. Solution of the model
: the solution to the model need not be the solution to the original problem
4. Sensitivity analysis
to determine the sensitivity of the solution to changes in the inpuit data
: how the solution will vary with the variation in the input data
5. Model evalution
to determine whether the answers produced by the model are realistic, acceptable and capable of implementation
6. Implementation of the study



The Structure of Mathematical Models

There are two  types of mathematical models: deterministic and probabilistic.

A deterministic model will always yield the same set of output values for a given set of input values, whereas a probabilistic model will typically yield many different sets of output values according to some probability distribution.


Decision variables or unknowns
to describe an optimal allocation of the scarce resources represented by the model
Parameters
: inputs, not adjustable but exactly or approximately known
Constraints
: conditions that limit the values that the decision variables can assume
Objective function
to measures the effectiveness of the system as a function of the decision variables
: The decisin variables are to be determined so that the objective function will be optimized.



Mathematical Techniques

mathematical programming

Mathematical models call for finding values of the decision variables that maximize or minimize the objective function subject to a system of inequality and equality constraints.

Linear programming

Integer programming

Stochastic programming

Nonlinear programming

network flow analysis

standard techniques vs. heuristic techniques (improvised for the particular problem and without firm mathematical basis)


Journals

Computer and Information Systems Abstracts Journal
Computer Journal
Decision Sciences
European Journal of Operational Research
IEEE Transactions on Automatic Control
Interfaces
International Abstracts in Operations Research
Journal of Computer and System Sciences
Journal of Research of the National Bureau of Standards
Journal of the ACM
Journal of the Canadian Operational Research Society
Management Science (published by The Institute for Management Science - TIMS)
Mathematical Programming
Mathematics  in Operations Research
Naval Research Logisics (published by the Office of Naval Research - ONR)
Operational Research Quarterly
Operations Research (published by the Operations Research Society of America - ORSA)
Operations Research Letters
OR/MS Today
ORSA Journal on Computing
SIAM Journal on Computing
Transportation Science
posted by maetel
4.1 Optimization Problems



eg.4.1


ezplot('x*log(x)')



* MATLAB에서 그래프 그리기

ezplot('x*log(x)')
로 해서 나왔으나, 일반적으로는 다음과 같은 식으로 한다.

x=0:0.1:10;
y=x.*log2(x);  ## 여기서, x 다음의 .은 element끼리의 scalar 연산을 지시함
plot(x,y)






http://en.wikipedia.org/wiki/Slack_variable


4.2 Convex Optimization




4.3 Linear Optimization Problems

http://en.wikipedia.org/wiki/Linear_program

 

 

4.4 Quadratic Optimization Problems

http://en.wikipedia.org/wiki/Quadratic_program


Markowitz portfolio optimization

 

4.4.2 Second-order cone programming

http://en.wikipedia.org/wiki/SOCP



http://en.wikipedia.org/wiki/Minimal_surface







posted by maetel
기고가 Lieven Vandenberghe
출판사: Cambridge University Press, 2004
ISBN 0521833787, 9780521833783
posted by maetel