it-swarm-korea.com

2D 점을 3D로 반전 투영하는 방법은 무엇입니까?

화면 공간에 4 개의 2D 점이 있으며 3D 공간으로 다시 투영해야합니다. 4 점 각각이 3D 회전 고정 사각형의 모서리라는 것을 알고 사각형의 크기를 알고 있습니다. 이로부터 3D 좌표를 얻으려면 어떻게해야합니까?

특정 API를 사용하지 않고 기존 프로젝션 매트릭스가 없습니다. 나는 이것을하기 위해 기본 수학을 찾고 있습니다. 물론 다른 참조가없는 단일 2D 점을 3D로 변환하기에는 데이터가 충분하지 않지만 4 점이 있으면 같은 평면에서 서로 직각을 이루고 있다는 것을 알고 있습니다. 그들 사이의 거리를 알면 거기서부터 알아낼 수 있어야합니다. 불행히도 나는 어떻게 해결할 수는 없습니다.

이것은 사진 측량의 우산에 해당 할 수 있지만 Google은 검색하여 유용한 정보를 얻지 못했습니다.

56
Joshua Carmody

좋아, 나는 답을 찾고 여기에 왔고 간단하고 직접적인 것을 찾지 못해서 멍청하지만 효과적이면서도 상대적으로 간단한 것을 몬테카를로 최적화했습니다.

간단히 말해서 알고리즘은 다음과 같습니다. 투영 된 행렬이 알려진 3D 좌표를 알려진 2D 좌표에 투영 할 때까지 투영 행렬을 무작위로 교란시킵니다.

다음은 Thomas the Tank Engine의 스틸 사진입니다.

Thomas the Tank Engine

GIMP를 사용하여지면에서 정사각형이라고 생각하는 것의 2D 좌표를 찾는다고 가정 해 봅시다 (실제로 정사각형인지 여부는 깊이에 대한 판단에 달려 있습니다).

With an outline of the square

2D 이미지에 (318, 247), (326, 312), (418, 241)(452, 303)의 4 가지 포인트가 있습니다.

일반적으로이 포인트는 3D 포인트 (0, 0, 0), (0, 0, 1), (1, 0, 0)(1, 0, 1)와 일치해야합니다. 즉, y = 0 평면의 단위 제곱입니다.

이 3D 좌표를 2D로 투영하는 것은 4D 벡터 [x, y, z, 1]에 4x4 프로젝션 매트릭스를 곱한 다음 x와 y 성분을 z로 나누어 원근 교정을 실제로 수행함으로써 수행됩니다. gluProject()도 현재 뷰포트를 고려하고 별도의 모델 뷰 매트릭스를 고려한다는 점을 제외하고 gluProject () 의 기능은 다소 다릅니다. 동일성 행렬). 실제로 OpenGL에서 작동하는 솔루션을 원하기 때문에 gluProject() 설명서를 보는 것이 매우 편리하지만 설명서에서 수식이 z로 나뉘어 있지 않다는 점에주의하십시오.

이 알고리즘은 일부 투영 행렬로 시작하여 원하는 투영이 나올 때까지 임의로 교란시키는 것입니다. 우리가 할 일은 4 개의 3D 점 각각을 투영하고 원하는 2D 점에 얼마나 가까운 지 확인하는 것입니다. 임의의 섭동으로 인해 투영 된 2D 점이 위에 표시된 점에 더 가까워지면 해당 행렬을 초기 (또는 이전) 추측보다 개선 된 것으로 유지합니다.

우리의 요점을 정의합시다 :

# Known 2D coordinates of our rectangle
i0 = Point2(318, 247)
i1 = Point2(326, 312)
i2 = Point2(418, 241)
i3 = Point2(452, 303)

# 3D coordinates corresponding to i0, i1, i2, i3
r0 = Point3(0, 0, 0)
r1 = Point3(0, 0, 1)
r2 = Point3(1, 0, 0)
r3 = Point3(1, 0, 1)

우리는 몇 가지 매트릭스로 시작해야합니다. 항등 매트릭스는 자연스러운 선택입니다

mat = [
    [1, 0, 0, 0],
    [0, 1, 0, 0],
    [0, 0, 1, 0],
    [0, 0, 0, 1],
]

투영 (기본적으로 행렬 곱셈)을 구현해야합니다.

def project(p, mat):
    x = mat[0][0] * p.x + mat[0][1] * p.y + mat[0][2] * p.z + mat[0][3] * 1
    y = mat[1][0] * p.x + mat[1][1] * p.y + mat[1][2] * p.z + mat[1][3] * 1
    w = mat[3][0] * p.x + mat[3][1] * p.y + mat[3][2] * p.z + mat[3][3] * 1
    return Point(720 * (x / w + 1) / 2., 576 - 576 * (y / w + 1) / 2.)

이것은 기본적으로 gluProject()이하는 것, 720과 576은 각각 이미지의 너비와 높이 (예 : 뷰포트)이며 576에서 빼서 y 좌표를 맨 위에서 세 었다는 사실을 계산합니다. OpenGL은 일반적으로 맨 아래부터 계산합니다. z를 계산하지 않는다는 것을 알 수 있습니다. 왜냐하면 여기서 실제로 필요하지 않기 때문입니다 (OpenGL이 깊이 버퍼에 사용하는 범위 내에 있도록하는 것이 편리 할 수는 있지만).

이제 올바른 솔루션에 얼마나 가까운 지 평가할 수있는 기능이 필요합니다. 이 함수가 반환하는 값은 한 행렬이 다른 행렬보다 나은지 확인하는 데 사용됩니다. 나는 제곱 거리의 합계로 가기로 결정했습니다.

# The squared distance between two points a and b
def norm2(a, b):
    dx = b.x - a.x
    dy = b.y - a.y
    return dx * dx + dy * dy

def evaluate(mat): 
    c0 = project(r0, mat)
    c1 = project(r1, mat)
    c2 = project(r2, mat)
    c3 = project(r3, mat)
    return norm2(i0, c0) + norm2(i1, c1) + norm2(i2, c2) + norm2(i3, c3)

행렬을 교란시키기 위해, 어떤 범위 내에서 임의의 양만큼 교란 할 요소를 선택하기 만하면됩니다 :

def perturb(amount):
    from copy import deepcopy
    from random import randrange, uniform
    mat2 = deepcopy(mat)
    mat2[randrange(4)][randrange(4)] += uniform(-amount, amount)

(우리는 project() 함수가 z를 계산하지 않기 때문에 실제로 mat[2]를 전혀 사용하지 않으며 모든 y 좌표가 0이므로 mat[*][1] 값도 무의미합니다.이 사실을 사용하여 그 값을 섭동하려고 시도하면 속도가 약간 떨어지지 만 운동으로 남습니다.)

편의상 지금까지 찾은 최고의 행렬에 대해 perturb()을 반복해서 호출하여 근사를 수행하는 함수를 추가해 보겠습니다.

def approximate(mat, amount, n=100000):
    est = evaluate(mat)

    for i in xrange(n):
        mat2 = perturb(mat, amount)
        est2 = evaluate(mat2)
        if est2 < est:
            mat = mat2
            est = est2

    return mat, est

이제 남은 일은 그것을 실행하는 것입니다.

for i in xrange(100):
    mat = approximate(mat, 1)
    mat = approximate(mat, .1)

나는 이것이 이미 꽤 정확한 답변을 제공한다는 것을 알았습니다. 잠시 동안 실행 한 후 찾은 매트릭스는 다음과 같습니다.

[
    [1.0836000765696232,  0,  0.16272110011060575, -0.44811064935115597],
    [0.09339193527789781, 1, -0.7990570384334473,   0.539087345090207  ],
    [0,                   0,  1,                    0                  ],
    [0.06700844759602216, 0, -0.8333379578853196,   3.875290562060915  ],
]

2.6e-5 오류가 발생했습니다. (우리가 계산에 사용되지 않은 요소가 실제로 초기 행렬에서 변경되지 않았다는 점에 유의하십시오. 이러한 항목을 변경해도 평가 결과가 변경되지 않으므로 변경 사항이 적용되지 않기 때문입니다.)

glLoadMatrix()을 사용하여 행렬을 OpenGL에 전달할 수 있습니다 (그러나 먼저 변환하고 모델 행렬을 항등 행렬과 함께로드하는 것을 잊지 마십시오).

def transpose(m):
    return [
        [m[0][0], m[1][0], m[2][0], m[3][0]],
        [m[0][1], m[1][1], m[2][1], m[3][1]],
        [m[0][2], m[1][2], m[2][2], m[3][2]],
        [m[0][3], m[1][3], m[2][3], m[3][3]],
    ]

glLoadMatrixf(transpose(mat))

이제 예를 들어 z 축을 따라 이동하여 트랙을 따라 다른 위치를 얻을 수 있습니다.

glTranslate(0, 0, frame)
frame = frame + 1

glBegin(GL_QUADS)
glVertex3f(0, 0, 0)
glVertex3f(0, 0, 1)
glVertex3f(1, 0, 1)
glVertex3f(1, 0, 0)
glEnd()

With 3D translation

확실히 이것은 수학적 관점에서 매우 우아하지 않습니다. 숫자를 꽂고 직접적이고 정확한 답변을 얻을 수있는 닫힌 형식 방정식을 얻지 못합니다. 그러나 방정식의 복잡성에 대해 걱정할 필요없이 제약 조건을 추가 할 수 있습니다. 예를 들어 높이를 통합하려면 집의 모서리를 사용하고 평가 기능에서지면에서 지붕까지의 거리가 적절해야한다고 말하고 알고리즘을 다시 실행하십시오. 예, 그것은 일종의 무차별적인 힘이지만 효과가 있으며 잘 작동합니다.

Choo choo!

67
Vegard

이것은 마커 기반 증강 현실에 대한 고전적인 문제입니다.

사각형 마커 (2D 바코드)가 있고 마커의 네 모서리를 찾은 후 포즈 (카메라와 관련된 변환 및 회전)를 찾으려고합니다. 개요-사진

나는이 분야에 대한 최신 공헌을 알지 못하지만 적어도 한 시점 (2009 년)까지 RPP는 위에서 언급 한 POSIT보다 성능이 뛰어나야한다고 생각했다. 소스를 제공하십시오.

(PS-약간 오래된 주제라는 것을 알고 있지만 어쨌든 게시물이 누군가에게 도움이 될 수 있습니다)

7
dim_tz

D. DeMenthon 객체의 모델을 알 때 2D 이미지의 특징점에서 객체의 pose (공간에서의 위치 및 방향)을 계산하는 알고리즘을 고안했습니다. 이것은 당신의 정확한 문제입니다 :

단일 이미지에서 객체의 포즈를 찾는 방법을 설명합니다. 우리는 이미지에서 물체의 4 개 이상의 비평면 특징점을 감지하고 일치시킬 수 있고 물체의 상대적인 지오메트리를 알고 있다고 가정합니다.

이 알고리즘은 Posit로 알려져 있으며, 고전 기사 인 "25 줄의 코드에서 모델 기반 객체 포즈"( 웹 사이트 , 섹션 4에서 사용 가능)에 설명되어 있습니다.

기사로 직접 연결되는 링크 : http://www.cfar.umd.edu/~daniel/daniel_papersfordownload/Pose25Lines.pdf OpenCV 구현 : http://opencv.willowgarage.com/ 위키/위치

아이디어는 정확한 자세로 수렴 될 때까지 scaled orthographic projection으로 원근 투영법을 반복적으로 근사화하는 것입니다.

5
Julien-L

2 차원 공간에서 2 개의 유효한 직사각형을 만들 수 있습니다. 원래의 매트릭스 투영을 알지 못하면 어느 것이 올바른지 알 수 없습니다. "상자"문제와 동일합니다. 두 개의 정사각형 (하나는 다른 내부에 있음)과 4 개의 내부 정점이 4 개의 외부 정점에 연결되어 있습니다. 위에서 아래로 또는 아래에서 위로 상자를보고 있습니까?

즉, 당신은 매트릭스 변환 T를 찾고 있습니다 ...

{{x1, y1, z1}, {x2, y2, z2}, {x3, y3, z3}, {x4, y4, z4}} x T = {{x1, y1}, {x2, y2}, { x3, y3}, {x4, y4}}

(4 x 3) x T = (4 x 2)

따라서 T는 (3 x 2) 행렬이어야합니다. 6 개의 미지수를 얻었습니다.

이제 T에 제약 시스템을 구축하고 Simplex로 해결하십시오. 구속 조건을 작성하려면 처음 두 점을 통과하는 선이 두 번째 두 점을 통과하는 선과 평행해야합니다. 점 1과 3을 통과하는 선은 점 2와 4를 통과하는 선과 평행해야한다는 것을 알고 있습니다. 1과 2를 통과하는 선은 점 2와 3을 통과하는 선과 직교해야한다는 것을 알고 있습니다. 1과 2의 줄 길이는 3과 4의 줄 길이와 같아야합니다. 1과 3의 줄 길이는 2와 4의 줄 길이와 같아야합니다.

이것을 더 쉽게하기 위해 직사각형에 대해 알고 있으므로 모든면의 길이를 알 수 있습니다.

그것은이 문제를 해결하기 위해 많은 제약을 줄 것입니다.

물론, 다시 돌아가려면 T-inverse를 찾을 수 있습니다.

@Rob : 예, 무한한 개수의 투영이 있지만 점이 사각형의 요구 사항을 충족해야하는 무한한 개수의 프로젝트는 없습니다.

@ nlucaroni : 예, 투영에 4 점이있는 경우에만 해결할 수 있습니다. 사각형이 2 포인트로만 투영되는 경우 (즉, 사각형의 평면이 투영 표면과 직교하는 경우)이를 해결할 수 없습니다.

흠 ... 집에 가서이 작은 보석을 써야 해요. 재미있을 것 같습니다.

업데이트 :

  1. 포인트 중 하나를 수정하지 않으면 무한한 투영이 있습니다. 원래 사각형의 점을 고정하면 두 가지 가능한 원래 사각형이 있습니다.
4
Jarrett Meyer

내 OpenGL 엔진의 경우 다음 코드 조각은 마우스/스크린 좌표를 3D 월드 좌표로 변환합니다. 진행 상황에 대한 실제 설명은 주석을 읽으십시오.

/* 기능 : YCamera :: CalculateWorldCoordinates 
 인수 : x 마우스 x 좌표 
 y 마우스 y 좌표 
 vec 좌표를 저장할 위치 
 반환 : 해당 없음 
 설명 : 마우스 좌표를 세계 좌표로 변환 
 */

void YCamera :: CalculateWorldCoordinates(float x, float y, YVector3 *vec) { // START GLint viewport[4]; GLdouble mvmatrix[16], projmatrix[16];

GLint real_y;
GLdouble mx, my, mz;

glGetIntegerv(GL_VIEWPORT, viewport);
glGetDoublev(GL_MODELVIEW_MATRIX, mvmatrix);
glGetDoublev(GL_PROJECTION_MATRIX, projmatrix);

real_y = viewport[3] - (GLint) y - 1;   // viewport[3] is height of window in pixels
gluUnProject((GLdouble) x, (GLdouble) real_y, 1.0, mvmatrix, projmatrix, viewport, &mx, &my, &mz);

/*  'mouse' is the point where mouse projection reaches FAR_PLANE.
    World coordinates is intersection of line(camera->mouse) with plane(z=0) (see LaMothe 306)

    Equation of line in 3D:
        (x-x0)/a = (y-y0)/b = (z-z0)/c      

    Intersection of line with plane:
        z = 0
        x-x0 = a(z-z0)/c  <=> x = x0+a(0-z0)/c  <=> x = x0 -a*z0/c
        y = y0 - b*z0/c

*/
double lx = fPosition.x - mx;
double ly = fPosition.y - my;
double lz = fPosition.z - mz;
double sum = lx*lx + ly*ly + lz*lz;
double normal = sqrt(sum);
double z0_c = fPosition.z / (lz/normal);

vec->x = (float) (fPosition.x - (lx/normal)*z0_c);
vec->y = (float) (fPosition.y - (ly/normal)*z0_c);
vec->z = 0.0f;
</ code>

}

4
user14208

Rons 접근 방식에 대한 후속 조치 : 사각형을 어떻게 회전했는지 알고 있으면 z 값을 찾을 수 있습니다.

요령은 투영을 한 투영 행렬을 찾는 것입니다. 다행히도 이것은 가능하고 심지어 저렴합니다. 관련 수학은 Paul Heckbert의 "이미지 워핑에 대한 투영 매핑"논문에서 찾을 수 있습니다.

http://pages.cs.wisc.edu/~dyer/cs766/readings/heckbert-proj.pdf

이렇게하면 투영 중에 손실 된 각 정점의 균일 한 부분을 복구 할 수 있습니다.

이제 Ron이 설명했듯이 여전히 포인트 대신 네 개의 선이 남아 있습니다. 그러나 원래 사각형의 크기를 알고 있기 때문에 아무것도 손실되지 않습니다. 이제 Ron의 방법과 2D 접근법의 데이터를 선형 방정식 솔버에 연결하고 z에 대해 풀 수 있습니다. 그런 식으로 각 정점의 정확한 z 값을 얻습니다.

참고 : 이것은 다음과 같은 이유로 작동합니다.

  1. 원래 모양은 사각형이었습니다
  2. 3D 공간에서 직사각형의 정확한 크기를 알고 있습니다.

정말 특별한 경우입니다.

도움이 되었으면 좋겠다, 닐스

2
Nils Pipenbrinck

점이 실제로 사각형의 일부라고 가정하면 일반적인 아이디어를 제공합니다.

최대 거리로 두 점을 찾으십시오. 대각선을 정의하는 것 같습니다 (예외 : 직사각형이 YZ 평면에 거의 평행하고 학생을 위해 남겨진 특수한 경우). 그것들을 A, C라고 부르십시오. BAD, BCD 각도를 계산하십시오. 이들은 직각과 비교하여 3D 공간에서 방향을 제공합니다. z 거리에 대해 알아 보려면 투영 된면을 알려진면과 상관시켜야하며, 3D 투영 방법 (1/z)을 기준으로 거리를 알기 위해 올바른 길을 가고 있습니다.

2
tzot

아무도 대답하지 않으면 집에 도착하면 선형 대수 책을 꺼낼 것입니다. 그러나 @ D G, 모든 행렬이 돌이킬 수있는 것은 아닙니다. 단일 행렬은 되돌릴 수 없음 (결정자 = 0 일 때). 투영 행렬 must 고유 값이 0과 1이므로 제곱이되기 때문에 실제로 항상 발생합니다.

쉬운 예는 [결정자 = 0이므로 [[0 1] [0 1]]입니다. 이는 선 x = y에 대한 투영입니다!

1
nlucaroni

2D 표면에 투영 한 투영에는 동일한 2D 모양으로 투영되는 3D 직사각형이 무한히 많습니다.

이런 식으로 생각하십시오. 3D 직사각형을 구성하는 4 개의 3D 점이 있습니다. (x0, y0, z0), (x1, y1, z1), (x2, y2, z2) 및 (x3, y3, z3)이라고합니다. 이러한 점을 x-y 평면에 투영하면 (x0, y0), (x1, y1), (x2, y2), (x3, y3)과 같은 z 좌표가 삭제됩니다.

이제 3D 공간으로 다시 투영하려면 z0, .., z3이 무엇인지 리버스 엔지니어링해야합니다. 그러나 a) 점 사이의 동일한 x-y 거리를 유지하고 b) 사각형이 작동하는 모양을 유지하는 모든 z 좌표 세트. 따라서이 (무한한) 세트의 멤버는 {(z0 + i, z1 + i, z2 + i, z3 + i) | i <-R}.

@Jarrett 편집 :이 문제를 해결하고 3D 공간에서 직사각형으로 끝났다고 상상해보십시오. 이제 사각형을 z 축 위아래로 슬라이드한다고 상상해보십시오. 변환 된 사각형의 무한한 양은 모두 동일한 x-y 투영을 갖습니다. "올바른"것을 찾았다는 것을 어떻게 알 수 있습니까?

편집 # 2 : 좋아, 이것은 내가이 질문에 대한 의견에서 나온 것입니다-이것에 대한 추론에 대한보다 직관적 인 접근법.

책상 위에 종이 한 장을 들고 있다고 상상해보십시오. 종이의 각 구석에는 책상을 향하는 무중력 레이저 포인터가 부착되어 있다고 가정하십시오. 종이는 3D 객체이고 책상의 레이저 포인터 도트는 2D 투영입니다.

이제 레이저 포인터 도트를 just를 보면 용지가 얼마나 높은지 알 수 있습니까?

당신은 할 수 없습니다. 용지를 똑바로 위아래로 움직입니다. 레이저 포인터는 용지 높이와 상관없이 책상 위의 동일한 지점에서 계속 빛납니다.

역 투영에서 z 좌표를 찾는 것은 책상 위의 레이저 포인터 점을 기준으로 용지 높이를 찾는 것과 같습니다.

1
Rob Dickerson
1
Ray Tayek

모양이 평면의 사각형이라는 것을 알고 있으면 문제를 더욱 제한 할 수 있습니다. 확실하게 "어떤"평면을 파악할 수 없으므로 z = 0이고 모서리 중 하나가 x = y = 0 인 평면에 있고 모서리가 x/y 축에 평행하도록 선택할 수 있습니다.

따라서 3d의 점은 {0,0,0}, {w, 0,0}, {w, h, 0} 및 {0, h, 0}입니다. 나는 절대 크기가 발견되지 않을 것이라고 확신합니다. 따라서 비율 w/h만이 해제되므로 이것은 알려지지 않은 것입니다.

이 평면과 관련하여 카메라는 공간에서 cx, cy, cz 지점에 있어야하며 nx, ny, nz 방향을 가리켜 야합니다 (길이가 1 인 벡터는 중복 됨). focal_length/image_width w의 요인. 이 숫자는 3x3 투영 행렬로 바뀝니다.

이는 총 7 개의 미지수 (w/h, cx, cy, cz, nx, ny 및 w)를 제공합니다.

4 개의 x + y 쌍으로 총 8 개의 알려진 항목이 있습니다.

따라서이 문제를 해결할 수 있습니다.

다음 단계는 Matlab 또는 Mathmatica를 사용하는 것입니다.

1
spitzak

네, 몬테 카를로가 효과가 있지만이 문제에 대한 더 나은 해결책을 찾았습니다. 이 코드는 완벽하게 작동하며 OpenCV를 사용합니다.

Cv2.CalibrateCamera(new List<List<Point3f>>() { points3d }, new List<List<Point2f>>() { points2d }, new Size(height, width), cameraMatrix, distCoefs, out rvecs, out tvecs, CalibrationFlags.ZeroTangentDist | CalibrationFlags.FixK1 | CalibrationFlags.FixK2 | CalibrationFlags.FixK3);

이 함수는 알려진 3d 및 2d 포인트, 화면 크기를 가져와 회전 (rvecs [0]), 변환 (tvecs [0]) 및 카메라 내장 값의 행렬을 반환합니다. 필요한 모든 것입니다.

1
Inflight

3D에서 2D로 투사하면 정보가 손실됩니다.

단일 점의 간단한 경우에 역투 영은 3D 공간을 통해 무한한 광선을 제공합니다.

입체 재구성은 일반적으로 2 개의 2D 이미지로 시작하여 3D로 다시 투영됩니다. 그런 다음 생성 된 두 3D 광선의 교차점을 찾으십시오.

투영은 다른 형태를 취할 수 있습니다. 직교 또는 원근법. 나는 당신이 직교 투영을 가정하고 있다고 추측하고 있습니까?

원래 매트릭스가 있다고 가정하면 3D 공간에 4 개의 광선이 생깁니다. 그런 다음 3D 직사각형 치수로 문제를 제한하고 해결할 수 있습니다.

2D 투영 평면에 평행 한 두 축을 중심으로 한 회전이 방향이 모호하므로 솔루션은 고유하지 않습니다. 다시 말해서 2d 이미지가 z 축에 수직 인 경우 3d 사각형을 x 축을 기준으로 시계 방향 또는 시계 반대 방향으로 회전하면 동일한 이미지가 생성됩니다. y 축도 마찬가지입니다.

사각형 평면이 z 축과 평행 한 경우 더 많은 솔루션이 있습니다.

원래 영사 행렬이 없으므로 모든 영사에 존재하는 임의의 스케일링 계수에 의해 모호성이 추가됩니다. 투영의 배율과 z 축 방향의 3d 평행 이동을 구분할 수 없습니다. 2D 투영 평면이 아닌 서로 관련 될 때 3D 공간에서 4 개의 점의 상대 위치에만 관심이있는 경우에는 문제가되지 않습니다.

투시 투영에서 일이 더 어려워집니다 ...

1
morechilli

훌륭한 답변을 주신 @Vegard에게 감사드립니다. 코드를 약간 정리했습니다.

import pandas as pd
import numpy as np

class Point2:
    def __init__(self,x,y):
        self.x = x
        self.y = y

class Point3:
    def __init__(self,x,y,z):
        self.x = x
        self.y = y
        self.z = z

# Known 2D coordinates of our rectangle
i0 = Point2(318, 247)
i1 = Point2(326, 312)
i2 = Point2(418, 241)
i3 = Point2(452, 303)

# 3D coordinates corresponding to i0, i1, i2, i3
r0 = Point3(0, 0, 0)
r1 = Point3(0, 0, 1)
r2 = Point3(1, 0, 0)
r3 = Point3(1, 0, 1)

mat = [
    [1, 0, 0, 0],
    [0, 1, 0, 0],
    [0, 0, 1, 0],
    [0, 0, 0, 1],
]

def project(p, mat):
    #print mat
    x = mat[0][0] * p.x + mat[0][1] * p.y + mat[0][2] * p.z + mat[0][3] * 1
    y = mat[1][0] * p.x + mat[1][1] * p.y + mat[1][2] * p.z + mat[1][3] * 1
    w = mat[3][0] * p.x + mat[3][1] * p.y + mat[3][2] * p.z + mat[3][3] * 1
    return Point2(720 * (x / w + 1) / 2., 576 - 576 * (y / w + 1) / 2.)

# The squared distance between two points a and b
def norm2(a, b):
    dx = b.x - a.x
    dy = b.y - a.y
    return dx * dx + dy * dy

def evaluate(mat): 
    c0 = project(r0, mat)
    c1 = project(r1, mat)
    c2 = project(r2, mat)
    c3 = project(r3, mat)
    return norm2(i0, c0) + norm2(i1, c1) + norm2(i2, c2) + norm2(i3, c3)    

def perturb(mat, amount):
    from copy import deepcopy
    from random import randrange, uniform
    mat2 = deepcopy(mat)
    mat2[randrange(4)][randrange(4)] += uniform(-amount, amount)
    return mat2

def approximate(mat, amount, n=1000):
    est = evaluate(mat)
    for i in xrange(n):
        mat2 = perturb(mat, amount)
        est2 = evaluate(mat2)
        if est2 < est:
            mat = mat2
            est = est2

    return mat, est

for i in xrange(1000):
    mat,est = approximate(mat, 1)
    print mat
    print est

.1을 사용한 대략적인 호출은 효과가 없었으므로 제거했습니다. 나는 그것을 잠시 동안 달렸고, 나는 그것이 마지막에 있는지 확인했다.

[[0.7576315397559887, 0, 0.11439449272592839, -0.314856490473439], 
[0.06440497208710227, 1, -0.5607502645413118, 0.38338196981556827], 
[0, 0, 1, 0], 
[0.05421620936883742, 0, -0.5673977598434641, 2.693116299312736]]

0.02 부근의 오류가 있습니다.

1
BBDynSys