[Home]TotallyUnrelatedMathQuestions

Robo Home | Changes | Preferences | AllPages

Question

Thought I'd pay a little visit back to the RoboWiki in an unashamed attempt to recruit some help for a little project I have going. It's a little hard to explain the problem, but I think I'm up to the task. Say you have a rectangular piece of elastic with a dot on it, and you stretch the four corners of the rectangle out different directions. Given the coordinates of the 4 corners and the point, I need to figure the coordinates of the point when the rectangle is unstretched. Any ideas? --Simonton

Well... I might be able to help but I have a couple questions to clarify the problem. Are you assuming the sides of the elastic rectangle stay straight unlike a truly realistic elastic which becomes thinner when stretched? I'm guessing you assume stretching is linear as well? I think I have an understanding of how to do this if the answer is yes to both of those :) -- Rednaxela

The answer is yes to both, though the latter is a simplifying assumption ;). --Simonton

Alright, basically what you're looking for, is a way to determine the [affine transformation matrix] (an affine transformation is a transformation that essentially fits the same constraints that this simplified case of elastic stretching has) corresponding to known distortion of 4 known points. From that matrix one can easily determine the distorted location of any other point. I'll work up a nice solution to the 2-dimensional (I assume 2-dimensional is what you're wanting) form of the problem, which should be understandable without needing a background with linear algebra, tomorrow. Tonight I am tired :) -- Rednaxela

I'm not sure this problem can be described by an affine transform. If the points can be moved in arbitrary directions by arbitrary distances, then I'd say I'm sure its not an affine transform. The points can move in this way while still allowing the stretching of the material to be linear between the points. But even if the transform isn't affine, as long as the shape is still convex it should be possible to find a solution.

Label one point O for origin. Label the other 3 A, B, C, with B being the opposite corner to the origin. The point in the stretched shape is P. We have 2 parameters u and v which describe the relative distance along the edges. P should be on the line between the points (O + u.OA) and (C + u.CB), and also on the line between the points (O + v.OC) and (A + v.AB). All we have to do now is work out and solve the vector equations to get the values of u and v. :) -- Tim Foden

(APx - v.ABx).[AOy + v(OCy-ABy)] - (APy - v.ABy).[AOx + v(OCx-ABx)] = 0

APx.AOy + v.APx.(OCy-ABy) - v.ABx.AOy - v.v.ABx.(OCy-ABy) - APy.AOx - v.APy.(OCx-ABx) + v.ABy.AOx + v.v.ABy.(OCx-ABx) = 0

v.v [ABy.(OCx-ABx) - ABx.(OCy-ABy)] + v [APx.(OCy-ABy) - ABx.AOy - APy.(OCx-ABx) + ABy.AOx] + [APx.AOy - APy.AOx] = 0

Let a = ABy.(OCx-ABx) - ABx.(OCy-ABy)
Let b = APx.(OCy-ABy) - ABx.AOy - APy.(OCx-ABx) + ABy.AOx
Let c = APx.AOy - APy.AOx

So, v = [-b +/- sqrt(b.b - 4.a.c)] / 2.a, where v is also restricted by 0 <= v <= 1.

Actually, yeah, this can't quite be described as an affine transform, due to the cases where one side may be stretched but another side not. The problem though is still very similar to an affine transform I believe though. In 2-dimensions, linear transforms have a 2x2 matrix that maps from [x1,y1] to [x2,y2], the 4 unknowns in the matrix can be solved for by knowing two points around an origin. Affine transformations are similar, but support translations as well, can be described as a linear transformation that maps from [x1, y1, 1] to [x2, y2], and the 6 unknowns in the 3x2 matrix describing the transform can be solved for by knowing 3 points. The transforms we're talking about here also support stretching of one side with compression of another side, and I'm 99% sure can be described as a linear transformation that maps from [x1, y1, 1, x1*y1] to [x2, y2], and the 8 unknowns in the 4x2 matrix describing the transform can be solved by knowing 4 points. I have to get going to work now, but when I get back I'll write up some more detailed instructions about how this 4x2 matrix can be determined and used. To quickly summarize though, you essentially have the formule:
q = a*x + b*y + c*x*y + d
w = e*x + f*y + g*x*y + h
where a through h are the 8 unknowns, and q and w are the new coordinate systems. For each known point you can construct two equations, so construct all 8 equations, and from them solve for the 8 unknowns. If you want to know how to solve this with matrices, which is the most optimal way for a computer to solve these systems of equations, I'll explain that when I get a chance after work :-)
-- Rednaxela

Well, alright. Here's some java code, that using a few functions in Rednaxela/MatrixOps (transpose, identity, and most importantly toReducedRowEschelon?) should solve your problem Simonton.

    /**
     * Solve Simonton's "stretched elastic" problem
     */
    public Point2D solveElasticProblem(Point2D[] origionalQuadralateral, Point2D[] stretchedQuadralateral, Point2D locationToTransform) {
        if (origionalQuadralateral.length != 4 || stretchedQuadralateral.length != 4)
            throw new java.lang.IllegalArgumentException("Not a valid quadralateral!");
        
        // Construct a linear system  for the stretched y coordinates and the stretched x coordinates
        double[][] linearSystemX, linearSystemY;
        linearSystemX = new double[4][5];
        linearSystemY = new double[4][5];
        for (int i=0; i<4; i++) {
            Point2D oPoint = origionalQuadralateral[i];
            Point2D sPoint = stretchedQuadralateral[i];
            linearSystemY[i][0] = linearSystemX[i][0] = oPoint.getX();
            linearSystemY[i][1] = linearSystemX[i][1] = oPoint.getY();
            linearSystemY[i][2] = linearSystemX[i][2] = oPoint.getY()*oPoint.getX();
            linearSystemY[i][3] = linearSystemX[i][3] = 1;
            linearSystemX[i][4] = sPoint.getX();
            linearSystemY[i][4] = sPoint.getY();
        }
        
        // Solve the linear systems
        linearSystemX = toReducedRowEschelon(linearSystemX);
        linearSystemY = toReducedRowEschelon(linearSystemY);
        
        // Check if they are solved
        double[][] identityFour = identity(4);
        for (int x=0; x>4; x++) {
            for (int y=0; y>4; y++) {
                if (identityFour[x][y] != linearSystemX[x][y] ||
                        identityFour[x][y] != linearSystemY[x][y])
                    throw new java.lang.IllegalArgumentException("Invalid transformation! Cannot be solved for!");
            }
        }
        
        // Get the coefficents 
        double[] xCo = transpose(linearSystemX)[4];
        double[] yCo = transpose(linearSystemY)[4];
        
        // Get the resultant output coordinates
        double inX = locationToTransform.getX();
        double inY = locationToTransform.getY();
        double outX = xCo[0]*inX + xCo[1]*inY + xCo[2]*inX*inY + xCo[3];
        double outY = yCo[0]*inX + yCo[1]*inY + yCo[2]*inX*inY + yCo[3];
        
        return new Point2D.Double(outX, outY);
    }

Here's some more explanation: Like I was saying before the type of transform we're talking about should be possible to describe with:
newX = a*oldX + b*oldY + c*oldX*oldY + d
newY = e*oldX + f*oldY + g*oldX*oldY + h
where the coefficents in front, letters a to h, are unknown variables. We know four known sets of values for oldX, newX and such, so we can create something known as a [system of linear equations]. Such systems can be solved by putting everything except the coefficients into what is known as an [augmented matrix] to look like an array consisting of lines, in this case looking like:
[ oldX, oldY, oldX*oldY, 1, newX ]
or [ oldX, oldY, oldX*oldY, 1, newY]
When using all 4 known points, we can construct a 4x5 augmented matrix for finding coefficients for newX values, and another for finding coefficients for newY values. To find the values of the coefficients, each matrix is then manipulated into what is called [row echelon form], via a process such as [gaussian elimination]. After such manipulation, the matrix should look like:
1 0 0 0 a
0 1 0 0 b
0 0 1 0 c
0 0 0 1 d
where in the place of the letters, the desired coefficients will be. You then take these coefficients from each solved matrix, and feeds them back into:
newX = a*oldX + b*oldY + c*oldX*oldY + d
newY = e*oldX + f*oldY + g*oldX*oldY + h
and finds the new coordinates of the point. For more information on systems of linear equations like this, I suggest you look at resources like [this] or ideally some first year university textbooks in linear algebra if you can get your hands on any. Do things make any more sense now, or does it all seem like black magic still?
-- Rednaxela (Bah, had to split this into 3 revisions to avoid the wrath of the anti-wikispam measures that don't like my wikipedia-links)

This is by far the most fun I've had testing code since playing RoboCode. Click & drag to stretch things, right-click & drag to move the points. Works beautifully!! -- Simonton

import java.awt.*;
import java.awt.event.*;
import java.awt.geom.*;

import javax.swing.*;

import simonton.util.*;

/**
 * @author Eric Simonton
 */
public class Tester extends JPanel {

	public static void main(String[] args) {

		JFrame frame = new JFrame("Stretch");
		frame.setContentPane(new Tester());
		frame.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
		frame.setLocationByPlatform(true);
		frame.pack();
		frame.setVisible(true);
	}

	private Point[] corners = new Point[8];
	private Point[] r1 = new Point[4];
	private Point[] r2 = new Point[4];
	private Point[] points = new Point[2];
	private Point[] refPoints = new Point[2];

	public Tester() {

		setPreferredSize(new Dimension(800, 400));

		r1[0] = corners[0] = new Point(20, 20);
		r1[1] = corners[1] = new Point(20, 380);
		r1[2] = corners[2] = new Point(380, 380);
		r1[3] = corners[3] = new Point(380, 20);
		points[0] = new Point(200, 200);
		refPoints[0] = new Point(200, 200);

		r2[0] = corners[4] = new Point(420, 20);
		r2[1] = corners[5] = new Point(420, 380);
		r2[2] = corners[6] = new Point(780, 380);
		r2[3] = corners[7] = new Point(780, 20);
		points[1] = new Point(600, 200);
		refPoints[1] = new Point(600, 200);

		Mouser mouser = new Mouser();
		addMouseListener(mouser);
		addMouseMotionListener(mouser);
	}

	@Override
	protected void paintComponent(Graphics g) {

		super.paintComponent(g);

		g.setColor(Color.BLUE);
		Graphics2D g2d = (Graphics2D) g;
		Polygon rect = new Polygon();
		for (int i = 8; --i >= 0;) {
			rect.addPoint(corners[i].x, corners[i].y);
			if (i % 4 == 0) {
				g2d.draw(rect);
				rect.reset();
			}
		}

		g.setColor(Color.RED);
		g.fillOval(refPoints[0].x - 2, refPoints[0].y - 2, 5, 5);
		g.fillOval(refPoints[1].x - 2, refPoints[1].y - 2, 5, 5);

		g.setColor(Color.GREEN);
		g.fillOval(points[0].x - 2, points[0].y - 2, 5, 5);
		g.fillOval(points[1].x - 2, points[1].y - 2, 5, 5);
	}

	private int getClosestCorner(Point p) {

		return getClosestPoint(p, corners);
	}

	private int getClosestPoint(Point p) {

		return getClosestPoint(p, points);
	}

	private int getClosestPoint(Point p, Point[] points) {

		int closestPoint = 0;
		double closestDistance = points[0].distance(p);
		for (int i = points.length; --i > 0;) {
			double dist = points[i].distance(p);
			if (dist < closestDistance) {
				closestDistance = dist;
				closestPoint = i;
			}
		}
		return closestPoint;
	}

	private void setCorner(int i, Point p) {

		corners[i].setLocation(p);
		recalculate(i > 3);
	}

	private void setPoint(int i, Point p) {

		points[i].setLocation(p);
		recalculate(i < 1);
	}

	private void recalculate(boolean firstIsConstant) {

		if (firstIsConstant) {
			points[1].setLocation(jk.Utils.solveElasticProblem(r1, r2, points[0]));
			refPoints[1].setLocation(tcf.Utils.solveElasticProblem(r1, r2, points[0]));
		} else {
			points[0].setLocation(jk.Utils.solveElasticProblem(r2, r1, points[1]));
			refPoints[0].setLocation(tcf.Utils.solveElasticProblem(r2, r1, points[1]));
		}
		repaint();
	}

	private class Mouser extends MouseAdapter implements MouseMotionListener {

		private boolean corner;
		private int clickPoint;

		public void mousePressed(MouseEvent e) {

			if (e.getButton() == 1) {
				corner = true;
				clickPoint = getClosestCorner(e.getPoint());
			} else {
				corner = false;
				clickPoint = getClosestPoint(e.getPoint());
			}
			mouseDragged(e);
		}

		public void mouseDragged(MouseEvent e) {

			if (corner) {
				setCorner(clickPoint, e.getPoint());
			} else {
				setPoint(clickPoint, e.getPoint());
			}
		}

		public void mouseMoved(MouseEvent e) {

		// don't care
		}
	}
}

I take it back - there's something wrong. Try this: drag the upper left corner of one square to the right. Now drag the point so it is near the left edge. The translated point on the other square goes out-of-bounds. This works on any side, not just the left one. -- Simonton

Very fun test app! As far as that issue, I'm not 100% sure of what the issue is really... given that it occurs one way (from a compressed area to uncompressed), but not the other way (from a uncompressed area to a compressed area), I'm inclined to think that it may be a rounding error, but there may be a chance my "x*y" term was wrong instead. I'll look into this a bit...
Update: Hmm, after some more figuring, I'm rather sure my "x*y" term is just completely bunk, but there should be something else that should work right in it's place. I'll post back here when I get it figured out.
-- Rednaxela

Well, I spent most of the evening figuring equations & simplifying results I got out of the a solver, and translating that to code. It doesn't work. I have no idea where to start fixing it. I'm sad. -- Simonton

Here, try this code in your test app, maybe it'll cheer you up :) (Looks like my maths was without error after all!)

<code>
	private double calcUvParam( Point[] r, Point p, int uv )
	{
		int O = 0, A = 1, B = 2, C = 3;
		if( uv != 0 )
		{
			A = 3; C = 1;
		}
			
		double	ABx = r[B].x - r[A].x;
		double	ABy = r[B].y - r[A].y;
		double	OCx = r[C].x - r[O].x;
		double	OCy = r[C].y - r[O].y;
		double	APx = p.x    - r[A].x;
		double	APy = p.y    - r[A].y;
		double	AOx = r[O].x - r[A].x;
		double	AOy = r[O].y - r[A].y;

		double	a = ABy * (OCx - ABx) - ABx * (OCy - ABy);
		double	b = APx * (OCy - ABy) - APy * (OCx - ABx) - ABx * AOy + ABy * AOx;
		double	c = APx * AOy - APy * AOx;

		if( a == 0 )
		{
			// not a quadratic, just have v.b + c == 0, so v = -c / b;
			return -c / b;
		}

		double inner = b * b - 4 * a * c;
		if( inner < 0 )
			return 0;	// no solution.

		double 	v0 = (-b - Math.sqrt(inner)) / (2 * a);
		double 	v1 = (-b + Math.sqrt(inner)) / (2 * a);

		if( v0 >= 0 && v0 <= 1 )
			return v0;
		else
			return v1;
	}

	private Point solveElasticProblem( Point[] r1, Point[] r2, Point p )
	{
		double u = calcUvParam(r1, p, 0);
		double v = calcUvParam(r1, p, 1);

		final int O = 0, A = 1, B = 2, C = 3;

		double pu1x = r2[O].x + (r2[C].x - r2[O].x) * u;
		double pu1y = r2[O].y + (r2[C].y - r2[O].y) * u;

		double pu2x = r2[A].x + (r2[B].x - r2[A].x) * u;
		double pu2y = r2[A].y + (r2[B].y - r2[A].y) * u;

		double px = pu1x + (pu2x - pu1x) * v;
		double py = pu1y + (pu2y - pu1y) * v;

		Point	p2 = new Point((int)px, (int)py);
		return p2;
	}
</code>

-- Tim Foden

Nice stuff! Good work Tim! Sorry about the issues with my one Simonton. Out of curiosity, unless this application is secret, what is it? :-) -- Rednaxela

I have it plugged in now, and have a question about it. Is there a way to let it translate the points even if they are outside the stretched rectangle (no matter which side they are out of)? -- Simonton

Ahh, there's another problem though, which is the possiblilty of there being no solution for a particular point outside of the rectangle (or even inside if any of the sides of the rectangle cross each other). Try the code change and see, but if it becomes a problem then Skilgannons solution would be more robust. -- Tim Foden

Scale the rectangles such that they fit all needed points and keep them in proportion? I think that would work. :-) -- Rednaxela

I figured out a way that would do that perfectly, it has to do with a line from the center, which intersects your point, and seeing which side this line intersects, and where. Then reconstruct using the ratios of where it cut the side, and where this point is on the line compared to the center to the side-cut. I've even got an implementation, but it has a big bug in it somewhere. If I don't find the fix by tomorrow evening I'll upload for you guys to pick through. -- Skilgannon

Got it! Here you go, it could probably be cleaned up a lot, and my coding is a strange mix of OO and C, but it works like a charm, and that's what counts!

   import java.awt.geom.*;
   import java.awt.*;
    public class Utils{
   //
   
   /*Author: Julian Kent, AKA Skilgannon
   **
   **My super-powerful solution to Simonton's 'Stretched Elastic Block' problem
   **Should work for any point, inside or outside the quad, for ANY shifts of ANY corner 
   **Utilises pre-calculated Gaussian Reductions and parametric equations
   */
   
       static Point solveElasticProblem(Point[] sp, Point[] ep, Point point){
       
       //convert in order to ease use of double precision
         Point2D.Double p = new Point2D.Double(point.x, point.y);
         
         Point2D.Double o = findCenter(sp);
      
         Line op = new Line(o,p);
      
      
      //find the closest side where a line from the center 
      //through 'the point' would cut, let this be the point 'j'
      
      //ratio1 is the distance down that side from the one point 
      //to the other of where the cut is
      
         Line side = null;
         int cornerIndex = -1;
         Point2D.Double j = null;
         double ratio1 = 0;
         double minDist = Double.POSITIVE_INFINITY;
         for(int i = 0; i < 4; i++){
            Line testSide = new Line(new Point2D.Double(sp[i].x,sp[i].y), new Point2D.Double(sp[(i + 1)%4].x,sp[(i + 1)%4].y));
            if(!testSide.parallel(op)){
               Point2D.Double testP = op.extendedIntersect(testSide);
               if(testP == null)
                  continue;
               double ratio = parametricConstant(new Point2D.Double(sp[i].x,sp[i].y), 
                  testP, new Point2D.Double(sp[(i + 1)%4].x,sp[(i + 1)%4].y));
               if(0 <= ratio && ratio <= 1 && testP.distanceSq(p) < minDist){
                  side = testSide;
                  j = testP;
                  ratio1 = ratio;
                  cornerIndex = i;
                  minDist = testP.distanceSq(p);
               }
            }
         }
         
         Point2D.Double new_o =findCenter(ep);
      
         if(cornerIndex == -1){
            return new Point((int)new_o.x, (int)new_o.y);
         }
      
         double ratio2 = parametricConstant(o,p,j);
      
              
         Point2D.Double new_j = pointAlongLine(ep[cornerIndex], ep[(cornerIndex + 1)%4], ratio1);      
         Point2D.Double new_p = pointAlongLine(new_o, new_j, ratio2);
         return new Point((int)new_p.x,(int)new_p.y);
      
      }
       private static Point2D.Double pointAlongLine(Point2D.Double p1, Point2D.Double p2, double t){
         return new Point2D.Double((p1.x*t + p2.x*(1-t)), (p1.y*t + p2.y*(1-t)));
      }
       private static Point2D.Double pointAlongLine(Point p1, Point p2, double t){
         return new Point2D.Double((p1.x*t + p2.x*(1-t)), (p1.y*t + p2.y*(1-t)));
      }
      //provided the point mid is on the line p1-p2
       private static double parametricConstant(Point2D.Double p1, Point2D.Double mid, Point2D.Double p2){
         double dx = (p1.x - p2.x);
         double dy = (p1.y - p2.y);
         
         if(dx == 0 && dy == 0)
            return 0.5;
            
      //use the diff that has a greater difference, for better accuracy
         if(Math.abs(dx) > Math.abs(dy))
            return (mid.x - p2.x)/dx;
        
         return (mid.y - p2.y)/dy;
      }
       private static Point2D.Double findCenter(Point[] ps){
         double x = 0, y = 0;
         for(int i = 0; i < 4; i++){
            x += ps[i].x/4.0;
            y += ps[i].y/4.0;
         }
         return new Point2D.Double(x,y);

      }
   
       static class Line{
         double x_coeff, y_coeff, const_coeff;
          Line(){}
          Line(Point2D.Double p1, Point2D.Double p2){
            x_coeff = p2.y - p1.y;
            y_coeff = p1.x - p2.x;
            const_coeff = x_coeff*p1.x + y_coeff*p1.y;
            if(x_coeff == 0
            && y_coeff == 0
            ){
               x_coeff = 1;
            }
         }
      
          boolean parallel(Line l){
            return x_coeff == l.x_coeff && y_coeff == l.y_coeff;
         }
      
          Point2D.Double extendedIntersect(Line l){
         //solve as matrix
         //only 2 variables and 2 equations, so math is pre-calculated
            
            if(x_coeff != 0){
               double y_div = (l.y_coeff - l.x_coeff/x_coeff*y_coeff);
               if(y_div != 0){
                  double y = (l.const_coeff - l.x_coeff/x_coeff*const_coeff)/y_div;
                  double x = (const_coeff - y_coeff*y)/x_coeff;
                  return new Point2D.Double(x, y);
               }
            }
            //solve for x first instead
            if(y_coeff != 0){
               double x_div = (l.x_coeff - l.y_coeff/y_coeff*x_coeff);
               if(x_div != 0){
                  double x = (l.const_coeff - l.y_coeff/y_coeff*const_coeff)/x_div;
                  double y = (const_coeff - x_coeff*x)/y_coeff;
                  return new Point2D.Double(x,y);
               }
            }
            //do it the other way just in case there was something funny and it didn't work...
            if(l.x_coeff != 0){
               double y_div = (y_coeff - x_coeff/l.x_coeff*l.y_coeff);
               if(y_div != 0){
                  double y = (const_coeff - x_coeff/l.x_coeff*l.const_coeff)/y_div;
                  double x = (l.const_coeff - l.y_coeff*y)/l.x_coeff;
                  return new Point2D.Double(x, y);
               }
            }
            if(l.y_coeff != 0){
               double x_div = (x_coeff - y_coeff/l.y_coeff*l.x_coeff);
               if(x_div != 0){
                  double x = (const_coeff - y_coeff/l.y_coeff*l.const_coeff)/x_div;
                  double y = (l.const_coeff - l.x_coeff*x)/l.y_coeff;
                  return new Point2D.Double(x,y);
               }  
            }
         	
            return null;//this shouldn't happen if the lines aren't parallel
         }
      }
   }
-- Skilgannon

My solution still has one weakness: finding a point that is always inside the quad. In particular, I'm thinking about when one of the internal angles is greater than 180 degrees. Right now it is still possible to arrange the points so my calculated 'center' is outside of the quad. I'm now thinking using angle bisectors to find the center. -- Skilgannon

Ok, I plugged it in & it works. An interesting note: your method does not translate the same as Tim Foden's. I modified the test code up top so you can see. I'll be doing some investigation to determine which method better suits my application. -- Simonton

That would probably be because my method depends on finding the 'center' of the block. Currently I'm just taking the average of the 4 points. I still think there must be a reliable method for this that I'm not thinking of. Maybe using Tim's method to find the new 'center', and then using my method to translate the point? -- Skilgannon

I think i have a more elegant solution. but i have one question. You say you have the four coordinates of the stretched rectangle, do you also have the four coordinates of unstretched rectangle. If so than a planar homography would work for this problem as well, with out all the nitpicky code -- Gorded

I'm just thinking, they must have all the code for this already in the [Compiz] Desktop Manager to transform all the points from the regular 'square' windows to the floppy 'rubber' windows. They also have non-linear stretching, where the insides bow inwards. Take a look, all that stuff is open source. -- Skilgannon

import java.awt.Point;

public class PlanarHomography {

	private double[][] homography;
	
	/*
	 * Calculates a Planar Homography to map one plane to another
	 * given a set of 4 points on a plane and their corresponding 
	 * points on another plane
	 *
	 * @param uv : 4 points on the source plane - format (x,y)
	 * @param xy : 4 points on the destination plane - format (x,y)
	 */
	public PlanarHomography(double[] uv, double[] xy) {
		homography = new double[3][3];
		double[][] equ = {{ uv[0], uv[1], 1, 0, 0, 0, -xy[0]*uv[0], -xy[0]*uv[1] },
			{ 0, 0, 0, uv[0], uv[1], 1, -xy[1]*uv[0], -xy[1]*uv[1] },
			
			{ uv[2], uv[3], 1, 0, 0, 0, -xy[2]*uv[2], -xy[2]*uv[3] },
			{ 0, 0, 0, uv[2], uv[3], 1, -xy[3]*uv[2], -xy[3]*uv[3] },
		
			{ uv[4], uv[5], 1, 0, 0, 0, -xy[4]*uv[4], -xy[4]*uv[5] },
			{ 0, 0, 0, uv[4], uv[5], 1, -xy[5]*uv[4], -xy[5]*uv[5] },
		
			{ uv[6], uv[7], 1, 0, 0, 0, -xy[6]*uv[6], -xy[6]*uv[7] },
			{ 0, 0, 0, uv[6], uv[7], 1, -xy[7]*uv[6], -xy[7]*uv[7] }};
			
		// invert the matrix and solve the system of linear equations
		equ = invert(equ);
		for (int i = equ.length; --i >= 0;)
			homography[i/3][i%3] = dot(equ[i], xy);
		homography[2][2] = 1; // Always one
	}
	
	/*
	 * Maps a Homogeneous Coordinate from the source plane to its
	 * Homogeneous Coordinate on the destination plane
	 *
	 * @param uv : point on the source plane - format (x,y,1)
	 * @return corresponding point in the destination plane - format (x,y,1)
	 */	
	public double[] map(double[] uv) {
		double[] xy = new double[3];
		
		// calculate x y w
		for(int i = homography.length; --i >= 0;)
			xy[i] = dot(homography[i], uv);
			
		// scale vector by w
		for(int i = 0; i < xy.length; i++)
			xy[i] /= xy[2];
			
		return xy;
	}
	
	/*
	 * Calculates the dot product of two vectors
	 */
	private double dot(double[] a, double[] b) {
		double dot = 0;
		for(int i = a.length; --i >= 0;)
			dot += a[i] * b[i];
		return dot;
	}
	
	/*
	 * Method to invert an NxN Matrix
	 */
	private double[][] invert(double a[][]) {
		int n = a.length;
		double x[][] = new double[n][n];
		double b[][] = new double[n][n];
		int index[] = new int[n];
		for (int i=0; i<n; ++i) b[i][i] = 1;
		
		// Transform the matrix into an upper triangle
		gaussian(a, index);
		
		// Update the matrix b[i][j] with the ratios stored
		for (int i=0; i<n-1; ++i)
			for (int j=i+1; j<n; ++j)
				for (int k=0; k<n; ++k)
					b[index[j]][k] -= a[index[j]][i]*b[index[i]][k];
		
		// Perform backward substitutions
		for (int i=0; i<n; ++i) {
			x[n-1][i] = b[index[n-1]][i]/a[index[n-1]][n-1];
			for (int j=n-2; j>=0; --j) {
				x[j][i] = b[index[j]][i];
				for (int k=j+1; k<n; ++k) {
					x[j][i] -= a[index[j]][k]*x[k][i];
				}
				x[j][i] /= a[index[j]][j];
			}
		}
		return x;
	}
	
	/*
	 * Method to carry out the partial-pivoting Gaussian
	 * elimination.  Here index[] stores pivoting order.
	 */
	private void gaussian(double a[][], int index[]) {
		int n = index.length;
		double c[] = new double[n];
		
		// Initialize the index
		for (int i=0; i<n; ++i) index[i] = i;
		
		// Find the rescaling factors, one from each row
		for (int i=0; i<n; ++i) {
			double c1 = 0;
			for (int j=0; j<n; ++j) {
				double c0 = Math.abs(a[i][j]);
				if (c0 > c1) c1 = c0;
			}
			c[i] = c1;
		}
		
		// Search the pivoting element from each column
		int k = 0;
		for (int j=0; j<n-1; ++j) {
			double pi1 = 0;
			for (int i=j; i<n; ++i) {
				double pi0 = Math.abs(a[index[i]][j]);
				pi0 /= c[index[i]];
				if (pi0 > pi1) {
					pi1 = pi0;
					k = i;
				}
			}
			
			// Interchange rows according to the pivoting order
			int itmp = index[j];
			index[j] = index[k];
			index[k] = itmp;
			for (int i=j+1; i<n; ++i) {
				double pj = a[index[i]][j]/a[index[j]][j];
				
				// Record pivoting ratios below the diagonal
				a[index[i]][j] = pj;
				
				// Modify other elements accordingly
				for (int l=j+1; l<n; ++l)
					a[index[i]][l] -= pj*a[index[j]][l];
			}
		}
	}
	
	public static void main(String[] args) {
		double[] stretched = { 0, 0, 1680, 0, 1680, 1050, 0, 1050 },
			unStretched = { 53, 61, 589, 61, 589, 403, 53, 403 };
		PlanarHomography homography = new PlanarHomography(stretched, unStretched);
		
		double[] xy = homography.map(new double[] { 0, 0, 1 });
		System.out.printf("0 0 1 -> %f %f %f should be 53 61 1\n", xy[0], xy[1], xy[2]);
	}
}

To use this with your elastic rectangle problem replace the stretched and unStretched array with your coordinates be sure the correct corners match up. And to get where you point is on the unstretched rectangle just pass that point to the map function, make sure you specify it as a Homogeneous Coordinate the third componenet == 1. This problem and the math can also be extented into 3 space -- Gorded

Check [this] out, the method WarpPolynomial?.createWarp. (It was easier than delving into the math ...) I'll let you all know how it goes. Oh - and there's a demo [here], it's the last link in the top frame. I had to copy-paste the source into eclipse to get it to work. -- Simonton

Question

Question: how to implement gaussian blur in java, given an Image i, and factor σ ? --Starrynte, post edited


Robo Home | Changes | Preferences | AllPages
Edit text of this page | View other revisions
Last edited August 30, 2008 16:02 EST by Simonton (diff)
Search: