1
\$\begingroup\$

I want to calculate the transform I need to render a 3D object in my world.

How I see my idea:

  • You open a camera
  • You point it at a QR code
  • The app detects the QR code and returns the 4 corners and the correct direction
  • The app calculates the transform that this QR code is currently at (I need help with this step)
  • The app uses that transform to render the 3D object in a fixed position over the QR code

I want something like Vuforia, but for QR codes (as they are used often).

I don't really know how I can extract the transform from the detected QR code (I will get the location of the 4 corners)

\$\endgroup\$
3
\$\begingroup\$

Deriving the transformation

Note that the transformation of the QR code captured by the camera is a perspective transformation. As such, no amount of affine transformation can emulate it (translating, scaling, rotating, shearing, and reflecting are not enough).

From the surface where the QR code is to the image, a perspective transformation is applied. Your task is to solve a projection matrix for that transformation.


Ok, we are looking for a matrix such that

image_point = M * world_point

This matrix needs to be a 3 by 3 matrix to be able to express the projection we need. Which also means we will augment the vectors with an additional value:

+-           -+   +-       -+   +-       -+
| image_x * w |   | a  b  c |   | world_x |
| image_y * w | = | d  e  f | * | world_y |
|       w     |   | g  h  1 |   |    1    |
+-           -+   +-       -+   +-       -+

Observe that:

image_x * w = a * x + b * y + c
image_y * w = d * x + e * y + f
w = gx + hy + 1

Therefore

image_x = (a * x + b * y + c) / (g * x + h * y + 1)
image_y = (d * x + e * y + f) / (g * x + h * y + 1)

Remember that the above formulas are the transformation of a point in the real world to a point in the image.

We also know the positions where they ended up in the camera... now you have to solve the equation system to find out the values of the constants a, b, c, d, e, f, g, and h that form the projection matrix.

For the positions in the real world, we will say that the first point is in position (0, 0) in the world, and that the next is at (S, 0), the other at (0, S) and the fourth point will be at (S, S). Where S is the desired size in the final image.

Those are seven incognitas, for eight equations (two for each one of the four points). This is doable.

Once you have the projection matrix, you can:

  • apply the inverse of the transformation to the image, and that will result in an image where the square is a square (within measurement error).

  • apply it to 3D models to match the perspective of the camera.


Ok, four points:

image_x1 = (a * x1 + b * y1 + c) / (g * x1 + h * y1 + 1)
image_y1 = (d * x1 + e * y1 + f) / (g * x1 + h * y1 + 1)

image_x2 = (a * x2 + b * y2 + c) / (g * x2 + h * y2 + 1)
image_y2 = (d * x2 + e * y2 + f) / (g * x2 + h * y2 + 1)

image_x3 = (a * x3 + b * y3 + c) / (g * x3 + h * y3 + 1)
image_y3 = (d * x3 + e * y3 + f) / (g * x3 + h * y3 + 1)

image_x4 = (a * x4 + b * y4 + c) / (g * x4 + h * y4 + 1)
image_y4 = (d * x4 + e * y4 + f) / (g * x4 + h * y4 + 1)

Replace with our invented world coordinates:

image_x1 = (a * 0 + b * 0 + c) / (g * 0 + h * 0 + 1)
image_y1 = (d * 0 + e * 0 + f) / (g * 0 + h * 0 + 1)

image_x2 = (a * S + b * 0 + c) / (g * S + h * 0 + 1)
image_y2 = (d * S + e * 0 + f) / (g * S + h * 0 + 1)

image_x3 = (a * 0 + b * S + c) / (g * 0 + h * S + 1)
image_y3 = (d * 0 + e * S + f) / (g * 0 + h * S + 1)

image_x4 = (a * S + b * S + c) / (g * S + h * S + 1)
image_y4 = (d * S + e * S + f) / (g * S + h * S + 1)

Remove the zero terms:

image_x1 = c
image_y1 = f

image_x2 = (a * S + c) / (g * S + 1)
image_y2 = (d * S + f) / (g * S + 1)

image_x3 = (b * S + c) / (h * S + 1)
image_y3 = (e * S + f) / (h * S + 1)

image_x4 = (a * S + b * S + c) / (g * S + h * S + 1)
image_y4 = (d * S + e * S + f) / (g * S + h * S + 1)

Look, we have c and f right off the bat (remember that you have all the coordinates in the image). That's the equation system to solve.

I will assume S = 1, to keep it simple. Will come back to this point.

image_x1 = c
image_y1 = f

image_x2 = (a + image_x1) / (g + 1)
image_y2 = (d + image_y1) / (g + 1)

image_x3 = (b + image_x1) / (h + 1)
image_y3 = (e + image_y1) / (h + 1)

image_x4 = (a + b + image_x1) / (g + h + 1)
image_y4 = (d + e + image_y1) / (g + h + 1)

Now, solve for a, d, b and e:

image_x2 * (g + 1) - image_x1 = a
image_y2 * (g + 1) - image_y1 = d

image_x3 * (h + 1) - image_x1 = b
image_y3 * (h + 1) - image_y1 = e

Replace:

image_x4 = ((image_x2 * (g + 1) - image_x1) + (image_x3 * (h + 1) - image_x1) + image_x1) / (g + h + 1)
image_y4 = ((image_y2 * (g + 1) - image_y1) + (image_y3 * (h + 1) - image_y1) + image_y1) / (g + h + 1)

Solve for g:

image_x4 = ((image_x2 * (g + 1) - image_x1) + (image_x3 * (h + 1) - image_x1) + image_x1) / (g + h + 1)

=>

g = (-image_x1 -(h + 1) * image_x4 + h * image_x3 + image_x3 + image_x2) / (image_x4 - image_x2)

Replace in the final equation:

image_y4 = ((image_y2 * (((-image_x1 -(h + 1) * image_x4 + h * image_x3 + image_x3 + image_x2) / (image_x4 - image_x2)) + 1) - image_y1) + (image_y3 * (h + 1) - image_y1) + image_y1) / (((-image_x1 -(h + 1) * image_x4 + h * image_x3 + image_x3 + image_x2) / (image_x4 - image_x2)) + h + 1)

Solve for h:

h = (image_x1 * (image_y2 - image_y4) + image_y1 * (image_x4 - image_x2) + image_x2 * image_y3 - image_y2 * image_x3 + image_x3 * image_y4 - image_y3 * image_x4)/(image_x2 * (image_y4 - image_y3) + image_y2 * (image_x3 - image_x4) - image_x3 * image_y4 + image_y3 * image_x4)

Finally, we have h in terms of positions in the image!

Go ahead, and replace it on g... having g and h you should be able to solve a, b, d, e. I am not posting any of that here because the expressions are too large.

Suffice to say that you may compute h, and plug it on the formula for g. Then plug them in the formulas for a, b, d, e that I used above. With that, you will have all the terms for the perspective matrix.


Since we got rid of S to keep it simple...

  • If you want to transform from world coordinates to image coordinates first scale by 1/S then apply the perspective transformation. This is the operation used to place an object in the image.

  • If you want to transform from image coordinates to world coordinates, apply the inverse transformation and then scale by S. This would be operation used, to get the image captured by the camera and transform it to true shape.

  • If you want to place 3D models, scale them to the factor S you wanted before applying the perspective transformation.

    Notes:

    • The origin is whichever you chose as first point, since we assumed it is at (0, 0) on the captured surface.

    • The transformation matrix we found is for 2D. Without knowledge of the real size of the QR (or some other reference information), we cannot really tell how far it is. You may have to add a scaling factor for the third dimension.

    • Similarly, the matrix has no concept of what side is up or down. You may have to add a transformation that flips the third dimension.

I have not tested placing 3D objects


Demonstration

When I implemented this, I discovered that there is a problem when image_x4 is equal to image_x2. The problem is that you get division by zero when computing g... I have tweaked this by replacing zero by epsilon.

Note: the inverse matrix for 2D non-affine transformation similar to the one presented above is the Adjugate matrix:

+- -+   +-     -+   +-                    -+
| x |   | a b c |   | x′*w′ = ax + by + c  |
| y | * | d e f | = | y′*w′ = dx + ey + f  |
| 1 |   | g h i |   |   w′  = gx + hy + i  |
+- -+   +-     -+   +-                    -+

=>

+-  -+   +-                 -+   +-   -+
| x′ |   | ei-fh ch-bi bf-ce |   | x*w |
| y′ | * | fg-di ai-cg cd-af | = | y*w |
| 1  |   | dh-eg bg-ah ae-bd |   |  w  |
+-  -+   +-                 -+   +-   -+

The code below does not include a method to compute an "inverse matrix" but to apply the inverse transformation of a given matrix.

/**
* @export
* */
var Matrix = (function () {
    /**
    * @constructor
    * @param {number} a
    * @param {number} b
    * @param {number} c
    * @param {number} d
    * @param {number} e
    * @param {number} f
    * @param {number} g
    * @param {number} h
    * @param {number} i
    */
    function Matrix(a, b, c, d, e, f, g, h, i) {
        if (typeof g === "undefined") { g = 0; }
        if (typeof h === "undefined") { h = 0; }
        if (typeof i === "undefined") { i = 1; }
        this.a = a;
        this.b = b;
        this.c = c;
        this.d = d;
        this.e = e;
        this.f = f;
        this.g = g;
        this.h = h;
        this.i = i;
    }
    /**
    * @param {number} x
    * @param {number} y
    * @return {Matrix}
    * @export
    */
    Matrix.translation = function (x, y) {
        return new Matrix(1, 0, x, 0, 1, y, 0, 0, 1);
    };

    /**
    * @param {number} theta
    * @return {Matrix}
    * @export
    */
    Matrix.rotation = function (theta) {
        return new Matrix(Math.cos(theta), -Math.sin(theta), 0, Math.sin(theta), Math.cos(theta), 0, 0, 0, 1);
    };

    /**
    * @param {number} scale
    * @return {Matrix}
    * @export
    */
    Matrix.scale = function (scale) {
        return new Matrix(scale, 0    , 0,
                          0    , scale, 0,
                          0    , 0    , 1);
    };

    /**
    * @param {Matrix} a
    * @param {Matrix} b
    * @return {Matrix}
    * @export
    */
    Matrix.product = function (a, b) {
        return new Matrix((a.a * b.a) + (a.b * b.d) + (a.c * b.g), (a.a * b.b) + (a.b * b.e) + (a.c * b.h), (a.a * b.c) + (a.b * b.f) + (a.c * b.i), (a.d * b.a) + (a.e * b.d) + (a.f * b.g), (a.d * b.b) + (a.e * b.e) + (a.f * b.h), (a.d * b.c) + (a.e * b.f) + (a.f * b.i), (a.g * b.a) + (a.h * b.d) + (a.i * b.g), (a.g * b.b) + (a.h * b.e) + (a.i * b.h), (a.g * b.c) + (a.h * b.f) + (a.i * b.i));
    };

    /**
    * @param {Matrix} m
    * @param v
    * @export
    */
    Matrix.transform = function (m, v) {
        var w = v.x * m.g + v.y * m.h + m.i;
        return {x: (v.x * m.a + v.y * m.b + m.c) / w, y:(v.x * m.d + v.y * m.e + m.f) / w};
    };

   /**
    * @param {Matrix} m
    * @param v
    * @export
    */
    Matrix.transformInverse = function (m, v) {
        var w = (v.x*(m.d*m.h-m.e*m.g)+v.y*(m.b*m.g-m.a*m.h)+m.a*m.e-m.b*m.d);
        var y = (v.x*(m.f*m.g-m.d*m.i)+v.y*(m.a*m.i-m.c*m.g)+m.c*m.d-m.a*m.f)/w;
        var x = (v.x*(m.e*m.i-m.f*m.h)+v.y*(m.c*m.h-m.b*m.i)+m.b*m.f-m.c*m.e)/w;
        return {x:x, y:y};
    };

    /**
    * @param a
    * @param b
    * @param c
    * @param d
    * @return {Matrix}
    * @export
    */
    Matrix.perspective = function (a, b, c, d) {
        var _hd = (b.x * (d.y - c.y) + b.y * (c.x - d.x) - c.x * d.y + c.y * d.x);
        if (_hd === 0) _hd = Number.EPSILON;
        var _h = (a.x * (b.y - d.y) + a.y * (d.x - b.x) + b.x * c.y - b.y * c.x + c.x * d.y - c.y * d.x) / _hd;
        var _gd = (d.x - b.x);
        if (_gd === 0) _hd = Number.EPSILON;
        var _g = (-a.x - (_h + 1) * d.x + _h * c.x + c.x + b.x) / _gd;
        var _a = b.x * (_g + 1) - a.x;
        var _d = b.y * (_g + 1) - a.y;
        var _b = c.x * (_h + 1) - a.x;
        var _e = c.y * (_h + 1) - a.y;
        var _c = a.x;
        var _f = a.y;
        var _i = 1;
        return new Matrix(_a, _b, _c, _d, _e, _f, _g, _h, _i);
    };
    return Matrix;
})();

var time = (new Date()).getTime(); // time in milliseconds
var polygon =[{x:60,y:20},{x:100,y:40},{x:100,y:80},{x:60,y:100},{x:20,y:80},{x:20,y:40}];
var angularSpeed = Math.PI / 2;
var source;
var target;
var a = null; var b = null; var c = null; var d = null;
function update()
{
  let new_time = (new Date()).getTime();
  let elapsed = (new_time - time) / 1000.0;
  let angle = angularSpeed * elapsed;
  let translation = Matrix.translation(-60, -60);
  let rotation = Matrix.rotation(angle);
  let matrix = Matrix.product(rotation, translation);
  let points = '';
  for (let index = polygon.length - 1; index >= 0; index--)
  {
      polygon[index] = Matrix.transformInverse(translation, Matrix.transform(matrix, polygon[index]));
      points += polygon[index].x + ", " + polygon[index].y + " ";
  }
  source.setAttribute("points", points);
  if (a !== null && b !== null && c !== null && d !== null)
  {
    points = '';
    let perspective = Matrix.perspective(a, b, d, c);
    for (let index = polygon.length - 1; index >= 0; index--)
    {
        var tmp = Matrix.transform(perspective, Matrix.transform(Matrix.scale(1/120), polygon[index]));
        points += tmp.x + ", " + tmp.y + " ";
    }
    target.setAttribute("points", points);
  }
  time = new_time;
}

// --------------
// Ugly UI code
// --------------
var world;
function getNode(n, v){
  n = document.createElementNS("http://www.w3.org/2000/svg", n);
  for (var p in v) n.setAttributeNS(null, p, v[p]);
  return n;
}
jQuery(document).ready(function(){
  world = document.getElementById("world");
  source = document.getElementById("source");
  target = document.getElementById("target");
  var arrow = null;
  var ab = null; var bc = null; var cd = null; var da = null;
  $(world).on('mousemove', function(e){
    let x = e.pageX * 1000 / $(world).width(); let y = e.pageY * 1000 / $(world).height();
    if (d === null)
    {
        if (c === null)
        {
            if(b === null)
            {
                if (a === null)
                {
                    /**/
                }
                else
                {
                    arrow.setAttribute("x1", a.x); arrow.setAttribute("y1", a.y);
                    arrow.setAttribute("x2", x); arrow.setAttribute("y2", y);
                }
            }
            else
            {
                arrow.setAttribute("x1", b.x); arrow.setAttribute("y1", b.y);
                arrow.setAttribute("x2", x); arrow.setAttribute("y2", y);
            }
        }
        else
        {
            arrow.setAttribute("x1", c.x); arrow.setAttribute("y1", c.y);
            arrow.setAttribute("x2", x); arrow.setAttribute("y2", y);
        }
    }
  });
  jQuery(world).click(function(e){
    let x = e.pageX * 1000 / $(world).width(); let y = e.pageY * 1000 / $(world).height();
    if (d === null)
    {
        if (c === null)
        {
            if(b === null)
            {
                if (a === null)
                {
                    a = {x: x, y:y};
                    arrow = getNode('line', {x1: x, x2: x, y1: y, y2: y, "stroke-width":1, stroke:"black"});
                    world.appendChild(arrow);
                }
                else
                {
                    b = {x: x, y:y};
                    ab = getNode('line', {x1: a.x, x2: b.x, y1: a.y, y2: b.y, "stroke-width":1, stroke:"red"});
                    world.appendChild(ab);
                }
            }
            else
            {
                c = {x: x, y:y};
                bc = getNode('line', {x1: b.x, x2: c.x, y1: b.y, y2: c.y, "stroke-width":1, stroke:"red"});
                world.appendChild(bc);
            }
        }
        else
        {
            d = {x: x, y:y};
            cd = getNode('line', {x1: c.x, x2: d.x, y1: c.y, y2: d.y, "stroke-width":1, stroke:"red"});
            world.appendChild(cd);
            da = getNode('line', {x1: d.x, x2: a.x, y1: d.y, y2: a.y, "stroke-width":1, stroke:"red"});
            world.appendChild(da);
            world.removeChild(arrow);
            arrow = null;
        }
    }
    else
    {
        a = null;
        b = null;
        c = null;
        d = null;
        world.removeChild(ab);
        world.removeChild(bc);
        world.removeChild(cd);
        world.removeChild(da);
        ab = null;
        bc = null;
        cd = null;
        da = null;
    }
  });
  setInterval(update, 40);
});
body{margin:0}
<script src="https://ajax.googleapis.com/ajax/libs/jquery/2.1.1/jquery.min.js"></script>
<svg id="world" viewBox="0 0 1000 1000">
<polygon id="target" points="60,20 100,40 100,80 60,100 20,80 20,40" fill="aqua" />
<rect id="ref", x="0" y="0" width="120" height="120" fill="white" stroke="black" stroke-width="1"></rect>
<polygon id="source" points="" fill="red"/>
</svg>

To use the demo, click in the background to draw a quadrilateral by clicking in the four corners. The code will transform the rotating hexagon to the quadrilateral you draw. Clicking again will remove the quadrilateral and allow you to draw a new one.

Note: it assumes that whatever you draw is a convex cuadrilateral, no validation is made on that.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.