Landmark Warp Tools Library User Guide

Overview

This is a library for estimating and applying transformations that are defined by landmark input/output pairs; that is, a warp defined by homologous sets of points. This framework allows the user to find a mapping that best maps the landmarks in one space to landmarks in another space. There are several types of warps that differ in how the mapping is estimated. This provides tools for computing, applying and implementing the warps in an extensible way. Currently there are the following types of functions: Euclidean, Similarity, Affine and non-linear spline with basis functions. There are the following function estimation methods: Procrustes, Generalized Procrustes, Dual Quaternion and linear least squares best fit. These differ in the constraints on the mapping that is estimated. For example, there are rigid transforms that do not necessarily map landmarks exactly and minimize the error of the mapping between landmarks and non-linear splines that map landmarks exactly and maximize the smoothness between landmarks.

Implementation

This library has an object oriented design that is supported by elements of the ShapeTools library.

A landmark is a mutable object that contains N-dimensional floating point data. An indexed set of landmarks is a LandmarkSet, which does not allow (approximately) duplicate elements. Furthermore, homologous landmark sets are contained in a LandmarkSetPair object, that requires the sets to have the same size and dimension. A landmark set pair is then used to estimate a warp object.

Each warp is an immutable object that can take a point from the input space and compute the output point. By convention, the input space is called the 'reference' space and the output space is called the 'target' space. Warp objects implement an interface IFunction. Like an IFunction, IFunction objects favor memory reuse by copying output to an argument as opposed to allocating and returning a result.

Although IFunction objects allow for only point based operations, given structures that contain multiple points, mapping objects can be used to apply the function to all the points in the structure without changing the transformation implementation. An example of a mapping function is one that implements IMap. This allows a minimal amount of implementation for the specific transformation. Furthermore, the process of applying the transformation to multiple points is easily reimplemented, i.e. for multi-threaded or distributed computing.

Each warp estimation method is an immutable object that takes a pair of landmark sets and returns the IFunction that best maps one to the other. The interface to an estimation method is IFunctionEstimator. There are factories for the estimation routines. The user specifies which estimation method should be returned using an enumeration of the implementing objects.

To allow IO operations, an object must have a parameterization object that implements IParameterization. This parameterization object allows the user to create a ParameterBlock that encodes the instance of the object. The class FunctionIO can be used to write a warp to disk automatically. The file will encode all the parameters necessary to recreate the warp object. This file is ASCII and human readable, although in general it is difficult to interpret the parameters.

There is also a file format for landmarks. This allows LandmarkSet objects to be read and written to disk. The format is ASCII and human readable and writable. Each point has space-delimited coordinates, and points are line-delimited. The order of the points in the stream implies the indices in the LandmarkSet. Comments are allowed after a '#' character, and arbitrary white-space and empty lines are allowed.

Mathematical Background

This section will address the mathematics of the warps.

For an affine warp, the transformation is the combination of an arbitrary linear transformation and a translation. This allows for translation, scaling, rotation and skewing. This is estimated by solving the linear least squares problem. Given two sets of landmarks that are related in a non-linear way, this will not map landmarks exactly. For an N-dimensional problem, at least N of the K landmarks must be in general position.

There is also a rigid warp, which is a combination of a translation and a rotation. There are multiple ways to estimate the parameters of this transformation. The Procrustes analysis estimates the best orthogonal rotation matrix and Dual Quaternion solves the problem in the domain of quaternions and translates the result to the euclidean domain. In addition, scaling may be considered when computing the warp. The generalized Procrustes method estimates the rigid with scaling warp. For a N-dimensional problem, at least N of the K landmarks must be in general position.

For the non-linear spline case, the function is described as the sum of an affine transform and the weighted sum of nonlinear basis functions. Each warp of this kind has a basis function of a single form that is modeled by an IBasisFunction. For a N-dimensional problem, at least N of the K landmarks must be in general position.

Specifically, consider two spaces A and B, and let there be two sets of points p_i from A and q.i from B for i = [0,K]. For our purposes, A and B will both be RN. We would like to find a function f:A->B such that f(p_i) = q_i for all i. This is the mapping we call the warp. Let G be a function G:RN->RNxN which maps a point to a NxN matrix-- this is the basis function. Let d be a function d:R^NxR^N->R^N that is the displacement function d(a,b) = b-a, and f(a) = a + d(a,b). We can define d as a linear combination of basis functions and a linear transformation:

d(p) = sum( G(p-p_i) * c_i ) + A * p + b

where the sum is from i = 0 to i = K. P is from R^N and is the point to be mapped. p_i from RN is the i-th landmark in the domain. c_i is from RN and is the i-th non-linear coefficient. A is from R^Nx^N and contains the affine coefficients. b is from R^N and contains the linear shift coefficients.

The coefficients are computed by solving a system of K*N+N*(N+1) linear equations which ensure the function maps landmarks to landmarks in addition to some affine transform.

Warp File Format

This is a file format that is compatible with all landmark-based warps that are IParameterizable that are recognized by the LandmarkWarpFactory class.

All data is written as ascii text. The format allows comments on any line and after data, as long as they are preceded by a '#' character (text after that symbol is not parsed). A valid file can have any number of comments. All arrays of numbers are space-delimited with a newline every 80 characters.

Format

The first line of data contains a magic string, "LandmarkWarpFramework", followed by the version of the framework. Next, there is a line indicating the number of parameter blocks that have been encoded. Each parameter block first encodes a header which is a single line containing a name, version and the number of parameters in the block. Next the parameter block encodes the parameters, adding a new line every 80 characeters. For most warps, only one parameter block is needed, but for spline warps, two are needed, one to specify the basis and one to specify the coefficients of the basis.

The meaning of the parameters is very specific to the type of warp. For details about the meaning of the parameters, the javadocs should be referenced.

The filenames have the form "*.lmw".

Example Warp File

Here is an example warp file:

# LandmarkWarpFile 2007-04-03
LandmarkWarpFramework 1.0 
ParameterBlock 2
SplineWarp 1.0 28
2.0 5.0 5.0 5.0 -5.0 -5.0 0.0 5.0 -5.0 -5.0 5.0 0.0 0.0060112295 -0.0 
0.0060112295 -0.0 0.0060112295 -0.0 0.0060112295 -0.0 -0.024044918 0.0 -8.0349545E-17 -0.0 
-1.3226087E-16 -0.0 7.2032137 -0.0 
ThinPlateBasis 1.0 1
2.0 

Landmark File Format

This is a file format that is compatible that encodes K (approximately) unique N-dimensional landmarks to a file.

All data is written as ascii text. The format allows comments on any line and after data, as long as they are preceded by a '#' character (text after that symbol is not parsed). The format also allows arbitrary white space and empty lines

The format is intended to be human readable and compatible with other applications.

Format

Each point is line delimited, where the order of the points in the file reflects the point indices in the LandmarkSet. The coordinates are space delimited.

The filenames have the form "*.lm".

Example Landmark File

Here is an example landmark file:

56.0 49.0
203.0 33.0
243.0 69.0
234.0 93.0
147.0 122.0
123.0 96.0
128.0 77.0
85.0 64.0
86.0 55.0

which is equivalent to

56.0 49.0 
203.0 33.0 # first comment
243.0 69.0

234.0 93.0
# second comment
147.0 122.0
  123.0 96.0
128.0  77.0
85.0 64.0
86.0 55.0

Example Use

Here is a way to make a simple 2-D spline warp which deforms the center of a square to the right:

import java.awt.image.BufferedImage;
import java.io.BufferedReader;
import java.io.File;
import java.io.FileInputStream;
import java.io.FileNotFoundException;
import java.io.FileOutputStream;
import java.io.FileReader;
import java.io.IOException;
import java.io.InputStreamReader;
import java.util.Random;

import edu.ucla.loni.ccb.discretedata.io.ArrayIO;
import edu.ucla.loni.ccb.functions.app.FunctionAppUtil;
import edu.ucla.loni.ccb.functions.basis.BasisElasticBodyHyper;
import edu.ucla.loni.ccb.functions.basis.BasisFunctionType;
import edu.ucla.loni.ccb.functions.estimation.FunctionEstimatorFactory;
import edu.ucla.loni.ccb.functions.estimation.LandmarkOperations;
import edu.ucla.loni.ccb.functions.estimation.LandmarkSet;
import edu.ucla.loni.ccb.functions.estimation.LandmarkSetPair;
import edu.ucla.loni.ccb.functions.estimation.SplineBasisEstimator;
import edu.ucla.loni.ccb.functions.io.FunctionIO;
import edu.ucla.loni.ccb.shape.Shape;
import edu.ucla.loni.ccb.shape.ShapeImpl;
import edu.ucla.loni.ccb.shape.geometry.CGALFaceSet;
import edu.ucla.loni.ccb.shape.geometry.IPoint;
import edu.ucla.loni.ccb.shape.geometry.IPointSet;
import edu.ucla.loni.ccb.shape.geometry.IrregularPointSet;
import edu.ucla.loni.ccb.shape.geometry.Point3D;
import edu.ucla.loni.ccb.shape.geometry.PointND;
import edu.ucla.loni.ccb.shape.math.VectorMath;
import edu.ucla.loni.ccb.shape.math.functions.IFunction;
import edu.ucla.loni.ccb.shape.math.functions.MapLoop;
import edu.ucla.loni.ccb.shapeio.duff.DuffIO;

public class Simple2DLandmarkWarpExample
  {
    /**
     * This demonstrates how to apply a warp to a single two dimensional point.
     * 
     * The program creates a simple warp, applies it to
     * a new point and writes the warp to disk.
     */
    public static void main(String [] args)
      {
        try {

          // Define example landmarks
          int dim = 2;
          int num = 5;
          float distance = 5;

          float[][] refData = new float[dim][num];
          float[][] targetData = new float[dim][num];

          float[] novelPoint = { .2f, 1.4f };
          float[] warpedNovelPoint = new float[dim];

          refData[0] = new float[] { distance, distance, -distance, -distance,
              0 };
          refData[1] = new float[] { distance, -distance, -distance, distance,
              0 };

          targetData[0] = new float[] { distance, distance, -distance,
              -distance, distance / 2 };
          targetData[1] = new float[] { distance, -distance, -distance,
              distance, 0 };

          // Create the landmark objects
          LandmarkSet refSet = LandmarkOperations
              .constructFromDimNumArray(refData);
          LandmarkSet targetSet = LandmarkOperations
              .constructFromDimNumArray(targetData);
          LandmarkSetPair landmarks = new LandmarkSetPair(refSet, targetSet);

          // Compute the warp

          IFunction warp = FunctionEstimatorFactory.createNonLinear(
              BasisFunctionType.THIN_PLATE, dim).estimate(landmarks);

          // Evaluate the warp (these are not used for anything here)
          warp.evaluate(novelPoint, warpedNovelPoint);
          System.out.println("Original warped novel point = ["
              + warpedNovelPoint[0] + "," + warpedNovelPoint[1] + "].");

          // Write out the warp file
          String filename = "output/warp" + FunctionIO.SUFFIX;
          FileOutputStream stream = new FileOutputStream(new File(filename));
          FunctionIO.write(stream, warp);

          // Read warp file
          BufferedReader reader = new BufferedReader(new FileReader(new File(
              filename)));
          IFunction readWarp = FunctionIO.read(reader);

          // Demonstrate that the read warp produces the same output at the
          // novel point
          readWarp.evaluate(novelPoint, warpedNovelPoint);
          System.out.println("Read and recomputed warped novel point = ["
              + warpedNovelPoint[0] + "," + warpedNovelPoint[1] + "].");

        } catch (Exception e) {
          System.err.println("Warp write failed.");
          e.printStackTrace();
        }
      }
  }

The above code creates two two-dimensional arrays, a set of source landmarks and a set of destination landmarks. The first four points are the same-- they are the corners of a square. The last point is different by distance/2 in the x-coordinate. Then the warp is created from these points. The warp is applied to a new point. This warp should keep points near the corners stationary and move points near the center to the right. Finally the warp is written to disk for later use.

Here is the resulting warp file:

# LandmarkWarpFile 2007-04-03
LandmarkWarpFramework 1.0 
ParameterBlock 2
SplineWarp 1.0 28
2.0 5.0 5.0 5.0 -5.0 -5.0 0.0 5.0 -5.0 -5.0 5.0 0.0 0.0060112295 -0.0 
0.0060112295 -0.0 0.0060112295 -0.0 0.0060112295 -0.0 -0.024044918 0.0 -8.0349545E-17 -0.0 
-1.3226087E-16 -0.0 7.2032137 -0.0 
ThinPlateBasis 1.0 1
2.0 

There is an example package in the library which provides landmarks for testing as well as main routines for applying the warps to existing shapes.