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).

femlib_NGrad.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);
          _                                   _
 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).

femlib_addRandomZ.m

 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_setRadialTemp.m

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)

 See also: gtlib_interp_tri_pts
 
 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

 See also: gtlib_interpTriangle

 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)


 See also: tpslib_ptswarp

 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

 See also: gtlib_plotGrowthTensorCross

 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

 See also: gtlib_plotGrowthTensor

 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

 See also: gtlib_plotGrowthTensor

 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:

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

 See also: gtlib_calcGrowthParameters

 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)
 
 See also: mctlib_CalcVel, mctlib_GaussGram, mctlib_normalise

 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
 
 See also: mctlib_CalcParam, mctlib_GaussGram, mctlib_normalise

 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))
 
 See also: mctlib_CalcParam, mctlibCalcVel, mctlib_normalise

 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

 See also: mctlib_CalcParam, mctlibCalcVel, mctlib_GaussGram

 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

 See also: mctlib_CalcParam, mctlibCalcVel, mctlib_GaussGram

 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)

 See also: shlib_Yml.m, shlib_decomp.m, shlib_gen_shfnc.m

 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(:));

 See also: shlib_B.m, shlib_decomp.m, shlib_gen_shfnc.m

 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]);

 See also: B_spharm.m, spharm_aih.m, gen_spherical_function.m

 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; 

 See also: shlib_B.m, shlib_Yml.m, shlib_decomp.m

 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.
 
 see also: find_tps_params

 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:
 I = imread('circuit.tif');
 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).