Michael Vogel wrote:
As I was not able to reverse the 3D–>2D conversion, I actually calculate the 2D coordinates for n×n geo coordinates, check which is the nearest to the mouse position, calculate n×n geo coordinates near this point and so on – not the best approach, but at least working.
I cleaned up this code a little but I think you could probably tune it to your individual case better than I could. It converts between the coordinates of the screen point (sx, sy) and a texture point (u, v). In your case the texture is the geo cordinates.
projection_map.pbiCode:
EnableExplicit
;Projection mapping routines translated from code presented in
;"Fundamentals of texture Mapping and Image Warping" by Paul S. Heckbert for his
;Master's Thesis at the University of California in 1989.
;
;Coded in PureBasic and adapted by Demivec.
Structure POINT_d
x.d
y.d
EndStructure
Structure poly_vert ;a polygon vetex
u.d ;texture space coordinates
v.d
sx.d ;screen space coordinates
sy.d
EndStructure
Structure poly_quad ;a polygon, specifically a quadrilateral
vert.poly_vert[4] ;vertices
EndStructure
;pmap return codes
Enumeration -1
#PMAP_BAD ;mapping doesn't exist (bad correspondence)
#PMAP_AFFINE ;mapping is affine
#PMAP_PROJECTIVE ;mapping is projective
EndEnumeration
;-matrix code
;Routines for manipulating 3 x 3 matrices
;These matrices are used to represent a 2-D projective mapping.
;Determinate for a 2x2 matrix.
Macro det2(a, b, c, d)
((a) * (d) - (b) * (c))
EndMacro
Procedure.d mx3d_adjoint(Array a.d(2), Array b.d(2))
;b = adjoint(a), returns determinant(a)
b(0, 0) = det2(a(1, 1), a(1, 2), a(2, 1), a(2, 2))
b(1, 0) = det2(a(1, 2), a(1, 0), a(2, 2), a(2, 0))
b(2, 0) = det2(a(1, 0), a(1, 1), a(2, 0), a(2, 1))
b(0, 1) = det2(a(2, 1), a(2, 2), a(0, 1), a(0, 2))
b(1, 1) = det2(a(2, 2), a(2, 0), a(0, 2), a(0, 0))
b(2, 1) = det2(a(2, 0), a(2, 1), a(0, 0), a(0, 1))
b(0, 2) = det2(a(0, 1), a(0, 2), a(1, 1), a(1, 2))
b(1, 2) = det2(a(0, 2), a(0, 0), a(1, 2), a(1, 0))
b(2, 2) = det2(a(0, 0), a(0, 1), a(1, 0), a(1, 1))
;If matrix is non-invertable the determinant will equal zero.
;The determinatnt may also equal zero due to floating point rounding errors.
ProcedureReturn a(0, 0) * b(0, 0) + a(0, 1) * b(0, 1) + a(0, 2) * b(0, 2)
EndProcedure
Procedure mx3d_mul(Array a.d(2), Array b.d(2), Array c.d(2))
;matrix multiply: c = a * b
c(0, 0) = a(0, 0) * b(0, 0) + a(0, 1) * b(1, 0) + a(0, 2) * b(2, 0)
c(0, 1) = a(0, 0) * b(0, 1) + a(0, 1) * b(1, 1) + a(0, 2) * b(2, 1)
c(0, 2) = a(0, 0) * b(0, 2) + a(0, 1) * b(1, 2) + a(0, 2) * b(2, 2)
c(1, 0) = a(1, 0) * b(0, 0) + a(1, 1) * b(1, 0) + a(1, 2) * b(2, 0)
c(1, 1) = a(1, 0) * b(0, 1) + a(1, 1) * b(1, 1) + a(1, 2) * b(2, 1)
c(1, 2) = a(1, 0) * b(0, 2) + a(1, 1) * b(1, 2) + a(1, 2) * b(2, 2)
c(2, 0) = a(2, 0) * b(0, 0) + a(2, 1) * b(1, 0) + a(2, 2) * b(2, 0)
c(2, 1) = a(2, 0) * b(0, 1) + a(2, 1) * b(1, 1) + a(2, 2) * b(2, 1)
c(2, 2) = a(2, 0) * b(0, 2) + a(2, 1) * b(1, 2) + a(2, 2) * b(2, 2)
EndProcedure
Procedure mx3d_transform(Array p.d(1), Array a_matrix.d(2), Array q.d(1)) ;the p & q matrices are 1x3
;transform point p by matrix a_Matrix: q = p * a_Matrix
q(0) = p(0) * a_matrix(0, 0) + p(1) * a_matrix(1, 0) + p(2) * a_matrix(2, 0)
q(1) = p(0) * a_matrix(0, 1) + p(1) * a_matrix(1, 1) + p(2) * a_matrix(2, 1)
q(2) = p(0) * a_matrix(0, 2) + p(1) * a_matrix(1, 2) + p(2) * a_matrix(2, 2)
EndProcedure
;-general projective mapping code
;General routines for 2-D projective mappings.
;These routines, pmap_quad_rect and pmap_square_quad,
;are independent of the poly_quad structure, so we do not
;think in terms of texture and screen space.
#TOLERANCE = 1e-13
; Macro zero(x): ((x) < #TOLERANCE And (x) > -#TOLERANCE): EndMacro ;used in conditionals
Procedure pmap_square_quad(Array quadVertex.POINT_d(1), Array sq_matrix.d(2))
;Find mapping between unit square and quadrilateral
;The correspondence is:
;
; (0, 0) --> quadVertex(0)
; (1, 0) --> quadVertex(1)
; (1, 1) --> quadVertex(2)
; (0, 1) --> quadVertex(3)
;quadVertex: vertices of quadrilateral
;sq : square->quad transform (returned)
Protected.d px, py
px = quadVertex(0)\x - quadVertex(1)\x + quadVertex(2)\x - quadVertex(3)\x
py = quadVertex(0)\y - quadVertex(1)\y + quadVertex(2)\y - quadVertex(3)\y
If (px < #TOLERANCE And px > -#TOLERANCE) And (py < #TOLERANCE And py > -#TOLERANCE) ;affine
sq_matrix(0, 0) = quadVertex(1)\x - quadVertex(0)\x
sq_matrix(1, 0) = quadVertex(2)\x - quadVertex(1)\x
sq_matrix(2, 0) = quadVertex(0)\x
sq_matrix(0, 1) = quadVertex(1)\y - quadVertex(0)\y
sq_matrix(1, 1) = quadVertex(2)\y - quadVertex(1)\y
sq_matrix(2, 1) = quadVertex(0)\y
sq_matrix(0, 2) = 0
sq_matrix(1, 2) = 0
sq_matrix(2, 2) = 1
ProcedureReturn #PMAP_AFFINE
Else ;projective
Protected.d dx1, dx2, dy1, dy2, del
dx1 = quadVertex(1)\x - quadVertex(2)\x
dx2 = quadVertex(3)\x - quadVertex(2)\x
dy1 = quadVertex(1)\y - quadVertex(2)\y
dy2 = quadVertex(3)\y - quadVertex(2)\y
del = det2(dx1, dx2, dy1, dy2)
If del = 0
MessageRequester("Error", "pmap_square_quad: bad mapping")
ProcedureReturn #PMAP_BAD
EndIf
sq_matrix(0, 2) = det2(px, dx2, py, dy2) / del
sq_matrix(1, 2) = det2(dx1, px, dy1, py) / del
sq_matrix(2, 2) = 1
sq_matrix(0, 0) = quadVertex(1)\x - quadVertex(0)\x + sq_matrix(0, 2) * quadVertex(1)\x
sq_matrix(1, 0) = quadVertex(3)\x - quadVertex(0)\x + sq_matrix(1, 2) * quadVertex(3)\x
sq_matrix(2, 0) = quadVertex(0)\x
sq_matrix(0, 1) = quadVertex(1)\y - quadVertex(0)\y + sq_matrix(0, 2) * quadVertex(1)\y
sq_matrix(1, 1) = quadVertex(3)\y - quadVertex(0)\y + sq_matrix(1, 2) * quadVertex(3)\y
sq_matrix(2, 1) = quadVertex(0)\y
ProcedureReturn #PMAP_PROJECTIVE
EndIf
EndProcedure
Procedure pmap_quad_rect(u0.d, v0.d, u1.d, v1.d, Array quadVertex.POINT_d(1), Array qr_matrix.d(2))
;Find mapping between quadrilateral and rectangle.
;The correspondence is:
;
; quadVertex(0) --> (u0, v0)
; quadVertex(1) --> (u1, v0)
; quadVertex(2) --> (u1, v1)
; quadVertex(3) --> (u0, v1)
;
;This method of computing the adjoint is numericaly cheaper than
;computing it symbolically.
;u0, v0, u1, v1: bounds of rectangle
;quadVertex(): vertices of quadrilateral
;qr_Matrix(): quad->rect transform (returned)
Protected ret, du.d, dv.d
Protected Dim rq_matrix.d(2, 2) ;rect->quad transform
du = u1 - u0
dv = v1 - v0
If du = 0 Or dv = 0
MessageRequester("Error", "pmap_quad_rect: null rectangle")
ProcedureReturn #PMAP_BAD
EndIf
;first find mapping from unit uv square to xy quadrilateral
ret = pmap_square_quad(quadVertex(), rq_matrix())
If ret = #PMAP_BAD: ProcedureReturn #PMAP_BAD: EndIf
;concatenate transform from uv rectangle (u0, v0, u1, v1) to unit square
rq_matrix(0, 0) / du
rq_matrix(1, 0) / dv
rq_matrix(2, 0) - rq_matrix(0, 0) * u0 + rq_matrix(1, 0) * v0
rq_matrix(0, 1) / du
rq_matrix(1, 1) / dv
rq_matrix(2, 1) - rq_matrix(0, 1) * u0 + rq_matrix(1, 1) * v0
rq_matrix(0, 2) / du
rq_matrix(1, 2) / dv
rq_matrix(2, 2) - rq_matrix(0, 2) * u0 + rq_matrix(1, 2) * v0
;now rq_Matrix is transform from uv rectangle to xy quadrilateral
;qr_Matrix = inverse transform, which maps xy to uv
If mx3d_adjoint(rq_matrix(), qr_matrix()) = 0
MessageRequester("Warning", "pmap_quad_rect: determinate = 0")
EndIf
ProcedureReturn ret
EndProcedure
;-projective mapping code
Procedure pmap_quad_quad(*p.poly_quad, Array st_matrix.d(2))
;Find the projective mapping st from screen to texture space
;given the screen and texture coordinates at the vertices of a polygon (specifically a quadrilateral).
;
;Method: concatenate screen_quad->mid_square and
;mid_square->texture_quad transform matrices.
;The alternative method, solving an 8x8 system of equations, takes longer.
Protected i, type1, type2
Protected Dim quadVertex.POINT_d(3) ;screen_quad vertices
Protected Dim ms_matrix.d(2, 2) ;mid_square->screen transform matrix
Protected Dim sm_matrix.d(2, 2) ;screen->mid transform matrix
Protected Dim mt_matrix.d(2, 2) ;mid->texture transform matrix
;prepare mid_square->screen transform matrix
For i = 0 To 3
quadVertex(i)\x = *p\vert[i]\sx: quadVertex(i)\y = *p\vert[i]\sy
Next
type1 = pmap_square_quad(quadVertex(), ms_matrix())
If mx3d_adjoint(ms_matrix(), sm_matrix()) = 0
MessageRequester("Warning", "pmap_quad_quad: determinant = 0")
EndIf
;prepare mid->texture transform matrix
For i = 0 To 3
quadVertex(i)\x = *p\vert[i]\u: quadVertex(i)\y = *p\vert[i]\v
Next
type2 = pmap_square_quad(quadVertex(), mt_matrix())
If type1 = #PMAP_BAD Or type2 = #PMAP_BAD
ProcedureReturn #PMAP_BAD
EndIf
mx3d_mul(sm_matrix(), mt_matrix(), st_matrix())
ProcedureReturn type1 | type2 ;#PMAP_PROJECTIVE prevails over #PMAP_AFFINE
EndProcedure
Procedure pmap_poly(*p.poly_quad, Array st_matrix.d(2))
;Find the projective mapping st from screen to texture space
;given the screen and texture coordinates at the vertices of a polygon (specifically a quadrilateral),
;which are passed in the fields (sx, sy) and (u, v) of p\vert[i], respectively.
Protected Dim screenQuad.POINT_d(3) ;vertices of screen quadrilateral
With *p
If \vert[0]\v = \vert[1]\v And \vert[2]\v = \vert[3]\v And \vert[1]\u = \vert[2]\u And \vert[3]\u = \vert[0]\u
;edges 0-1 and 2-3 are horz, so 1-2 and 3-0 are vert
screenQuad(0)\x = \vert[0]\sx: screenQuad(0)\y = \vert[0]\sy
screenQuad(1)\x = \vert[1]\sx: screenQuad(1)\y = \vert[1]\sy
screenQuad(2)\x = \vert[2]\sx: screenQuad(2)\y = \vert[2]\sy
screenQuad(3)\x = \vert[3]\sx: screenQuad(3)\y = \vert[3]\sy
ProcedureReturn pmap_quad_rect(\vert[0]\u, \vert[0]\v, \vert[2]\u, \vert[2]\v, screenQuad(), st_matrix())
ElseIf \vert[0]\u = \vert[1]\u And \vert[2]\u = \vert[3]\u And \vert[1]\v = \vert[2]\v And \vert[3]\v = \vert[0]\v
;edges 0-1 and 2-3 are vert, so 1-2 and 3-0 are horz
screenQuad(0)\x = \vert[1]\sx: screenQuad(0)\y = \vert[1]\sy
screenQuad(1)\x = \vert[2]\sx: screenQuad(1)\y = \vert[2]\sy
screenQuad(2)\x = \vert[3]\sx: screenQuad(2)\y = \vert[3]\sy
screenQuad(3)\x = \vert[0]\sx: screenQuad(3)\y = \vert[0]\sy
ProcedureReturn pmap_quad_rect(\vert[1]\u, \vert[1]\v, \vert[3]\u, \vert[3]\v, screenQuad(), st_matrix())
Else
;texture is not an orthogonally-oriented rectangle
ProcedureReturn pmap_quad_quad(*p, st_matrix())
EndIf
EndWith
EndProcedure
DisableExplicit
Test code:
Code:
; Define
#mcolor = $E0E0E0 ;$8A8A8A
#ncolor = $FFD0D0 ;$4C804C
#pcolor = $D0D0FF ;$343434
#W=640
#H=480
XIncludeFile "projection_map.pbi"
Global Dim polygon(7)
Global Dim Real.d(7)
Global screenQuad.poly_quad
Global Dim screenToGeo.d(2, 2)
Global Dim geoToScreen.d(2, 2)
Global Dim screenCoord.d(2)
Global Dim screenCoord2.d(2)
Global Dim geoCoord.d(2)
Global Dim geoCoord2.d(2)
Procedure InitValues()
polygon(0)=40: polygon(1)=40
polygon(2)=340: polygon(3)=30
polygon(4)=550: polygon(5)=300
polygon(6)=100: polygon(7)=440
;test case 1 for another arbitrary polygon
; polygon(0)=66: polygon(1)=39
; polygon(2)=444: polygon(3)=33
; polygon(4)=300: polygon(5)=457
; polygon(6)=17: polygon(7)=440
Real(0)=-100: Real(1)=-100
Real(2)=100: Real(3)=-100
Real(4)=100: Real(5)=100
Real(6)=-100: Real(7)=100
;test case 1
;Other arbitrary rectangles used as geo coordinates, to test robustness of methods.
; Real(0)=-20: Real(1)=-30
; Real(2)=60: Real(3)=-30
; Real(4)=60: Real(5)=40
; Real(6)=-20: Real(7)=40
;test case 2
; Real(0)=-20: Real(1)=30
; Real(2)=60: Real(3)=30
; Real(4)=60: Real(5)=40
; Real(6)=-20: Real(7)=40
;set (sx, sy) for each vertex of screen quadrilateral along with the corresponding (u, v) geo points
With screenQuad
;The following values is used for a reflection transformation to adjust for y axis being flipped
;between the screen and the geo cordinates. The value has only been tested based on a rectangular
;or square quadrilateral that is being used for the geo cordinates.
Protected dy.d = - Real(3) * 2 ;The least y value of the rectangle needs to be used for the reflection adjustment.
\vert[0]\u = Real(0): \vert[0]\v = Real(1) + dy: \vert[0]\sx = polygon(0): \vert[0]\sy = polygon(1)
\vert[1]\u = Real(2): \vert[1]\v = Real(3) + dy: \vert[1]\sx = polygon(2): \vert[1]\sy = polygon(3)
\vert[2]\u = Real(4): \vert[2]\v = Real(5) + dy: \vert[2]\sx = polygon(4): \vert[2]\sy = polygon(5)
\vert[3]\u = Real(6): \vert[3]\v = Real(7) + dy: \vert[3]\sx = polygon(6): \vert[3]\sy = polygon(7)
EndWith
EndProcedure
Procedure.i IsPointInsideTriangle2D(x0.d,y0.d,x1.d,y1.d,x2.d,y2.d,px.d,py.d)
Protected V1x.d=x1-x0
Protected V1y.d=y1-y0
Protected V2x.d=x2-x1
Protected V2y.d=y2-y1
Protected V3x.d=x0-x2
Protected V3y.d=y0-y2
If V1y*(px-x0)>=V1x*(py-y0)
If V2y*(px-x1)>=V2x*(py-y1) And V3y*(px-x2)>=V3x*(py-y2)
ProcedureReturn #True
EndIf
Else
If V2y*(px-x1)<V2x*(py-y1) And V3y*(px-x2)<V3x*(py-y2)
ProcedureReturn #True
EndIf
EndIf
ProcedureReturn #False
EndProcedure
Procedure.i IsPointInsideRectancle2D(x0.d,y0.d,x1.d,y1.d,x2.d,y2.d,x3.d,y3.d,x.d,y.d)
LineXY(x0,y0,x1,y1,#ncolor)
LineXY(x0,y0,x3,y3,#ncolor)
LineXY(x2,y2,x1,y1,#ncolor)
LineXY(x2,y2,x3,y3,#ncolor)
If IsPointInsideTriangle2D(x0,y0,x1,y1,x2,y2,x,y)
ProcedureReturn 1
EndIf
If IsPointInsideTriangle2D(x2,y2,x3,y3,x0,y0,x,y)
ProcedureReturn 2
EndIf
ProcedureReturn 0
EndProcedure
Procedure Cross(x,y,size,color)
LineXY(x-size,y-size,x+size,y+size,color)
LineXY(x+size,y-size,x-size,y+size,color)
EndProcedure
Procedure UpdateImage()
Protected x,y,color
x=WindowMouseX(0)
y=WindowMouseY(0)
;InitValues()
StartDrawing(ImageOutput(0))
Box(0,0,#W,#H,#Black)
DrawingMode(#PB_2DDrawing_Transparent)
For i=0 To 3
LineXY(polygon(i<<1&7),polygon(i<<1&7+1),polygon((i<<1+2)&7),polygon((i<<1+3)&7),#White)
DrawText(polygon(i<<1&7),polygon(i<<1&7+1),Str(Real(i<<1&7))+"|"+Str(Real(i<<1&7+1)),#Yellow)
Next i
;center point
geoCoord(0) = (Real(0) + Real(4)) / 2: geoCoord(1) = (Real(3) + Real(7)) / 2: geoCoord(2) = 1
mx3d_transform(geoCoord(), geoToScreen(), screenCoord())
Cross(screenCoord(0) / screenCoord(2), screenCoord(1) / screenCoord(2), 5, #pcolor)
Select IsPointInsideRectancle2D(polygon(0),polygon(1),polygon(2),polygon(3),polygon(4),polygon(5),polygon(6),polygon(7),x,y)
Case 0
color=#Red
Case 1
color=#Green
Case 2
color=#Yellow
EndSelect
;cross hairs
LineXY(x-29,y,x+29,y,color)
LineXY(x,y-29,x,y+29,color)
If color<>#Red
; x_ and y_ are calculated here...
screenCoord(0) = x: screenCoord(1) = y: screenCoord(2) = 1
mx3d_transform(screenCoord(), screenToGeo(), geoCoord())
Protected.d x_, y_
x_ = geoCoord(0) / geoCoord(2)
y_ = geoCoord(1) / geoCoord(2)
;horizontal line
geoCoord(0) = Real(0): geoCoord(1) = y_: geoCoord(2) = 1
mx3d_transform(geoCoord(), geoToScreen(), screenCoord())
geoCoord2(0) = Real(4): geoCoord2(1) = y_: geoCoord2(2) = 1
mx3d_transform(geoCoord2(), geoToScreen(), screenCoord2())
LineXY(screenCoord(0) / screenCoord(2), screenCoord(1) / screenCoord(2), screenCoord2(0) / screenCoord2(2), screenCoord2(1) / screenCoord2(2), #White)
;vertical line
geoCoord(0) = x_: geoCoord(1) = Real(3): geoCoord(2) = 1
mx3d_transform(geoCoord(), geoToScreen(), screenCoord())
geoCoord2(0) = x_: geoCoord2(1) = Real(7): geoCoord2(2) = 1
mx3d_transform(geoCoord2(), geoToScreen(), screenCoord2())
LineXY(screenCoord(0) / screenCoord(2), screenCoord(1) / screenCoord(2), screenCoord2(0) / screenCoord2(2), screenCoord2(1) / screenCoord2(2), #White)
DrawText(x+10,y-25,Str(x_)+" | "+Str(y_), #Yellow)
EndIf
StopDrawing()
SetGadgetState(0,ImageID(0))
EndProcedure
OpenWindow(0,0,0,#W,#H,"")
ImageGadget(0,0,0,#W,#H,0)
CreateImage(0,#W,#H,24)
InitValues()
;compute screen to geo transform
pmap_poly(screenQuad, screenToGeo())
If Not mx3d_adjoint(screenToGeo(), geoToScreen()): End: EndIf
UpdateImage()
Repeat
Select WaitWindowEvent()
Case #WM_MOUSEMOVE
UpdateImage()
Case #WM_CHAR
End
EndSelect
ForEver
I believe this does what you wanted.
I had one small problem that required some additional work to resolve. I setup two quadrilaterals with the coordinates you gave (Real(), and Polygon()). I then use the routines to create the transform matrices to convert between the two sets of coordinates. The problem that occurs required me to reflect the Real() coordinates to account for the flipping of the Y axis betwen the quadrilaterals defined by the Polygon() and the Real() coordinates. Either set of coordinates could probably have been used and so I chose the Real() ones because they were rectangular and thereby easier to do.
I'm going to see if I can add the reflective transformation, in a separate procedure call, that combines it as a matrix with the transformation matrix that has already been produced. This would allow referencing the original coordinates without them becoming distorted by the way the reflection transformation is in the current code.
This change distorts the copy of the coordinates because the reflected Y values are not the same as the originals. This means if you want to reference them, such as for displaying them for instance, you need to look to another copy of them (i.e. the ones stored in Real()). I do this when positioning things on the screen polygon based on the geo coordinates when displaying the center point and lines traversing the polygon based on the mouse position.
I did enough testing to recommend this as a good solution. I am going to explore other useful aspects of the projection mapping routines and possibly post them in the Tips & Tricks forum because of their usefullness.
@Edit: updated code, and explanatory test
FAQ
Search
Register
Login
