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
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\MCTLib\mctlibj
C:\Users\AB\DArT_Toolshed\Attachments\MCTLib\mctlibj_dist
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
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).