Quad warping / homography without opencv


#1

hey hey

just managed to solve how to do the quad warping without using opencv. just pass the original and destination 4 points to the findHomography function and it will return the correct 4x4 matrix to use with glMultMatrix.

the function in opencv can work with more than 4 points that is useful for noisy reads of the points but if you use only 4 points this works the same and its indeed way faster since there’s no need to convert from ofPoint to cvPoint and cv memory structures.

there’s also a function in the file that can be useful for other purposes, a linear equation solver for any number of terms using gaussian elimination.

homography.h.zip


Quad warping point to point
#2

oh sweet, I’m curious to see this in action.

is this a bug?

void findHomography(ofPoint src[4], ofPoint dst[4], float homography[9]){

for(int i=0;i<16;i++) homography[i] = aux_H[i];
}

take care,
zach


#3

oh, yes, at first i was using a 3x3 matrix and forgot to update that. have updated the attachment.


#4

very cool!
that is really great - it felt ridiculous to include OpenCv just for that function.

Wanna cross-post to here - so people know to use your function?
http://forum.openframeworks.cc/t/quad-warping-an-entire-opengl-view-solved/509/24


#5

sure!

btw. have added another function to get the transformation matrix as an ofxMatrix4x4 so it can be used also to transform individual vectors instead of the whole gl matrix:

  
ofxMatrix4x4 H = findHomography(src,dst);  
ofxVec2f v(x,y);  
ofxVec2f t = v*H;  
  

where t is v transformed to be the same as if you were drawing it after using glMultMatrix


quad warping an entire OpenGL view (solved)
homography for camera calibration
quad warping an entire OpenGL view (solved)
#6

Hi!!

I’ve been trying this for a project that i’m working on. The thing is that i need to warp the openframeworks screen. But not only in the corners (4 points), i want to warp it across a grid 4x4 for example (16 regions).

I tried using this like that:

  
  
ofxDraggable rectDragger[5][5];  
ofPoint src[16][4];  
ofPoint dst[16][4];  
GLfloat myMatrix[16];  
ofRectangle rects[16];  
ofxMatrix4x4 matrix;  
  
for(int i=0;i<4;i++){  
		for(int j=0;j<4;j++){  
				ofNoFill();  
				ofSetColor(0xC5FF63);  
				ofSetPolyMode(OF_POLY_WINDING_NONZERO);  
				ofBeginShape();  
					ofVertex(rectDragger[i][j].x, rectDragger[i][j].y);  
					dst[j+(i*4)][0].set(rectDragger[i][j].x, rectDragger[i][j].y, 0);  
					  
					ofVertex(rectDragger[i][j+1].x, rectDragger[i][j+1].y);  
					dst[j+(i*4)][1].set(rectDragger[i][j+1].x, rectDragger[i][j+1].y, 0);  
					  
					ofVertex(rectDragger[i+1][j+1].x, rectDragger[i+1][j+1].y);  
					dst[j+(i*4)][2].set(rectDragger[i+1][j+1].x, rectDragger[i+1][j+1].y, 0);  
					  
					ofVertex(rectDragger[i+1][j].x, rectDragger[i+1][j].y);  
					dst[j+(i*4)][3].set(rectDragger[i+1][j].x, rectDragger[i+1][j].y, 0);  
					  
				ofEndShape(true);  
			  
				homograf.findHomography(src[j+(i*4)], dst[j+(i*4)],myMatrix);  
				ofFill();  
				ofSetColor(0xFF7D4F);  
				glPushMatrix();	  
				glMultMatrixf(myMatrix);  
				glPushMatrix();  
					glTranslatef(rects[j+(i*4)].x+40,rects[j+(i*4)].y+40,0);  
					ofCircle(0,0, 40);  
				glPopMatrix();  
				glPopMatrix();  
				  
		}  
	}  
  
  

but i’m confused how to use this to draw what is showing the OF screen.

thanks!


#7

Playing a little bit I came up with this.

[attachment=0:1qr6s8fs]Captura de pantalla 2010-05-14 a las 17.21.52.png[/attachment:1qr6s8fs]

the only problem that I have now is that in the neighborhoods seems the image does have correspondence with the other cells, as you could see for example in the knuckles.

  
  
class Warper : public ofBaseApp {  
	public:  
		ofxDraggable **rectDragger;  
		ofPoint (*src)[4];  
		ofPoint (*dst)[4];  
		GLfloat myMatrix[16];  
		ofRectangle * rects;  
		ofTexture *partitions;  
		int numR,numC;  
		int texW,texH;   
		int index;  
		homography homograf;  
	  
	Warper(){  
		index = 0;  
		numC = numR = 0;  
	}  
	void init(int numRows,int numColumns,int posX,int posY,int totalW,int totalH){  
		numR = numRows;  
		numC = numColumns;  
		  
		rects = new ofRectangle[numRows*numColumns];  
		src = new ofPoint[numRows*numColumns][4];  
		dst = new ofPoint[numRows*numColumns][4];  
		partitions = new ofTexture[numRows*numColumns];  
		  
		rectDragger = new ofxDraggable*[(numRows+1)];  
		for(int i=0;i<numRows+1;i++){  
			rectDragger[i] = new ofxDraggable[(numColumns+1)];  
		}  
		  
		int offSetX = totalW/numColumns;  
		int offSetY = totalH/numRows;  
		  
		  
		for(int i=0;i<numRows+1;i++){  
			for(int j=0;j<numColumns+1;j++){  
				rectDragger[i][j].setSize(20,20);  
				rectDragger[i][j].setPos((i*(totalW/numColumns))+posX,(j*(totalH/numRows))+posY);  
				rectDragger[i][j].ident = index;  
				rectDragger[i][j].enableMouseEvents();  
				rectDragger[i][j].enableKeyEvents();  
				index++;  
				cout << index << endl;  
			}  
		}  
				  
		  
		for(int i=0;i<numRows;i++){  
			for(int j=0;j<numColumns;j++){  
				  
				rects[j+(i*numR)].x = rectDragger[i][j].x;  
				rects[j+(i*numR)].y = rectDragger[i][j].y;  
				  
				src[j+(i*numR)][0].set(rectDragger[i][j].x,rectDragger[i][j].y,0);  
				src[j+(i*numR)][1].set(rectDragger[i][j+1].x,rectDragger[i][j+1].y,0);  
				src[j+(i*numR)][2].set(rectDragger[i+1][j+1].x,rectDragger[i+1][j+1].y,0);	  
				src[j+(i*numR)][3].set(rectDragger[i+1][j].x,rectDragger[i+1][j].y,0);  
				  
			}  
		}  
		texW = totalW/numColumns;  
		texH = totalH/numRows;  
		for(int i=0;i<numRows*numColumns;i++){	  
			  
			partitions[i].allocate(texW,texH, GL_RGB);  
		}  
	};  
	  
	void draw(){  
		glEnable(GL_TEXTURE_2D);  
		for (int i=0; i<16; i++) {  
			partitions[i].loadScreenData(((i/numR)*(texW)), (i%numC)*(texH), (texW), (texH));  
		}  
		  
		ofSetColor(0x000000);  
		ofRect(0, 0, ofGetWidth(),ofGetHeight());  
		for(int i=0;i<numR;i++){  
			for(int j=0;j<numC;j++){  
					  
					ofNoFill();  
					ofSetColor(0xC5FF63);  
					ofSetPolyMode(OF_POLY_WINDING_NONZERO);  
					ofBeginShape();  
				  
					ofVertex(rectDragger[i][j].x, rectDragger[i][j].y);  
					ofVertex(rectDragger[i][j+1].x, rectDragger[i][j+1].y);  
					ofVertex(rectDragger[i+1][j+1].x, rectDragger[i+1][j+1].y);  
					ofVertex(rectDragger[i+1][j].x, rectDragger[i+1][j].y);  
					ofEndShape(true);  
  
				dst[j+(i*numR)][0].set(rectDragger[i][j].x, rectDragger[i][j].y, 0);  
				dst[j+(i*numR)][1].set(rectDragger[i][j+1].x, rectDragger[i][j+1].y, 0);  
				dst[j+(i*numR)][2].set(rectDragger[i+1][j+1].x, rectDragger[i+1][j+1].y, 0);  
				dst[j+(i*numR)][3].set(rectDragger[i+1][j].x, rectDragger[i+1][j].y, 0);  
				  
				homograf.findHomography(src[j+(i*numR)], dst[j+(i*numR)],myMatrix);  
				ofFill();  
				ofSetColor(0xffffff);  
				glPushMatrix();	  
					glMultMatrixf(myMatrix);  
					glTranslatef(rects[j+(i*numR)].x,rects[j+(i*numR)].y,0);  
					partitions[j+(i*numR)].draw(0,0);  
				glPopMatrix();  
				  
			}  
		}  
		  
	}  
	  
};  
  

Anyone has some suggestions?

thanks!

![](http://forum.openframeworks.cc/uploads/default/727/Captura de pantalla 2010-05-14 a las 17.21.52.png)


Grid / Mesh warping
#8

you could subdivide your mesh more finely.

alternatively, you could do pixel warping instead, but that’s very slow.


#9

I have done something like this. The basic process is to draw to an FBO and then map that FBO to a bezier shape in OpenGL. The bezier shape can then be warped to any size/shape you would like by editing individual points. I used this method to do a pretty complicated interactive wall on a curved screen, so it works pretty well.

Here is my post, or you could say journal entry :-), on how I did it: http://forum.openframeworks.cc/t/ofxfbo-texture-mapping-–-cylinder/4002/0 skip to the bottom to get the best code example.

Also, if you are using Alpha blending, you have to change the blending mode depending on if you are drawing to FBO, or mapping the texture to the bezier shape. see http://forum.openframeworks.cc/t/fbo-problems-with-alpha/1643/10 for details.

DISCLAIMER: I am very much a beginner at openGL, so the code is probably a bit sloppy, but works. (As we say in the South, Redneck programmin’: Don’t know why it works, but it does!)

*edit* Here is a screen shot of what can be accomplished. The pink circles are control points of the mesh.


#10

blerg
i’ve been hitting my head against some bricks for 20hours now…

I keep getting duplicate symbol errors with this file:

ld: duplicate symbol homography::findHomography(ofPoint*, ofPoint*, float*)in /Users/elliot/openFrameworks 0.6/apps/elliots/bricks pc/build/bricks pc.build/Debug/bricks pc.build/Objects-normal/i386/testApp.o and /Users/elliot/openFrameworks 0.6/apps/elliots/bricks pc/build/bricks pc.build/Debug/bricks pc.build/Objects-normal/i386/main.o
collect2: ld returned 1 exit status
Command /Developer/usr/bin/g+±4.2 failed with exit code 1

I tried putting it in a class to see if that got around it. no luck. then i tried Xcode 4 (different compiler). Still same error.

brick wall. brick wall. smack…smack…

think i might try the opencv approach


#11

the .h file you are using elliot needs a include gaurd, since it’s being included multiple times and you are getting a multiple definition error. Try adding:

#pragma once

at the top or,

#ifndef _H_HOMOGRAPHY
#define _H_HOMOGRAPHY

// content

#endif

at the top and bottom.

that prevents it from being included multiple times.

or just include it from one .cpp file.

hope that helps,
zach


#12

have edited the original file to add the include guard


#13

hi

i am using the matrix warping successfully.
i am noticing that my memory use slowly goes up.
would it make sense to do something like this at the end of the warping:

  
if(matrix != NULL) delete matrix;  
	if(src_mat != NULL) delete src_mat;  
	if(dst_mat != NULL) delete dst_mat;  
	if(translate != NULL) delete translate;  

i am not too familiar with that so please excuse my lack of knowledge.

stephan.


Project tracker position back to screen fixing the distortion
#14

Hi - thanks very much for sharing this code - just what I was looking for - here’s the simplest working example I could make - you click four times to place the corners for the destination of the image - see below.

I’ve been playing with projection surfaces and using this for a soft keystoning. What I’d really like now is to use the excellent work Kyle McDonald and others have been doing with Structured Light techniques to discover the matrix automatically - tricky I know!

Cheers,
Dave

  
  
---testApp.h  
  
#ifndef _TEST_APP  
#define _TEST_APP  
  
  
#include "ofMain.h"  
  
#include "ofxVectorMath.h"  
  
class testApp : public ofBaseApp{  
  
	public:  
		void setup();  
		void draw();  
  
		void mousePressed(int x, int y, int button);  
		void findHomography(ofPoint src[4], ofPoint dst[4], float homography[16]);  
		void gaussian_elimination(float *input, int n);  
	  
		ofImage image;  
		ofPoint des[4];  
		int clickCount;  
};  
  
#endif  
  
---testApp.cpp  
  
#include "testApp.h"  
  
void testApp::setup(){  
	clickCount=0;  
	image.loadImage("image.jpg");  
}  
  
void testApp::draw(){  
	if(clickCount==4){  
		ofPoint src[]={ofPoint(0,0),ofPoint(image.width,0),ofPoint(image.width,image.height),ofPoint(0,image.height)};  
		GLfloat matrix[16];  
		findHomography(src,des,matrix);  
		  
		glPushMatrix();     
		glMultMatrixf(matrix);  
		glPushMatrix();  
		image.draw(0,0);  
		glPopMatrix();  
		glPopMatrix();  
	}  
}  
  
void testApp::mousePressed(int x, int y, int button){  
	//click to place the destination of corners of the image - starting in the top left and moving clockwise  
	if(clickCount<4) {  
		des[clickCount]=ofPoint(x,y);  
		++clickCount;  
	}  
}  
  
void testApp::findHomography(ofPoint src[4], ofPoint dst[4], float homography[16]){  
	// arturo castro - 08/01/2010  
	//  
	// create the equation system to be solved  
	//  
	// from: Multiple View Geometry in Computer Vision 2ed  
	//       Hartley R. and Zisserman A.  
	//  
	// x' = xH  
	// where H is the homography: a 3 by 3 matrix  
	// that transformed to inhomogeneous coordinates for each point  
	// gives the following equations for each point:  
	//  
	// x' * (h31*x + h32*y + h33) = h11*x + h12*y + h13  
	// y' * (h31*x + h32*y + h33) = h21*x + h22*y + h23  
	//  
	// as the homography is scale independent we can let h33 be 1 (indeed any of the terms)  
	// so for 4 points we have 8 equations for 8 terms to solve: h11 - h32  
	// after ordering the terms it gives the following matrix  
	// that can be solved with gaussian elimination:  
	  
	float P[8][9]={  
		{-src[0].x, -src[0].y, -1,   0,   0,  0, src[0].x*dst[0].x, src[0].y*dst[0].x, -dst[0].x }, // h11  
		{  0,   0,  0, -src[0].x, -src[0].y, -1, src[0].x*dst[0].y, src[0].y*dst[0].y, -dst[0].y }, // h12  
		  
		{-src[1].x, -src[1].y, -1,   0,   0,  0, src[1].x*dst[1].x, src[1].y*dst[1].x, -dst[1].x }, // h13  
		{  0,   0,  0, -src[1].x, -src[1].y, -1, src[1].x*dst[1].y, src[1].y*dst[1].y, -dst[1].y }, // h21  
		  
		{-src[2].x, -src[2].y, -1,   0,   0,  0, src[2].x*dst[2].x, src[2].y*dst[2].x, -dst[2].x }, // h22  
		{  0,   0,  0, -src[2].x, -src[2].y, -1, src[2].x*dst[2].y, src[2].y*dst[2].y, -dst[2].y }, // h23  
		  
		{-src[3].x, -src[3].y, -1,   0,   0,  0, src[3].x*dst[3].x, src[3].y*dst[3].x, -dst[3].x }, // h31  
		{  0,   0,  0, -src[3].x, -src[3].y, -1, src[3].x*dst[3].y, src[3].y*dst[3].y, -dst[3].y }, // h32  
	};  
	  
	gaussian_elimination(&P[0][0],9);  
	  
	// gaussian elimination gives the results of the equation system  
	// in the last column of the original matrix.  
	// opengl needs the transposed 4x4 matrix:  
	float aux_H[]={ P[0][8],P[3][8],0,P[6][8], // h11  h21 0 h31  
		P[1][8],P[4][8],0,P[7][8], // h12  h22 0 h32  
		0      ,      0,0,0,       // 0    0   0 0  
		P[2][8],P[5][8],0,1};      // h13  h23 0 h33  
	  
	for(int i=0;i<16;i++) homography[i] = aux_H[i];  
}  
  
void testApp::gaussian_elimination(float *input, int n){  
	// arturo castro - 08/01/2010  
	//  
	// ported to c from pseudocode in  
	// [http://en.wikipedia.org/wiki/Gaussian-elimination](http://en.wikipedia.org/wiki/Gaussian-elimination)  
	  
	float * A = input;  
	int i = 0;  
	int j = 0;  
	int m = n-1;  
	while (i < m && j < n){  
		// Find pivot in column j, starting in row i:  
		int maxi = i;  
		for(int k = i+1; k<m; k++){  
			if(fabs(A[k*n+j]) > fabs(A[maxi*n+j])){  
				maxi = k;  
			}  
		}  
		if (A[maxi*n+j] != 0){  
			//swap rows i and maxi, but do not change the value of i  
			if(i!=maxi)  
				for(int k=0;k<n;k++){  
					float aux = A[i*n+k];  
					A[i*n+k]=A[maxi*n+k];  
					A[maxi*n+k]=aux;  
				}  
			//Now A[i,j] will contain the old value of A[maxi,j].  
			//divide each entry in row i by A[i,j]  
			float A_ij=A[i*n+j];  
			for(int k=0;k<n;k++){  
				A[i*n+k]/=A_ij;  
			}  
			//Now A[i,j] will have the value 1.  
			for(int u = i+1; u< m; u++){  
				//subtract A[u,j] * row i from row u  
				float A_uj = A[u*n+j];  
				for(int k=0;k<n;k++){  
					A[u*n+k]-=A_uj*A[i*n+k];  
				}  
				//Now A[u,j] will be 0, since A[u,j] - A[i,j] * A[u,j] = A[u,j] - 1 * A[u,j] = 0.  
			}  
			  
			i++;  
		}  
		j++;  
	}  
	  
	//back substitution  
	for(int i=m-2;i>=0;i--){  
		for(int j=i+1;j<n-1;j++){  
			A[i*n+m]-=A[i*n+j]*A[j*n+m];  
			//A[i*n+j]=0;  
		}  
	}  
}  
  


#15

@zach
seems so simple. i cant think why that didn’t cross my mind.
i was (as may be obvious from the above wording) a bit paniced at the time and was starting to overlook things a little… :smiley:
i went for cv in the end as it was only 1 matrix that needed calculating

I find homography to be such a useful feature (when in vvvv, i’m using it in almost every project), that i’d put a vote in for moving this into ofxVectorMath or something else that comes with oF.

again thanks for the hint
will try this out (with the guarded code) when xcode finishes reinstalling
(upgraded iPad, didn’t know i now have to upgrade xcode to match, which is several hours out…)


#16

ok, xcode back on
yes i had tried that
and still got the error

my current arrangement (though before it was just all in testapp) is:

app>class A>class B>class C>include homography
class C calls homography

get duplicate symbol on gaussian

will figure it out this time


#17

p.s. i also have to include ofxVectorMath, else i cant use ofxMatrix4x4
maybe that’s a hint to where I’m going wrong?


#18

yep, wrapped it up and now that issues gone
will post when properly tested…
not sure where that was coming from

anyway, making a class with static functions


#19

i’m probably the only one retarded enough to need this,
but here it is wrapped up

the functions are static so you can call

  
  
ofxHomographyHelper::findHomography(verticies_source, verticies_destination, your_ofxMatrix4x4.getPtr());  
  

without instantiating the class

http://kimchiandchips.googlecode.com/files/ofxHomographyHelper.zip


Homography on points
#20

I’m completely new to C++ and OpenFrameworks, and I was trying to experiment a bit with this interesting quad-warp examples I’ve found here on the topic.
I’m playing with several-quads setup, for something that could be used for projection-mapping experiments.
I have a quad class that seems to work, and I can initialize several quads independently, select an ‘active’ quad and warp and modify its vertices position. I even can add new quads on-the-fly (I’m using a std::vector class to store them) but I encounter problems if I try to use a keypress event to delete an existing quad and I get a segfault. I suspect it has to do with matrix push-and-pop sequence, but I was not able to debug the code due to my total lack of experience.
I was wondering if someone more experienced and skilled could have a look at my (rather naive) code and give me a hint or tip.

best to all, francesco

------- testApp.h

  
#ifndef _TEST_APP  
#define _TEST_APP  
  
#include "ofMain.h"  
#include "ofAddons.h"  
#include <vector>  
  
  
//stupid ball class  
class ball{  
  
	public:  
		ball(){  
  
		}  
  
		void setup(float x, float y, float r){  
			pos.x = x;  
			pos.y = y;  
			pos.z = r;  
  
			vel.x = 2.0;  
			vel.y = 1.6;  
		}  
  
		void update(int width, int height){  
			if( pos.x + pos.z >= width){  
				pos.x = width - pos.z;  
				vel.x *= -1;  
			}else if( pos.x - pos.z <= 0){  
				pos.x = pos.z;  
				vel.x *= -1;  
			}  
  
			if( pos.y + pos.z >= height){  
				pos.y = height - pos.z;  
				vel.y *= -1;  
			}else if( pos.y - pos.z <= 0){  
				pos.y = pos.z;  
				vel.y *= -1;  
			}  
  
			pos.x += vel.x;  
			pos.y += vel.y;  
		}  
  
		void draw(){  
			ofSetRectMode(OF_RECTMODE_CENTER);  
			ofCircle(pos.x, pos.y,  pos.z);  
			ofSetRectMode(OF_RECTMODE_CORNER);  
		}  
  
		ofPoint pos;  
		ofPoint vel;  
  
};  
  
  
class quad{  
  
	public:  
		quad(){  
		}  
  
		ofPoint corners[4];  
  
		ofPoint src[4];  
		ofPoint dst[4];  
		//lets make a matrix for openGL  
		//this will be the matrix that peforms the transformation  
		GLfloat matrix[16];  
		ofTrueTypeFont ttf;  
		ofTrueTypeFont ttf2;  
		ofImage img;  
		ball balls[80];  
		int borderColor;  
		bool initialized;  
  
  
		void setup(float x1, float y1, float x2, float y2, float x3, float y3, float x4, float y4) {  
  
            initialized = True;  
			//loads load in some truetype fonts  
			ttf.loadFont("type/frabk.ttf", 22);  
			ttf2.loadFont("type/frabk.ttf", 14);  
			//lets load a test image too  
			img.loadImage("car.jpg");  
  
  
				//this is just for our gui / mouse handles  
			//this will end up being the destination quad we are warping too  
			//corners[0].x = 0.0;  
			//corners[0].y = 0.0;  
  
			//corners[1].x = 1.0;  
			//corners[1].y = 0.0;  
  
			//corners[2].x = 1.0;  
			//corners[2].y = 1.0;  
  
			//corners[3].x = 0.0;  
			//corners[3].y = 1.0;  
  
            corners[0].x = x1;  
			corners[0].y = y1;  
  
			corners[1].x = x2;  
			corners[1].y = y2;  
  
			corners[2].x = x3;  
			corners[2].y = y3;  
  
			corners[3].x = x4;  
			corners[3].y = y4;  
  
  
			//lets setup some stupid particles  
			for(int i = 0; i < 80; i++){  
				balls[i].setup(ofRandom(10, ofGetWidth() - 10), ofRandom(10, ofGetHeight()-10), ofRandom(5, 25));  
				balls[i].vel.x = ofRandom(1.5, 2.8);  
				balls[i].vel.y = ofRandom(1.5, 2.8);  
			}  
		}  
  
		void update() {  
  
            borderColor = 0xFF6600;  
			for(int i = 0; i < 80; i++){  
				balls[i].update(ofGetWidth(), ofGetHeight());  
			}  
  
  
//we set matrix to the default - 0 translation  
			//and 1.0 scale for x y z and w  
			for(int i = 0; i < 16; i++){  
				if(i % 5 != 0) matrix[i] = 0.0;  
				else matrix[i] = 1.0;  
			}  
  
			//we set the warp coordinates  
			//source coordinates as the dimensions of our window  
			src[0].x = 0;  
			src[0].y = 0;  
			src[1].x = ofGetWidth();  
			src[1].y = 0;  
			src[2].x = ofGetWidth();  
			src[2].y = ofGetHeight();  
			src[3].x = 0;  
			src[3].y = ofGetHeight();  
  
			//corners are in 0.0 - 1.0 range  
			//so we scale up so that they are at the window's scale  
			for(int i = 0; i < 4; i++){  
				dst[i].x = corners[i].x  * (float)ofGetWidth();  
				dst[i].y = corners[i].y * (float) ofGetHeight();  
			}  
  
  
  
		}  
  
		void draw() {  
			ofPushMatrix();  
  
  
  
		    	// find transformation matrix  
    			findHomography(src, dst, matrix);  
  
  
				//finally lets multiply our matrix  
				//wooooo hoooo!  
				glMultMatrixf(matrix);  
  
  
				// -- NOW LETS DRAW!!!!!!  -----  
  
				//test an image  
				ofSetColor(0xAAAAAA);  
				img.draw(20, 60);  
  
				//lets draw a bounding box  
				ofNoFill();  
				ofSetColor(borderColor);  
				ofRect(1, 1, ofGetWidth()-2, ofGetHeight()-2);  
  
				//our particles  
 				ofEnableAlphaBlending();  
				ofSetColor(255, 120, 0, 130);  
				ofFill();  
				for(int i = 0; i < 40; i++)balls[i].draw();  
 				ofDisableAlphaBlending();  
  
				//some text  
				ofSetColor(0x000000);  
				ttf.drawString("grab corners to warp openGL graphics", 28, 33);  
  
				ofSetColor(0xFFFFFF);  
				ttf.drawString("grab corners to warp openGL graphics", 30, 30);  
  
				ofSetColor(0x000000);  
				ttf2.drawString("warps images nicely too!", 558, 533);  
  
				ofSetColor(0xFF6600);  
				ttf2.drawString("warps images nicely too!", 560, 530);  
  
				ofSetColor(0x000000);  
				ttf2.drawString("press 'r' to reset", 558, 558);  
  
        			ofSetColor(0xFFFFFF);  
				ttf2.drawString("press 'r' to reset", 560, 555);  
  
				ofPopMatrix();  
  
		}  
  
  
  
void gaussian_elimination(float *input, int n){  
	// ported to c from pseudocode in  
	// [http://en.wikipedia.org/wiki/Gaussian-elimination](http://en.wikipedia.org/wiki/Gaussian-elimination)  
  
	float * A = input;  
	int i = 0;  
	int j = 0;  
	int m = n-1;  
	while (i < m && j < n){  
	  // Find pivot in column j, starting in row i:  
	  int maxi = i;  
	  for(int k = i+1; k<m; k++){  
	    if(fabs(A[k*n+j]) > fabs(A[maxi*n+j])){  
	      maxi = k;  
	    }  
	  }  
	  if (A[maxi*n+j] != 0){  
	    //swap rows i and maxi, but do not change the value of i  
		if(i!=maxi)  
		for(int k=0;k<n;k++){  
			float aux = A[i*n+k];  
			A[i*n+k]=A[maxi*n+k];  
			A[maxi*n+k]=aux;  
		}  
	    //Now A[i,j] will contain the old value of A[maxi,j].  
	    //divide each entry in row i by A[i,j]  
		float A_ij=A[i*n+j];  
		for(int k=0;k<n;k++){  
			A[i*n+k]/=A_ij;  
		}  
	    //Now A[i,j] will have the value 1.  
	    for(int u = i+1; u< m; u++){  
    		//subtract A[u,j] * row i from row u  
	    	float A_uj = A[u*n+j];  
	    	for(int k=0;k<n;k++){  
	    		A[u*n+k]-=A_uj*A[i*n+k];  
	    	}  
	      //Now A[u,j] will be 0, since A[u,j] - A[i,j] * A[u,j] = A[u,j] - 1 * A[u,j] = 0.  
	    }  
  
	    i++;  
	  }  
	  j++;  
	}  
  
	//back substitution  
	for(int i=m-2;i>=0;i--){  
		for(int j=i+1;j<n-1;j++){  
			A[i*n+m]-=A[i*n+j]*A[j*n+m];  
			//A[i*n+j]=0;  
		}  
	}  
}  
  
void findHomography(ofPoint src[4], ofPoint dst[4], float homography[16]){  
  
	// create the equation system to be solved  
	//  
	// from: Multiple View Geometry in Computer Vision 2ed  
	//       Hartley R. and Zisserman A.  
	//  
	// x' = xH  
	// where H is the homography: a 3 by 3 matrix  
	// that transformed to inhomogeneous coordinates for each point  
	// gives the following equations for each point:  
	//  
	// x' * (h31*x + h32*y + h33) = h11*x + h12*y + h13  
	// y' * (h31*x + h32*y + h33) = h21*x + h22*y + h23  
	//  
	// as the homography is scale independent we can let h33 be 1 (indeed any of the terms)  
	// so for 4 points we have 8 equations for 8 terms to solve: h11 - h32  
	// after ordering the terms it gives the following matrix  
	// that can be solved with gaussian elimination:  
  
	float P[8][9]={  
			{-src[0].x, -src[0].y, -1,   0,   0,  0, src[0].x*dst[0].x, src[0].y*dst[0].x, -dst[0].x }, // h11  
			{  0,   0,  0, -src[0].x, -src[0].y, -1, src[0].x*dst[0].y, src[0].y*dst[0].y, -dst[0].y }, // h12  
  
			{-src[1].x, -src[1].y, -1,   0,   0,  0, src[1].x*dst[1].x, src[1].y*dst[1].x, -dst[1].x }, // h13  
		    {  0,   0,  0, -src[1].x, -src[1].y, -1, src[1].x*dst[1].y, src[1].y*dst[1].y, -dst[1].y }, // h21  
  
			{-src[2].x, -src[2].y, -1,   0,   0,  0, src[2].x*dst[2].x, src[2].y*dst[2].x, -dst[2].x }, // h22  
		    {  0,   0,  0, -src[2].x, -src[2].y, -1, src[2].x*dst[2].y, src[2].y*dst[2].y, -dst[2].y }, // h23  
  
			{-src[3].x, -src[3].y, -1,   0,   0,  0, src[3].x*dst[3].x, src[3].y*dst[3].x, -dst[3].x }, // h31  
		    {  0,   0,  0, -src[3].x, -src[3].y, -1, src[3].x*dst[3].y, src[3].y*dst[3].y, -dst[3].y }, // h32  
	};  
  
	gaussian_elimination(&P[0][0],9);  
  
	// gaussian elimination gives the results of the equation system  
	// in the last column of the original matrix.  
	// opengl needs the transposed 4x4 matrix:  
	float aux_H[]={ P[0][8],P[3][8],0,P[6][8], // h11  h21 0 h31  
					P[1][8],P[4][8],0,P[7][8], // h12  h22 0 h32  
					0      ,      0,0,0,       // 0    0   0 0  
					P[2][8],P[5][8],0,1};      // h13  h23 0 h33  
  
	for(int i=0;i<16;i++) homography[i] = aux_H[i];  
}  
  
  
  
  
  
  
  
  
};  
  
  
  
  
  
class testApp : public ofSimpleApp{  
  
	public:  
  
		void setup();  
		void update();  
		void draw();  
  
		void keyPressed(int key);  
		void keyReleased(int key);  
		void mouseMoved(int x, int y );  
		void mouseDragged(int x, int y, int button);  
		void mousePressed(int x, int y, int button);  
		void mouseReleased();  
		int whichCorner;  
		ofTrueTypeFont ttf;  
  
        std::vector<quad> quads;  
		//quad quads[10];  
		int nOfQuads;  
		int activeQuad;  
  
};  
  
#endif  

---------------------- testApp.cpp

  
#include "testApp.h"  
#include "stdio.h"  
#include <iostream>  
  
  
//--------------------------------------------------------------  
void testApp::setup(){  
	//we run at 60 fps!  
	ofSetVerticalSync(true);  
  
  
  
    ttf.loadFont("type/frabk.ttf", 11);  
  
	//for(int i = 0; i<10; i++){  
	//quads[i].setup();  
	//}  
	//quads[0].setup(0.0,0.0,0.5,0.0,0.5,0.5,0.0,0.5);  
    //quads[1].setup(0.5,0.0,1.0,0.0,1.0,0.5,0.5,0.5);  
    //quads[2].setup(0.0,0.5,0.5,0.5,0.5,1.0,0.0,1.0);  
    //quads[3].setup(0.5,0.5,1.0,0.5,1.0,1.0,0.5,1.0);  
    quads.reserve(16);  
    quad one;  
    one.setup(0.0,0.0,0.5,0.0,0.5,0.5,0.0,0.5);  
    quad two;  
    two.setup(0.5,0.0,1.0,0.0,1.0,0.5,0.5,0.5);  
    quad three;  
    three.setup(0.0,0.5,0.5,0.5,0.5,1.0,0.0,1.0);  
    quad four;  
    four.setup(0.5,0.5,1.0,0.5,1.0,1.0,0.5,1.0);  
    quads.push_back(one);  
    quads.push_back(two);  
    quads.push_back(three);  
    quads.push_back(four);  
	activeQuad = 0;  
}  
  
//--------------------------------------------------------------  
void testApp::update(){  
    nOfQuads = quads.size();  
	ofBackground(20, 20, 20);  
	ofSetWindowShape(800, 600);  
    for(int i = 0; i < nOfQuads; i++){  
		quads[i].update();  
	}  
}  
  
//--------------------------------------------------------------  
void testApp::draw(){  
  
    nOfQuads = quads.size();  
    quads[activeQuad].borderColor = 0xFFFFFF;  
  
	for(int i = 0; i < nOfQuads; i++){  
        quads[i].draw();  
	}  
  
    ofSetColor(0xFFFFFF);  
    ttf.drawString("active quad: "+ofToString(activeQuad), 30, 570);  
  
}  
  
//--------------------------------------------------------------  
void testApp::keyPressed(int key){  
  
if ( key =='r' || key == 'R' )  
	{  
	quads[activeQuad].corners[0].x = 0.0;  
	quads[activeQuad].corners[0].y = 0.0;  
  
	quads[activeQuad].corners[1].x = 1.0;  
	quads[activeQuad].corners[1].y = 0.0;  
  
	quads[activeQuad].corners[2].x = 1.0;  
	quads[activeQuad].corners[2].y = 1.0;  
  
	quads[activeQuad].corners[3].x = 0.0;  
	quads[activeQuad].corners[3].y = 1.0;  
	}  
  
if ( key =='>' )  
	{  
	activeQuad += 1;  
	if (activeQuad > quads.size()-1) {activeQuad = 0;}  
	}  
  
if ( key =='<' )  
	{  
	activeQuad -= 1;  
	if (activeQuad < 0) {activeQuad = quads.size()-1;}  
	}  
  
if ( key =='a' )  
	{  
    if (quads.size() < 16) {  
	quad i;  
	i.setup(0.25,0.25,0.75,0.25,0.75,0.75,0.25,0.75);  
    quads.push_back(i);  
    activeQuad = quads.size()-1;  
	}  
	}  
  
if ( key =='d' )  
	{  
    if (quads.size() > 1) {  
	quads.erase(quads.begin()+activeQuad);  
    activeQuad -= 1;  
    nOfQuads = quads.size();  
	}  
	}  
  
  
}  
  
//--------------------------------------------------------------  
void testApp::keyReleased(int key){  
}  
  
  
//--------------------------------------------------------------  
void testApp::mouseMoved(int x, int y ){  
  
  
}  
  
//--------------------------------------------------------------  
void testApp::mouseDragged(int x, int y, int button){  
  
		float scaleX = (float)x / ofGetWidth();  
		float scaleY = (float)y / ofGetHeight();  
  
		if(whichCorner >= 0){  
			quads[activeQuad].corners[whichCorner].x = scaleX;  
			quads[activeQuad].corners[whichCorner].y = scaleY;  
		}  
  
  
}  
  
//--------------------------------------------------------------  
void testApp::mousePressed(int x, int y, int button){  
  
	float smallestDist = 1.0;  
	whichCorner = -1;  
  
	for(int i = 0; i < 4; i++){  
		float distx = quads[activeQuad].corners[i].x - (float)x/ofGetWidth();  
		float disty = quads[activeQuad].corners[i].y - (float)y/ofGetHeight();  
		float dist  = sqrt( distx * distx + disty * disty);  
  
		if(dist < smallestDist && dist < 0.1){  
			whichCorner = i;  
			smallestDist = dist;  
		}  
	}  
}  
  
//--------------------------------------------------------------  
void testApp::mouseReleased(){  
	whichCorner = -1;  
}