Monday, February 5, 2018

Perspective Transform




Among the various geometric transformations, the Perspective transform is very useful. By default, this transform is not available in Image > Transform. However, it is very easy to implement such a tool...


UPDATE: Perspective Transform was rewritten using the latest version of JavaScript (ECMAScript 2015+).

1. Input image


Fig.1: Distorted image

2. How to unwarp the image

  1. Download the script es6_perspective_transform.js at the end of this post and open it in ImageJ with File > Open...
  2. Then, select four points with the Multi-Points Tool in a zig-zag order as shown in Fig.2.
    3. Fig.2: Selection of four points with the Multi-Points Tool
  3. Finally, run the script, set the size of the output squared image and  the interpolation mode (Fig. 3).
    4. Fig.3: Dialog


  4. Et voilà! Depending of the interpolation mode (None or Bilinear), the rendering may be different.
Fig. 4: Unwarped images with interpolation modes: None (left) and Bilinear (Right). The image calculated with a bilinear mode appears smoother than those with no interpolation.

3. The script

With this new version, the script is split in two parts:
  • (i) The core functions implementing the perspective transform  and located in the file es6_perspective_transform.js.
  • (ii) A main script — warping_main.js — adding a dialog box and calling the core functions.

+++ Javascript:es6_perspective_transform.js +++
/**
* 2016-02-15
* 2018-02-15 - UPDATE for ECMAScript 2015+
* Jean-Christophe Taveau
*
* Adapted From OpenVG
* Taken from http://read.pudn.com/downloads159/sourcecode/graph/713580/ri_package_1.1/ri/src/riVGU.cpp__.htm
* VGUErrorCode RI_APIENTRY vguComputeWarpSquareToQuad(
* VGfloat dx0, VGfloat dy0,
* VGfloat dx1, VGfloat dy1,
* VGfloat dx2, VGfloat dy2,
* VGfloat dx3, VGfloat dy3,
* VGfloat * matrix
*
* from Heckbert:Fundamentals of Texture Mapping and Image Warping
* Note that his mapping of vertices is different from OpenVG's
* (0,0) => (dx0,dy0)
* (1,0) => (dx1,dy1)
* (0,1) => (dx2,dy2)
* (1,1) => (dx3,dy3)
*/
const WARP = {};
WARP.interpolateNone = (img,x,y) => {
return img.get(Math.floor(x),Math.floor(y));
}
WARP.interpolateBilinear = (img,x,y) => {
let dx = x - Math.floor(x);
let dy = y - Math.floor(y);
let pix00 = img.get(Math.floor(x),Math.floor(y));
let pix01 = img.get(Math.floor(x),Math.floor(y)+1);
let pix10 = img.get(Math.floor(x)+1,Math.floor(y));
let pix11 = img.get(Math.floor(x)+1,Math.floor(y)+1);
let R1 = (1 - dx) * pix00 + dx * pix10;
let R2 = (1 - dx) * pix01 + dx * pix11;
return (1 - dy) * R1 + dy * R2;
}
WARP.mult = (matrix,x,y) => {
let xo=matrix[0]*x+matrix[3]*y + matrix[6];
let yo=matrix[1]*x+matrix[4]*y + matrix[7];
let w =matrix[2]*x+matrix[5]*y + matrix[8];
let arr=[];
arr[0]=xo/w;
arr[1]=yo/w;
return arr;
}
WARP.computeSquareToQuad = (x,y) => {
let matrix = new Array(9);
let dx0 = x[0]; var dy0 = y[0];
let dx1 = x[1]; var dy1 = y[1];
let dx2 = x[2]; var dy2 = y[2];
let dx3 = x[3]; var dy3 = y[3];
let diffx1 = dx1 - dx3;
let diffy1 = dy1 - dy3;
let diffx2 = dx2 - dx3;
let diffy2 = dy2 - dy3;
let det = diffx1*diffy2 - diffx2*diffy1;
if (det == 0.0) {
return null;
}
let sumx = dx0 - dx1 + dx3 - dx2;
let sumy = dy0 - dy1 + dy3 - dy2;
if (sumx === 0.0 && sumy === 0.0) {
//affine mapping
matrix[0] = dx1 - dx0;
matrix[1] = dy1 - dy0;
matrix[2] = 0.0;
matrix[3] = dx3 - dx1;
matrix[4] = dy3 - dy1;
matrix[5] = 0.0;
matrix[6] = dx0;
matrix[7] = dy0;
matrix[8] = 1.0;
return matrix;
}
let oodet = 1.0 / det;
let g = (sumx*diffy2 - diffx2*sumy) * oodet;
let h = (diffx1*sumy - sumx*diffy1) * oodet;
matrix[0] = dx1-dx0+g*dx1;
matrix[1] = dy1-dy0+g*dy1;
matrix[2] = g;
matrix[3] = dx2-dx0+h*dx2;
matrix[4] = dy2-dy0+h*dy2;
matrix[5] = h;
matrix[6] = dx0;
matrix[7] = dy0;
matrix[8] = 1.0;
return matrix;
}
WARP.warp = (inp,outp,matrix,size) => {
for (let y=0;y<size;y++) {
for (let x=0;x<size;x++)
{
let xy = WARP.mult(matrix,x/size,y/size);
// Calc pixel value depending of interpolation mode
let pix = interpolate(inp,xy[0],xy[1]);
outp.set(x,y,pix);
}
}
}
view raw es6_perspective_transform.js hosted with ❤ by GitHub
+++ End of Script: es6_perspective_transform.js +++

To use the core functions, we assume that you created a directory 0Homemade/ in the ImageJ/plugins folder.

... And, here is an example of use...

+++ Javascript:warping_main.js +++
/**
* Warping image
* Jean-Christophe Taveau
* Requires the file es6_perspective_transform.js
*/
// Don't like Java `load`, rename it by `require`
const require = load;
// import module(s)
require('./plugins/0Homemade/es6_perspective_transform.js');
/************* F U N C T I O N S *************/
function dialog() {
const gd = new GenericDialog("Perspective...");
gd.addNumericField("Size of the output image: ", SIZE, 0);
gd.addChoice("Interpolation: ", ["None","Bilinear"], 0);
gd.showDialog();
if (gd.wasCanceled()) {
throw("-- End of Script --");
return;
}
SIZE = gd.getNextNumber();
let mode = gd.getNextChoiceIndex();
interpolate = (mode === 0) ? WARP.interpolateNone : WARP.interpolateBilinear;
}
/************* M A I N *************/
//1- Get information about input image
const input = IJ.getImage();
const nframes=input.getNSlices();
// Global variables
let SIZE = 240; // Default value
let interpolate;
dialog();
// 2- Get the first four corners-points
const roi = input.getRoi();
const p = roi.getPolygon();
let xpoints = [0,1,2,3].map( i => p.xpoints[i]);
let ypoints = [0,1,2,3].map( i => p.ypoints[i]);
// 3- Compute the matrix
var matrix = WARP.computeSquareToQuad(xpoints, ypoints);
//4- Apply the transform to the slices of the image/stack
const output = IJ.createImage("unwarp_"+input.getTitle(), "8-bit Black", SIZE, SIZE, nframes);
if (nframes > 1) {
const inStack = input.getImageStack();
const outStack = output.getImageStack();
for (let i=1; i<=nframes; i++) {
WARP.warp(inStack.getProcessor(i),outStack.getProcessor(i),matrix,SIZE);
if ( (i % 10) === 0) {
IJ.showProgress(i,nframes);
}
}
}
else {
WARP.warp(input.getProcessor(),output.getProcessor(),matrix,SIZE);
}
output.show();
IJ.showMessage("End of Script");
throw "End of Script";
view raw warping_main.js hosted with ❤ by GitHub
+++ End of Script: warping_main.js +++

The old implementation for archive...
+++ Javascript:perspective_transform.js +++
/******************************
* 2016-02-15
* Jean-Christophe Taveau
*
* Adapted From OpenVG
* Taken from http://read.pudn.com/downloads159/sourcecode/graph/713580/ri_package_1.1/ri/src/riVGU.cpp__.htm
* VGUErrorCode RI_APIENTRY vguComputeWarpSquareToQuad(
* VGfloat dx0, VGfloat dy0,
* VGfloat dx1, VGfloat dy1,
* VGfloat dx2, VGfloat dy2,
* VGfloat dx3, VGfloat dy3,
* VGfloat * matrix
* )
*
* from Heckbert:Fundamentals of Texture Mapping and Image Warping
* Note that his mapping of vertices is different from OpenVG's
* (0,0) => (dx0,dy0)
* (1,0) => (dx1,dy1)
* (0,1) => (dx2,dy2)
* (1,1) => (dx3,dy3)
*
*
*******************************/
//1- Get information about input image
var input = IJ.getImage();
var nframes=input.getNSlices();
var SIZE = 240;
dialog();
// 2- Get the four corners-points
var roi = input.getRoi();
var p = roi.getPolygon();
for (var i=0; i<p.npoints; i++) {
IJ.log("X:"+p.xpoints[i]+",Y:"+p.ypoints[i]);
}
// 3- Compute the matrix
var matrix=computeSquareToQuad(p.xpoints, p.ypoints);
IJ.log(matrix[0]+" "+matrix[3]+" "+matrix[6]);
IJ.log(matrix[1]+" "+matrix[4]+" "+matrix[7]);
IJ.log(matrix[2]+" "+matrix[5]+" "+matrix[8]);
//4- Apply the transform to the slices of the image/stack
var output = IJ.createImage("unwarp_"+input.getTitle(), "8-bit Black", SIZE, SIZE, nframes);
if (nframes > 1) {
var inStack = input.getImageStack();
var outStack = output.getImageStack();
for (var i=1; i<=nframes; i++) {
warp(inStack.getProcessor(i),outStack.getProcessor(i),matrix,SIZE);
if ( (i%10) == 0) {
IJ.showProgress(i,nframes);
}
}
}
else {
warp(input.getProcessor(),output.getProcessor(),matrix,SIZE);
}
output.show();
IJ.showMessage("End of Script");
throw "End of Script";
// F U N C T I O N S
function dialog() {
var gd = new GenericDialog("Perspective...");
gd.addNumericField("Size of the output image: ", SIZE, 0);
gd.addChoice("Interpolation: ", ["None","Bilinear"], 0);
gd.showDialog();
if (gd.wasCanceled()) {
throw("-- End of Script --");
return;
}
SIZE = gd.getNextNumber();
var mode = gd.getNextChoiceIndex();
if (mode == 0) {
interpolate = function (img,v0,v1) {
return interpolateNone(img,v0,v1);
}
}
else {
interpolate = function (img,v0,v1) {
return interpolateBilinear(img,v0,v1);
}
}
}
function warp(inp,outp,matrix,size) {
for (var y=0;y<size;y++) {
for (var x=0;x<size;x++)
{
var xy=mult(matrix,x/size,y/size);
// NEAREST PIXEL VALUE
var pix = interpolate(inp,xy[0],xy[1]);
outp.set(x,y,pix);
}
}
}
function interpolateNone(img,x,y) {
return img.get(Math.floor(x),Math.floor(y));
}
function interpolateBilinear(img,x,y) {
var dx = x - Math.floor(x);
var dy = y - Math.floor(y);
var pix00 = img.get(Math.floor(x),Math.floor(y));
var pix01 = img.get(Math.floor(x),Math.floor(y)+1);
var pix10 = img.get(Math.floor(x)+1,Math.floor(y));
var pix11 = img.get(Math.floor(x)+1,Math.floor(y)+1);
var R1 = (1 - dx) * pix00 + dx * pix10;
var R2 = (1 - dx) * pix01 + dx * pix11;
return (1 - dy) * R1 + dy * R2;
}
function mult(matrix,x,y) {
xo=matrix[0]*x+matrix[3]*y + matrix[6];
yo=matrix[1]*x+matrix[4]*y + matrix[7];
w =matrix[2]*x+matrix[5]*y + matrix[8];
var arr=[];
arr[0]=xo/w;
arr[1]=yo/w;
return arr;
}
function computeSquareToQuad(x,y) {
var matrix = new Array(9);
var dx0 = x[0]; var dy0 = y[0];
var dx1 = x[1]; var dy1 = y[1];
var dx2 = x[2]; var dy2 = y[2];
var dx3 = x[3]; var dy3 = y[3];
var diffx1 = dx1 - dx3;
var diffy1 = dy1 - dy3;
var diffx2 = dx2 - dx3;
var diffy2 = dy2 - dy3;
var det = diffx1*diffy2 - diffx2*diffy1;
if (det == 0.0)
return null;
var sumx = dx0 - dx1 + dx3 - dx2;
var sumy = dy0 - dy1 + dy3 - dy2;
if (sumx == 0.0 && sumy == 0.0)
{ //affine mapping
matrix[0] = dx1 - dx0;
matrix[1] = dy1 - dy0;
matrix[2] = 0.0;
matrix[3] = dx3 - dx1;
matrix[4] = dy3 - dy1;
matrix[5] = 0.0;
matrix[6] = dx0;
matrix[7] = dy0;
matrix[8] = 1.0;
return matrix;
}
var oodet = 1.0 / det;
var g = (sumx*diffy2 - diffx2*sumy) * oodet;
var h = (diffx1*sumy - sumx*diffy1) * oodet;
matrix[0] = dx1-dx0+g*dx1;
matrix[1] = dy1-dy0+g*dy1;
matrix[2] = g;
matrix[3] = dx2-dx0+h*dx2;
matrix[4] = dy2-dy0+h*dy2;
matrix[5] = h;
matrix[6] = dx0;
matrix[7] = dy0;
matrix[8] = 1.0;
return matrix;
}
view raw perspective_transform.js hosted with ❤ by GitHub
+++ End of Script: perspective_transform.js +++
Posted by Jean-Christophe T. at 17:49
Labels: , , ,

No comments:

Post a Comment

Subscribe to: Post Comments (Atom)