# Attachments file level help

## C:\Users\AB\DArT_Toolshed\Attachments\FEMLib

femlib_GQP.m

 function [GQP, w_gauss] = femlib_GQP

Format  each row is an xi, eta, zeta component

Dr. A. I. Hanna (2006)


femlib_GQPwedge.m

 function [GQP, w_gauss] = femlib_GQPwedge

Format  each row is an xi, eta, zeta component

Dr. A. I. Hanna (2006)


femlib_N.m

 function Nprime = femlib_NGrad(xi, eta, zeta)

Calculates the derivative matrix of shape functions with respect to each
natrual co-ordinate.

The trial and test functions are given by

N1 = 0.5*xi*(1+zeta);
N2 = 0.5*eta*(1+zeta);
N3 = 0.5*(1-xi-eta)*(1+zeta);
N4 = 0.5*xi*(1-zeta);
N5 = 0.5*eta*(1-zeta);
N6 = 0.5*(1-xi-eta)*(1-zeta);
_                                        _
N = |N1 0 0 N2 0 0 N3 0 0 N4 0 0 N5 0 0 N6 0 0 |
|0 N1 0 0 N2 0 0 N3 0 0 N4 0 0 N5 0 0 N6 0 |
|0  0 N1 0 0 N2 0 0 N3 0 0 N4 0 0 N5 0 0 N6|
-                                        -

Dr. A. I. Hanna (2006).


 function Nprime = femlib_NGrad(xi, eta, zeta)

Calculates the derivative matrix of shape functions with respect to each
natrual co-ordinate.

The trial and test functions are given by

N1 = 0.5*xi*(1+zeta);
N2 = 0.5*eta*(1+zeta);
N3 = 0.5*(1-xi-eta)*(1+zeta);
N4 = 0.5*xi*(1-zeta);
N5 = 0.5*eta*(1-zeta);
N6 = 0.5*(1-xi-eta)*(1-zeta);
_                                   _
Nprime = |dN1/d_xi   dN2/d_xi ....  dN6/d_xi  |
|dN1/d_eta  dN2/d_eta      dN6/d_eta |
|dN1/d_zeta dN2/d_zeta     dN6/d_zeta|
-                                  -

Dr. A. I. Hanna (2006).


 function mesh = femlib_addRandomZ(mesh)

Adds random noise to the z-axis to allow bending into this plane.

Dr. A. I. Hanna (2006)


femlib_calcResidualStress.m

 function varargout = femlib_calcResidualStress(varargin)

Dr. A. I. Hanna (2006)


femlib_constructB.m

 function B = femlib_constructB(G)

Here we are constructing the B matrix used in the integral (B'DB)

Recall:
-                                   -
|dN1/d_xi   dN2/d_xi   .... dN6/d_xi  |
Nprime =  |dN1/d_eta  dN2/d_eta  .... dN6/d_eta |
|dN1/d_zeta dN2/d_zeta .... dN6/d_zeta|
-                                   -
and

J = Nprime*Xe;
G = inv(J)*Nprime;

so
-                            -
|dN1/d_x  dN2/d_x ....  dN6/d_x|
G = |dN1/d_y  dN2/d_y ....  dN6/d_y|
|dN1/d_z  dN2/d_z ....  dN6/d_z|
-                            -

-                                                                                       -
|dN1/d_x    0       0     dN2/d_x       0        0        .... dN6/d_x      0          0  |
|   0    dN1/d_y    0         0     dN2/d_y      0        ....    0      dN6/d_y       0  |
B = |   0       0    dN1/dz       0         0     dN2/d_z     ....    0         0      dN6/dz |
|   0    dN1/dz  dN1/dy       0      dN2/dz    dN2/dy     ....    0      dN6/dz    dN6/dy |
|dN1/dz    0     dN1/dx     dN2/dz      0      dN2/dx     .... dN6/dz       0      dN6/dx |
|dN1/dy  dN1/dx     0       dN2/dy   dN2/dx       0       .... dN6/dy     dN6/dx        0 |
-                                                                                       -

and our strain vector has elements [e_x, e_y, e_z, e_yz, e_xz, e_xy]'

Dr. A. I. Hanna (2006).


femlib_createBasicMesh.m

 function mesh = femlib_createBasicMesh(pts, T)

Generates a default mesh structure

Dr. A. I. Hanna (2006)


femlib_eliminateEquations.m

 function [K,f, usedEquations] = femlib_eliminateEquations( K, f, fixedNodes )

Here we remove the equations associated with fixedNodes from the global
stiffness matrix and global force vector.

Dr. A. I. Hanna (2006).


femlib_findLongestEdge.m

femlib_findNeighbour.m

femlib_fixTriMesh.m

 function T = femlib_fixTriMesh(T, X)

femlib_fixTriMesh takes in a triangulation and the set of vertices, then for each
triangle it reorders the indices to make sure that they are ordered
clockwise.

Dr. A. I. Hanna (2006)


femlib_getElements.m

 function element = femlib_getElements(pts)

Here we need to do a few things, firstly we need to get the
triangulation, and then we need to make sure that all the triangles are
ordered clockwise. Then we are set.

Dr. A. I. Hanna (2006).


femlib_getFixedEquations.m

 function fixedEquations = femlib_getFixedEquations(dorsalDivide)

We we are returning the indices to the equations we want to keep fixed

Dr. A. I. Hanna (2006).


femlib_getMaterialStiffnessMatrix.m

 function D = femlib_getMaterialStiffnessMatrix(E, v)

This function takes the Youngs modulus (E) and Poissons ratio (v) and
constructs a stiffness matrix. It assumes that the materials in the model
are make of homogenous material.

Dr. A. I. Hanna


femlib_getPreStrain.m

 function preStress = femlib_getPreStrain

Dr. A. I. Hanna (2007)
recall that the stress is defined by
stress = [du/dx, dv/dy, dw/dz, du/dy + dv/dx, du/dz + dw/dx, dv/dz +dw/dy];
preStrain = [-1 0 0 0 0 0]'; % lateral growth
preStrain = [0 -1 0 0 0 0]'; % longitudal growth
preStrain = [0 0 -1 0 0 0]'; % vertical growth
preStrain = [0 0 0 -1 0 0]'; % pure shear yz-plane
preStrain = [0 0 0 0 -1 0]'; % pure shear xz-plane
preStrain = [0 0 0 0 0 -1]'; % pure shear xy-plane


femlib_getPreStress.m

 recall that the stress is defined by
stress = [du/dx, dv/dy, dw/dz, du/dy + dv/dx, du/dz + dw/dx, dv/dz + dw/dy];


femlib_orderVertices.m

 function [rpts] = femlib_orderVertices(pts)

Orders the vertices given in pts in a clockwise order.

Dr. A. I. Hanna (2006).


femlib_plotMesh.m

 function ph = femlib_plotMesh(X, element, dorsalDivide)

Inputs:
X - a Nx2 matrix of global vertices
element - an array of structures of elements
dorsalDivide - the index that seperates the dorsal from the ventral

Outputs:
ph - the cell array of plot handles

Dr. A. I. Hanna (2006).


femlib_plotWedge.m

 function varargout = femlib_plotWedge(varargin)

Inputs:
'pts' - a matrix of x, y, z, data points
'edgecolor' - the color of the wedge edges
'axish' - the specified axis handle (default = gca).

A drawing routine for displaying wedge shaped finite elements. The matrix
associated with the parameter 'pts' must but in a certain order. we have
the upper triangle points specified then the lower triangle points. See
the example below.

Example:
X = [1 0 1;
0 1 1;
0 0 1;
1 0 -1;
0 1 -1;
0 0 -1];

femlib_plotWedge('PTS', X, 'edgecolor', 'g');

Dr. A. I. Hanna (2006)


femlib_setUniformTemp.m

femlib_splitElements.m

 function [mesh] = femlib_splitElements(mesh)

Dr. A. I. Hanna (2007).


femlib_splitWedgeElement.m

 function mesh = femlib_splitWedgeElement(mesh, T_target, E_target)

Dr. A. I. Hanna (2007)


femlib_stressWedge.m

 function varargout = femlib_stressWedge(varargin)

Dr. A. I. Hanna (2006)


femlib_subdivideWedgeElement.m

femlib_tridata2obj.m

 function femlib_tridata2obj(varargin)

Dr. A. I. Hanna (2007)


femlib_updateMeshPlot.m

 function femlib_updateMeshPlot(X, elements, dorsalDivide, ph)

Updates the plot created by femlib_plotMesh, rather than replot the mesh
we simply replace the X, Y, Z, data in the plot.

Dr. A. I. Hanna (2006).


femlib_wedgeArea.m

femlib_wedgeIteration.m

 function X = femlib_wedgeIteration(X, meshInfo)

A script to perform a single iteration of the FEM procedure for growing a
body defined by wedge elements.

Dr. A. I. Hanna (2006)


## C:\Users\AB\DArT_Toolshed\Attachments\GTLib

gtlib_calcGrowthParameters.m

function [m1, m2, tri, growthparams] = gtlib_calcGrowthParameters(m1, m2, N, tri)

Dr. A. I. Hanna (2007)


gtlib_demo_GT_example.m

 function gtlib_demo_GT_example

This is a script that takes a direction of growth together with principle
and minor growth rate and computes the growth tensor.

Dr. A. I. Hanna (2007)
Calculate the growth tensor


gtlib_demo_GT_example2.m

 function gtlib_demo_GT_example2

This is a script that takes a direction of growth together with principle
and minor growth rate and computes the growth tensor.

Dr. A. I. Hanna (2007)
Calculate the growth tensor


gtlib_demo_GT_linearity_example.m

 function gtlib_demo_GT_linearity_example

This is a script that takes a direction of growth together with principle
and minor growth rate and computes the growth tensor.

Dr. A. I. Hanna (2007)


gtlib_growth3dTriangle.m

GROWTH3DTRIANGLE Compute the growth of a triangle in 3D space
Inputs are:
p = 3*3 matrix where each line is a point at time t0
q = 3*3 matrix where each line is a point at time t1
Outputs:
kmax, kmix - Growth along the major and minor axis
umax, umin - Unit vectors defining the major and minor axis

Dr. P. Barbier de Reuille (2009)


gtlib_growthParams2Tensor.m

 function T = gtlib_growthParams2Tensor(r1, r2, alpha)

A method that takes the direction of principle growth, together with the
growth rates along that direction, and along the direction perpendicular
to that direction and produces a growth tensor T.

Example
Calculate the growth tensor
GT = gtlib_growthParams2Tensor(2, 4, (3*pi)/9);

For rigour we perform the reverse process
[major, minor, theta] = gtlib_growthTensor2Params(GT)

To visualize the growth tensor we make the unit circle
theta = linspace(0, 2*pi, 200);
X = [cos(theta); sin(theta); zeros(size(theta))];
And apply the growth tensor to the data
D = GT*X;
%The Art Corner

figure(1); clf; hold on;
plot3(X(1,:), X(2,:), X(3,:), 'b', 'LineWidth', 2);
plot3(D(1,:), D(2,:), D(3,:), 'r', 'LineWidth', 2);
V = D-X;
N = size(X,2);
m = 50;
ind = 1:N/m:N;
quiver3(X(1,ind), X(2,ind), X(3,ind), V(1,ind), V(2,ind), V(3,ind),0);
grid on;
view(2);
axis image;
xlabel('X'); ylabel('Y')
title('Example of an Anisotropic Growth Tensor')

Dr. A. I. Hanna (2007)


gtlib_growthTensor2Params.m

 function [major, minor, theta] = gtlib_growthTensor2Params(GT)

A method that takes a growth tensor T and returns the direction of
principle growth together with the growth rate along that direction and
the growth rate perpendicular to that direction.

Example
Calculate the growth tensor
GT = gtlib_growthParams2Tensor(2, 4, (3*pi)/9);

For rigour we perform the reverse process
[major, minor, theta] = gtlib_growthTensor2Params(GT)

To visualize the growth tensor we make the unit circle
theta = linspace(0, 2*pi, 200);
X = [cos(theta); sin(theta); zeros(size(theta))];
And apply the growth tensor to the data
D = GT*X;
%%%%%%%%%%%%%%%%%%%%

The Art Corner

%%%%%%%%%%%%%%%%%%%%
figure(1); clf; hold on;
plot3(X(1,:), X(2,:), X(3,:), 'b', 'LineWidth', 2);
plot3(D(1,:), D(2,:), D(3,:), 'r', 'LineWidth', 2);
V = D-X;
N = size(X,2);
m = 50;
ind = 1:N/m:N;
quiver3(X(1,ind), X(2,ind), X(3,ind), V(1,ind), V(2,ind), V(3,ind),0);
grid on;
view(2);
axis image;
xlabel('X'); ylabel('Y')
title('Example of an Anisotropic Growth Tensor')

Dr. A. I. Hanna (2007)


gtlib_interpTriangle.m

 function t = gtlib_interpTriangle(t)

Dr. A. I. Hanna (2007)


gtlib_interp_tri_pts.m

 function pts = gtlib_interp_tri_pts(pts)

A function that takes a 3x2 matrix that defines a triangle and
interpolates along the edges.

Input Params:
pts - a 3x2 matrix of points

Output Params:
pts - a 6x2 matrix of points

Dr. A. I. Hanna (2006).

Used to calculated warps that need at least 6 points for a particular triangle


gtlib_interpolateDisplacementField.m

function [X, w] = gtlib_interpolateDisplacementField(m1, m2, Nx, Ny)

Dr. A. I. Hanna (2007)


gtlib_plotAAMTriangleGrowth.m

 version 2


gtlib_plotAAMTriangleGrowth_MidveinPlane3.m

 version 2


gtlib_plotAAMTriangleGrowth_kmaxkminkarea.m

 version 2


gtlib_plotGrowthTensor.m

 function varargout = gtlib_plotGrowthTensor(varargin)

A function that draws the growth tensor, input parameters are

'growth_tensor' - a 3x3 growth tensor
'offset' - a 2x1 vector with an x & y offset
'colour' - the colour of the edge of the tensor
'Parent' - the handle of the axis where to draw the tensor
'scale' - the scale of the tensor
'line_width' - the width of the plot lines

Returns the plot handle to the tensor

Dr. A. I. Hanna (2007)


gtlib_plotGrowthTensorCross.m

 function varargout = gtlib_plotGrowthTensorCross(varargin)

A function that draws the growth tensor cross, input parameters are

'growth_tensor' - a 3x3 growth tensor
'offset' - a 2x1 vector with an x & y offset
'colour' - the colour of the edge of the tensor
'Parent' - the handle of the axis where to draw the tensor
'scale' - the scale of the tensor

Returns the plot handle to the tensor

Dr. A. I. Hanna (2007)


gtlib_plotPrincipalDirection.m

 function varargout = gtlib_plotPrincipalDirection(varargin)

A function that draws the growth tensor cross, input parameters are

'growth_tensor' - a 3x3 growth tensor
'offset' - a 2x1 vector with an x & y offset
'colour' - the colour of the edge of the tensor
'Parent' - the handle of the axis where to draw the tensor
'scale' - the scale of the tensor
'linwidth' - the width of the line

Returns the plot handle to the tensor

Dr. A. I. Hanna (2007)


gtlib_svd2growthparams.m

 function growthparams = gtlib_svd2growthparams(g1, g2, t)

before we would have done
smax = 1/S(1,1);
smin = 1/S(2,2);
H = (1./smax).*(1./smin);
G = t.*((log(2))./(log(H)));
A = ((1./smax)./(1./smin)).^(G./t);

and then returned [G, A, theta], and then we would have used

doublex = sqrt(2*A);
doubley = doublex/A;
rate_x = log2(doublex)/G;
rate_y = log2(doubley)/G;

to calculate the stretch ratios we now calculate these stretch ratios from the principle strains (e1, e2) by the fact that
log2(e1)/t = log2(sqrt(2*A))/G

Proof:

log2(e1)/t = log2(sqrt(2A))/G
= (log2(sqrt(2)) + log2(sqrt(A)))/G
= (0.5 + 0.5*log2(A))/G
= (0.5 + 0.5*log2((e1/e2)^(G/t)))/G
= (0.5 + (G/(2*t))*log2(e1/e2))/G
= 1/(2*G) + (1/(2*t))*log2(e1/e2)
= 1/(2*G) + (1/(2*t))*(log2(e1) - log2(e2))
= (log2(H) + log2(e1) - log2(e2))/(2*t)
= (log2(e1) + log2(e2) + log2(e1) - log2(e2))/(2*t)
= log2(e1)/t.

In a similar way we can prove that

log2(e2)/t = log2(sqrt(2*A)/A)/G

and so our stretch ratios are given as

[log2(S(1,1))/t, log2(S(2,2))/t, theta]

Dr. A. I. Hanna (2007)


gtlib_vertexCellData2GrowthParams.m

function varargout = gtlib_vertexCellData2GrowthParams(varargin)

Example:

X = data.X;
X = {X{1:5:end}};
GROWTHSTRUCT = gtlib_vertexCellData2GrowthParams('X', X);
save('testgrowthdata.mat', 'GROWTHSTRUCT');

Dr. A. I. Hanna (2006)


IsoGaussn.m

 function G = IsoGaussn(varargin)

Generates an isotropic gaussian of dimension n.

Inputs:
'offset' - the center of the gaussian
'sigma' - the variance of the gaussian
'X' - the input data where rows are the number of samples and each column
is a dimension.

Outputs
'G' - the value of the gaussian

Dr. A. I. Hanna (2008)


## C:\Users\AB\DArT_Toolshed\Attachments\MCTLib\Data

mctj_demo.m

 function mctj_demo(id)

A method that demonstrates the use of the motion coherent velocity field
estimation methods. Here we generate

Inputs:
id - 1, 2, 3 (various demos)

Example 1:
mct_demo;

Example 2:
mct_demo('id', 1);

Example 3:
mct_demo('id', 2);

Example 4:
mct_demo('id', 3);

Dr. A. I. Hanna, CMP, UEA, 2008.


mctj_initialize_java.m

 function mctj_initialize_java

Initalize java for motioncoherence

Inputs:
mctroot - the root directory, (it must contain motioncoherence/motioncoherence.jar)

Dr. A. I. Hanna, CMP & JIC, 2008


mctlib_CalcParam.m

 function varargout = mctlib_CalcParam(varargin)

Calculates the parameters for the smooth regularized velocity field using
in MCT.

Inputs:
'x0' - the resting points
'f' - the known values of the function being approximates
'sigma' - controls the high pass filtering (large values make for a
smoother field (default = 1))
'lambda' - the regularization parameters (default = 1)

Dr. A. I. Hanna (2007)


mctlib_CalcVel.m

 function varargout = mctlib_CalcVel(varargin)

Calculates the velocity for the smooth regularized velocity field using
in MCT.

Inputs:
'x0' - the resting points
'y0' - the points after the velocity field has been applied
'sigma' - controls the high pass filtering (large values make for a
smoother field (default = 1))
'beta' - the parameters calculated from mctlib_CalcParam

Dr. A. I. Hanna (2007)


mctlib_GaussGram.m

 function varargout = mctlib_GaussGram(varargin)

Calculates the Gaussian Gram matrix for the two vectors V1 and V2 where
V1 is a NxD matrix and V2 is a MxD matrix and D is the dimension.
The output G is a NxM matrix where each element G(i,j) = dot(V1(i,:), V2(j,:))

Inputs:
'v1' - the first set of points
'v2' - the second set of points
'sigma' - controls the high pass filtering (large values make for a
smoother field (default = 1))

Dr. A. I. Hanna (2007)


mctlib_denormalise.m

 function varargout = mctlib_normalise(varargin)

Normalises a data set so it has zero mean and unit scale.

Inputs:
'pts' - the points to normalise

Dr. A. I. Hanna (2007)


mctlib_findClosest.m

 function cind = mctlib_findClosest(X, Y)

For each point in the warped template point set, find the closest point
in the reference set and make the correspondance.

Dr. A. I. Hanna, CMP, 2008.


mctlib_normalise.m

 function varargout = mctlib_normalise(varargin)

Normalises a data set so it has zero mean and unit scale.

Inputs:
'pts' - the points to normalise

Dr. A. I. Hanna (2007)


mctlib_plotAssignment.m

 function ph = mctlib_plotAssignment(varargin)

Plots the assignment between two point sets

Dr. A. I. Hanna, CMP, 2008.


mctlib_regression_demo.m

 function mctlib_regression_demo

A function that shows the method of regularization using a smoothness
term based on a high-pass filtering of the function under question. We
are searching for a function that is smooth in the sense that its high
frequency component has low energy.

Refs: Motion Coherence Theory (MCT)

Dr. A. I. Hanna (2007)


mctlib_warp_demo.m

 function mctlib_warp_demo

This simple matlab script loads some leaf point models and warps one
pointmodel onto the other. The script also generates a mesh covering the
pointmodels and shows the warped grid.

Dr. A. I. Hanna, CMP, UEA, Norwich, UK, 2008.


mctlib_warpptsxy.m

 function xw = mctlib_warpptsxy(varargin)

Inputs:
'xx' - the NxD matrix of data to warp, where D is the number of
dimensions, and N is the number of samples.
'x0' - the MxD matrix of reference points
'y0' - the MxD matrix of template points
'lambda' - the regularization parameters (default = 1)
'sigma' - the smoothness parameters (default = 1)

Outputs:
'xw' - the warped input 'xx' matrix

Dr. A. I. Hanna, CMP, UEA, Norwich, UK, 2008


## C:\Users\AB\DArT_Toolshed\Attachments\PCALib

pcalib_GPA.m

 function [X] = pcalib_GPA(varargin)

Generalized procrustes analysis

figure(1); clf;
subplot(1,2,1);
template = [0 1 1 0; 0 0 1 1];
template = template(:);
X = repmat(template, 1, 100);
X = X + rand(size(X))/10;
D = pcalib_randDeformMatrix2D(X);
pcalib_plotshapematrix(D, 'g', gca);
opts.scaling = 0;
opts.rotation = 1;
opts.translation = 1;
[X, template] = pcalib_GPA('data',D, 'template', template, 'opts', opts);
subplot(1,2,2);
pcalib_plotshapematrix(X, 'r', gca);

Dr. A. I. Hanna (2006)


pcalib_GPAMeanScale.m

 function [X, template, meanScale] = pcalib_GPAMeanScale(varargin)

Generalized procrustes analysis

figure(1); clf;
subplot(1,2,1);
template = [0 1 1 0; 0 0 1 1];
template = template(:);
X = repmat(template, 1, 100);
X = X + rand(size(X))/10;
D = pcalib_randDeformMatrix2D(X);
pcalib_plotshapematrix(D, 'g', gca);
opts.scaling = 0;
opts.rotation = 1;
opts.translation = 1;
[X, template] = pcalib_GPA('data',D, 'template', template, 'opts', opts);
subplot(1,2,2);
pcalib_plotshapematrix(X, 'r', gca);

Paul Southam 2011


pcalib_calc_deviation.m

pcalib_centershape.m

 function [x, mu] = pcalib_centershape(x, dim)

Dr. A. I. Hanna (2007)


pcalib_demo_001.m

pcalib_estimate_template.m

 function [template, temp_centroid] = estimate_template(X, dim)

Dr. A. I. Hanna (2007);


pcalib_pca.m

 function [Xm, P, b, pcaDat] = pcalib_pca(X, v)
A function that performs PCA on a matrix of data, where each column is an
experiment, and each row corresponds to points in that experiment.

re-arranged Dr. A. I. Hanna (2005).


pcalib_plot3Vector.m

%%%%%

%%%%%%


pcalib_plotshape.m

 function pcalib_plotshape(x, c, axis_h)

A simple script to plot the shape given in the column vector x.

Inputs:
x - a column vector of the form x = [x_1, y_1, x_2, y_2, ..., x_N, y_N]';
c - the color to plot the shape (default: c = [0 0 1])
axis_h - the handle of the axis to plot the shape (default: axis_h = gca)

Dr. A. I. Hanna (2006)


pcalib_plotshapematrix.m

 function pcalib_plotshapematrix(X, c, axis_h)

A simple script to plot the shapes given in the matrix X.

Inputs:
X - a matrix where each column is of the form x = [x_1, y_1, x_2, y_2, ..., x_N, y_N]';
c - the color to plot the shape (default: c = [0 0 1])
axis_h - the handle of the axis to plot the shape (default: axis_h = gca)

Dr. A. I. Hanna (2006)


pcalib_procrustes.m

 function [xaligned,ProcrustesAlignment, theta]=pcalib_procrustes(x,opts, template_pts)

This function takes in the raw point models and then aligns them
according to the options in opts. These options are a 1x3 array
containing zeros or ones. For example opts = [1 1 1] will align the
shapes according the scale, rotation and translation, so the structure of
opts is as follows, opts = [toggle_scale, toggle_rotation, toggle_translation];
If you want to align your shapes to make rotations and translations
invariant but you want to capture scaling, then use opts = [0 1 1];

The function returns 2 parameters, xaligned which are the aligned shapes,
and ProcrustesAlignment which is a 4x1 array of numbers, with the
following format

ProcrustesAlignment = [scale*cos(theta), scale*sin(theta), trans_x, trans_y]';

Inputs:
x - an NxM matrix where N is the number of parameters and M is the
number of samples. REMEMBER: each column of x is a sample, this
column representation of the sample is of the form [x1, y1, x2, y2,...]';

opts - this is the options vector (default = [1 1 1])

template_pts - this is a column vector of the form [x1, y2, ... etc
that define a template specified by the user.

Dr. A. I. Hanna (2005);


pcalib_project2tangentspace.m

 function [X] = pcalib_project2tangentspace(X, meanshape)

A function that maps the aligned shapes in X to the tangent space given
by meanshape.

Inputs:
X - a MxN matrix where M is the dimension of the shape, and N is the number of samples.
NOTE: X is of the form X = [x_1, y_1, x_2, y_2, ..., x_N, y_N]'.
meanshape - this is the meanshape in question that defines the
projection onto the tangent space. It is a column vector with the same
structure as the columns in X.
Outputs:
X - the same structure as the input X, but the shapes have now
been projected into the tangent space.

Dr. A. I. Hanna (2006).


pcalib_randDeformMatrix2D.m

 Dr. A. I. Hanna (2006)


pcalib_rotateshape.m

 function [x, R] = pcalib_rotateshape(x, template)

Takes a general shape x and a template shape t. The function returns a shape y and
a matrix R such that arg min ||y - template||

Dr. A. I. Hanna (2007)


pcalib_scaleshape.m

 function [x] = pcalib_scaleshape(x, temp_scale)

metric here is sqrt{ sum[ (x-xbar)^2 + (y-ybar)^2 ] }
since we have already centered the shape, (i.e. xbar = ybar = 0),
we can do the following.

Dr. A. I. Hanna (2007)


pcalib_shapesizemetric.m

 function [S] = pcalib_shapesizemetric(x, dim)

this is S = sqrt( sum (x(j) - mu_x)^2 + (y(j) - mu_y)^2 )

Dr. A. I. Hanna (2007);


pcalib_transshape.m

 function [x] = pcalib_transshape(x, temp_centroid)

Dr. A. I. Hanna (2007)


pcalib_walkpc.m

 function pcalib_walkpc(varargin)

A simple script to walk along a specific principle component.

Inputs:
Xm - the mean shape as a column vector of the form x = [x_1, y_1, x_2, y_2, ..., x_N, y_N]';
P - a matrix whose columns represent the principle components
pcind - an scalar which denotes the princple axes to move along
stdevlimits - number of standard deviations to move [-1 1];
axis_h - the handle to the axis to show the walk

Dr. A. I. Hanna (2008)


pcalib_walkpc3d.m

 function pcalib_walkpc(Xm, P, b_model, pcind, stdevlimits, dstep, axis_h)

A simple script to walk along a specific principle component.

Inputs:
Xm - the mean shape as a column vector of the form x = [x_1, y_1, x_2, y_2, ..., x_N, y_N]';
P - a matrix whose columns represent the principle components
pcind - an scalar which denotes the princple axes to move along
stdevlimits - number of standard deviations to move [-1 1];
axis_h - the handle to the axis to show the walk

Dr. A. I. Hanna (2006)


## C:\Users\AB\DArT_Toolshed\Attachments\SHLib

shlib_B.m

 function B = shlib_B(nl, theta, phi)

Constructs the matrix of basis functions
where b(i,l^2 + l + m) = Yml(theta, phi)

Dr. A. I. Hanna (2006).


shlib_Yml.m

 function Yml = shlib_Yml(l, m, theta, phi)

A matlab function that takes a given order and degree, and the matrix of
theta and phi and constructs a spherical harmonic from these. The
analogue in the 1D case would be to give a particular frequency.

Inputs:
m - order of the spherical harmonic
l - degree of the spherical harmonic
theta - matrix of polar coordinates \theta \in [0, \pi]
phi - matrix of azimuthal cooridinates \phi \in [0, 2\pi)

Example:

[x, y] = meshgrid(-1:.1:1);
z = x.^2 + y.^2;
[phi,theta,R] = cart2sph(x,y,z);
Yml = spharm_aih(2,2, theta(:), phi(:));

Dr. A. I. Hanna (2006)


shlib_decomp.m

 function [C, xs, ys, zs] = spharm_decomp(x, y, z, theta, phi, nl)

This function takes a surface defined by x, y, z together with a maximum
order of decomposition. The function returns the spherical harmonic
coefficients, C, together with the approximated surface given by xs, ys,
zs.
Example:

%Make some test data
[x, y, z] = gen_spherical_function(6, 1, [30 30]);
%Decompose into spherical harmonics
[C, xs, ys, zs] = spharm_decomp(x, y, z, 10);
%Draw the result (this is where all the time is spent)
figure(1); clf;
set(gcf, 'Color', [0 0 0],'InvertHardCopy','off');
subplot(1,2,1);
%Plot the surface
surf(x,y,z)
light; lighting phong
axis square tight; view(40,30);
title('Original Surface', 'Color', 'w');
set(gca, 'XColor', 'w', 'YColor', 'w', 'ZColor', 'w', 'Color', [0 0 0]);
%plot the estimated surface
subplot(1,2,2);
surf(xs, ys, zs);
light; lighting phong
axis square tight;
title('Estimated Surface', 'Color', 'w');
view(40,30)
set(gca, 'XColor', 'w', 'YColor', 'w', 'ZColor', 'w', 'Color', [0 0 0]);

Dr. A. I. Hanna (2006).


shlib_demo_001.m

 This is a really simple script that takes a known surface, together with
its azimuth and altitude angles and decomposes it in to a set of
spherical harmonic functions.

Dr. A. I. Hanna (2006).


shlib_demo_002.m

 This is a really simple script that takes a known surface, together with
its azimuth and altitude angles and decomposes it in to a set of
spherical harmonic functions.


shlib_demo_003.m

function shlib_demo_003
This function demonstrates the rotation invariance properties of the
spherical harmonic decomposition.

Dr. A. I. Hanna (2006).


shlib_gen_shfnc.m

 function [x, y, z, theta, phi] = shlib_gen_shfnc(degree, order, delta)

This function constructs a specific surface from spherical harmonics
with a given degree and order.

Example:

[x, y, z] = gen_spherical_function(6, 1, [10 10]);
surf(x, y, z); axis tight vis3d;

Dr. A. I. Hanna (2006)


shlib_recon.m

 [x, y, z] = meshgrid(limits(1, 1):dstep(1):limits(1, 2), limits(2, 1):dstep(2):limits(2, 2), limits(3, 1):dstep(3):limits(3, 2));
[phi,theta,R] = cart2sph(x,y,z);
theta = theta + pi/2;
phi = phi + pi;
[ROWS, COLS] = size(x);
T = theta(:);
P = phi(:);


shlib_showHarmonic.m

 function shlib_showHarmonic(degree, order)

Dr. A. I. Hanna (2007)


shlib_sph_recon.m

shlib_werner_boys_surface.m

## C:\Users\AB\DArT_Toolshed\Attachments\TPSLib\MatlabTPS

tpslib_K.m

tpslib_L.m

tpslib_P.m

tpslib_U_rbs.m

 function [U_r] = tpslib_U_rbs(r)

This function is a radial basis function used in thin-plate spline
parameter estimation.

Dr. A. I. Hanna (2005)


tpslib_matdistance.m

 DISTANCE - computes Euclidean distance matrix

E = distance(A,B)

A - (DxM) matrix
B - (DxN) matrix

Returns:
E - (MxN) Euclidean distances between vectors in A and B


tpslib_psi_tps.m

 function [tps] = tpslib_psi_tps(u, a, w, opts)

This function is the approximate warping function given by the thin plate
spline method.

Input Params:
u - the set of points you wish to warp
a, w - the parameters from pts2TPS_param
opts - the set of points that you are warping to

Output Params:
tps - a Mx2 matrix of points defined by the warping function.

Dr. A. I. Hanna (2006).


tpslib_pts2TPS_param.m

 function [w, a, K] =tpslib_ pts2TPS_param(ipts, opts)

This function calculates the parameters needed in the thin plate spline
method.

Input params:
ipts - the set of base points for an image;
opts - the set of points you want to warp an image to

Output Params:
w, a - parameters used by psi_tps
K - the Euclidean distance matrix used by the radial basis function.

Dr A. I. Hanna (2006).


tpslib_ptswarp.m

 function [wpts] = tpslib_ptswarp(pts, ipts, opts, c)

This function implements the THIN PLATE SPLINE method of point warping.

Input Params:
pts - the set of points to be warped
ipts - the set of input pts that the image is to be warped to
opts - the set of input pts that the image is already registered to

Dr. A. I. Hanna (2006).


tpslib_tpsimwarp.m

 function [wimg, w, a] = tpslib_tpsimwarp(img, ipts, opts)

This function implements the THIN PLATE SPLINE method of image warping.

Input Params:
img - the image to be warped
ipts - the set of input pts that the image is to be warped to
opts - the set of input pts that the image is already registered to

Dr. A. I. Hanna (2006).


tpswarp_c.m

gen_corr_data.m

 function [Y, E] = gen_corr_data(R, N)

So we want a transform so that Y = Q*X where E[X*X'] = I and E[Y*Y'] = R
where R is some predefined correlation matrix. With this reasoning we now
have the following.

Constraints: E[Y*Y'] = R, E[X*X'] = I, Q = Q'

Proof:

E[Y*Y'] = E[Q*X*X'*Q'] = Q*E[X*X']*Q' = Q*Q' = Q^2 = R => Q = R^(1/2)

we can now write R = V*L*V' where the columns of V are the eigenvalues of
R, and L is a diagonal matrix of correspoing eigenvalues. Therefore, we
can now write R^(1/2) = V*L^(1/2)*V', if we get some complex values then we
should just take the real part.

Inputs:
R - the desired covariance matrix
N - the number of desired samples

Ouptuts:
Y - the output data
E - mean square error of actual estimated covariance matrix and
desired covariance matrix.

Example:
R = [1 .8; 0.8 1];
[Y, E] = gen_corr_data(R, 100);
fprintf('MSE: %f\n', E);
figure;
plot(Y(:,1), Y(:,2), '.');
Dr. A. I. Hanna (2005).


harris_corner_detector.m

 function varargout = harris_corner_detector(varargin)

A method that calculates the intensity image for the harris corner
detector.

Default parameters:

'image' - the input image (default = 'pout.tif')
'method' {'harris', 'aih'} (default = 'aih')
'sigma' - the square root of the variance of the gaussian window (default = 4)
(note: if window_type='uniform', sigma controls the size of the window)
'window_type' - {'gaussian', 'uniform'} (default = 'gaussian')

Example:
corner_im = harris_corner_detector;

Example:
corner_im = harris_corner_detector('image', I);

Dr. A. I. Hanna (2007)


mk_2D_lattice.m

 MK_2D_LATTICE Return adjacency matrix for nearest neighbor connected 2D lattice
G = mk_2D_lattice(nrows, ncols, con)
G(k1, k2) = 1 iff k1=(i1,j1) is a neighbor of k2=(i2,j2)
(Two pixels are neighbors if their Euclidean distance is less than r.)
We assume no wrap around.

This is the neighborhood as a function of con:

con=4,r=1  con=8,r=sqrt(2)   con=12,r=2   con=24,r=sqrt(8)
x         x x x x x
x          x x x            x x x       x x x x x
x o x        x o x          x x o x x     x x o x x
x          x x x            x x x       x x x x x
x         x x x x x

Examples:
Consider a 3x4 grid
1  4  7  10
2  5  8  11
3  6  9  12

4-connected:
G=mk_2D_lattice(3,4,4);
find(G(1,:)) = [2 4]
find(G(5,:)) = [2 4 6 8]

8-connected:
G=mk_2D_lattice(3,4,8);
find(G(1,:)) = [2 4 5]
find(G(5,:)) = [1 2 3 4 6 7 8 9]


tslib_pts2gtfield.m

function [Gfull, Gasym, X] = tslib_pts2gtfield(varargin)

A matlab function that takes two sets of points, then it estimates the velocity
field using the motion coherence theory (MCT) for a range of X and Y values
with specific dx and dy spacing. Finally, it estimates the growth tensors
at those regularly spaced intervals by taking the derivative of the
velocity field.

Inputs:
'materialPts1' - the initial set of measured landmarks
'materialPts2' - the set of measured landmarks after the velocity field
has been applied
'sigma' - the size of the Gaussian used in the interation of the MCT (default = 1)
'lambda' - controls the smoothing of the interpolation (default = 0.01)
'xlims' - a 1x2 vector giving the limits of the velocity field in the
x-direction (default = [-1 1])
'ylims' - a 1x2 vector giving the limits of the velocity field in the
y-direction (default = [-1 1])
'dx' - the spacing between the x components (default = 0.5)
'dy' - the spacing between the y components (default = 0.5)
'verbose' - plots the interpolted field + growth tensors (default = 0)

Outputs:
'Gsym' - a 2x2xM matrix of symmetric growth tensors
'Gasym' - a 2x2xM matrix of asymmetric growth tensors (rigid body)
'X' - a 2xM matrix of (x, y) points associated with each growth tensor.

Example 1:

[Gsym, Gasym, X] = tslib_pts2gtfield('demo');

Example 2:

matPts1 = randn(20, 2);
matPts2 = matPts1.*(ones(size(matPts1,1),1)*(ones(1,2) + rand(1,2)/10 + 0.01));
[Gsym, Gasym, X] = tslib_pts2gtfield('materialPts1', matPts1, 'materialPts2', matPts2, 'verbose', 1);

(Part of the ToolShed Library)

Dr. A. I. Hanna, CMP, 2008


tslib_vel2gtfield.m

function [Gsym, Gasym] = tslib_vel2gtfield(varargin)

A Matlab function that takes a velocity matrix for the x-axis, and a
velocity matrix for the y-axis and estimates growth tensors for each
point in the velocity matrices.

Inputs:
'Vx' - the NxM velocity matrix for the x-axis
'Vy' - the NxM velocity matrix for the y-axis
'dx' - the spacing used along the x-axis
'dy' - the spacing used along the y-axis

Outputs:
'Gsym' - the symmetric growth tensor such that Vxy = Vyx,
i.e Gsym = (G + G')/2
'Gasym' - the asymmetric growth tensor such that Vxy = -Vyx and Vxx = Vyy = 0
i.e Gasym = (G - G')/2, corresponds to rigid body rotation

(Part of the ToolShed Library)

Dr. A. I. Hanna, CMP, 2008


unirndrange.m

 function r = unirndrange(x1, x2, N, M)

A Simple script that takes the range input and size of the desired
matrix and returns that size matrix whose entries are uniformly
distributed between x1 and x2.

Inputs:
x1 - the start range
x2 - the end range
N - the number of rows
M - the number of cols
Outputs:
r - the NxM matrix with uniform entries

Example 1:
r = unirndrange(10, 100, 10, 2);

Example 2:
r = unirndrange(10, 100, 100);

Example 3:
r = unirndrange(10, 100, 1, 200);

Dr. A. I. Hanna (2006).