![[Home]](/images/iso2mesh_2015_banner.png)
A registration toolbox for FEM surface meshes
Author: Qianqian Fang <fangq at nmr.mgh.harvard.edu>
- Photon Migration Lab
- Martinos Center for Biomedical Imaging
- Massachusetts General Hospital (Harvard Medical School)
- Bldg. 149, 13th St., Charlestown, MA 02148
Table of Content
- 1. Introduction
- 2. List of functions
- 3. Acknowledgement
1. Introduction
"Metch", coined from "mesh" and "match", is a Matlab/Octave-based mesh/volume registration toolbox. It provides straightforward functions to register point clouds (or surfaces) to a triangular/cubic surface mesh by calculating an optimal affine transformation (in terms of matrix A for scaling and rotation, and b for translation). It also allows one to project a point cloud onto the surface using surface norms and guarantee the conformity of the points to the surface.
2. List of functions
metchgui.m
alldata = metchgui(node,elem,points) or metchgui(volume,points)
A GUI to register a point cloud to a mesh or volumetric image
parameters:
node: node coordinate of the surface mesh (nn x 3)
elem: element list of the surface mesh (3 columns for
triangular mesh, 4 columns for cubic surface mesh)
the input can also be two parameters in form of metchgui(volume,points),
where volume is a 3D image (array).
outputs:
alldata: a structrure containing all processing outputs
the fields include:
.node: the input node
.elem: the input surface mesh elements
.volume: if the input volumetric image
.A0: the affine rotation for selected point pairs (after Initialize)
.b0: the affine translation for selected point pairs (after Initialize)
.A: the affine rotation for the point cloud (after Optimize)
.b: the affine translation for the point cloud (after Optimize)
.points: the input point cloud
.pointsinit: the point cloud after initialization
.pointsopt: the point cloud after optimization
.pointsproj: the point cloud after projecting to the surface
.initplot: the handle to the point cloud plot after init
.optplot: the handle to the point cloud plot after optimization
.projplot: the handle to the point cloud plot after projection
If user supplys an output variable, the GUI will not return until the user hits the "close" button or close the window; if user does not supply any output, the call will return immediately; any data user intends to save, he has to click on "Save Session" button and provides a mat-file file name. A single structure named "metchsession" will be stored in this file.
This function is matlab-only.
metchgui_one.m
alldata = metchgui_one(node,elem,points) or metchgui_one(volume,points)
this is the same as metchgui.m except it does not need metchgui.fig in order
to run. It can also be run on lower versions of matlab.
regpt2surf.m
[A,b,newpos]=regpt2surf(node,elem,p,pmask,A0,b0,cmask,maxiter)Perform point cloud registration to a triangular surface (surface can be either triangular or cubic), Gauss-Newton method is used for the calculation
parameters:
node: node coordinate of the surface mesh (nn x 3)
elem: element list of the surface mesh (3 columns for
triangular mesh, 4 columns for cubic surface mesh)
p: points to be registered, 3 columns for x,y and z respectively
pmask: a mask vector with the same length as p, determines the
method to handle the point, if pmask(i)=-1, the point is a free
node and can be move by the optimization, if pmask(i)=0, the
point is fixed; if pmask(i)=n>0, the distance between p(i,:)
and node(n,:) will be part of the object function and be optimized
A0: a 3x3 matrix, as the initial guess for the affine A matrix (rotation&scaling)
b0: a 3x1 vector, as the initial guess for the affine b vector (translation)
cmask: a binary 12x1 vector, determines which element of [A(:);b] will be optimized
maxiter: a integer, specifying the optimization iterations
outputs:
newpos: the registered positions for p, newpos=A*p'+b
proj2mesh.m
[newpt elemid weight]=proj2mesh(v,f,pt,nv,cn)
project a point cloud on to the surface mesh (surface can only be triangular)
parameters:
v: node coordinate of the surface mesh (nn x 3)
f: element list of the surface mesh (3 columns for
triangular mesh, 4 columns for cubic surface mesh)
pt: points to be projected, 3 columns for x,y and z respectively
nv: nodal norms (vector) calculated from nodesurfnorm.m
with dimensions of (size(v,1),3)
cn: a integer vector with the length of p, denoting the closest
surface nodes (indices of v) for each point in p. this
value can be calculated from dist2surf.m
if nv and cn are not supplied, proj2mesh will project the point
cloud onto the surface by the direction pointing to the centroid
of the mesh
outputs:
newpt: the projected points from p
elemid: a vector of length of p, denotes which surface trangle (in elem)
contains the projected point
weight: the barycentric coordinate for each projected points, these are
the weights
affinemap.m
[A,b]=affinemap(pfrom,pto)
calculate an affine transform (A matrix and b vector) to map n vertices from one space to the other using least square solutions
parameters:
pfrom: nx3 matrix, each row is a 3d point in original space
pto: nx3 matrix, each row is a 3d point in the mapped space
outputs:
A: 3x3 matrix, the calculated affine A matrix
b: 3x1 vector, the calculated affine b vector
the solution will satisfy the following equation: A*pfrom'+b=pto
dist2surf.m
[d2surf,cn]=dist2surf(node,nv,p)
calculate the distances from a point cloud to a surface, and return the indices of the closest surface node
parameters:
node: node coordinate of the surface mesh (nn x 3)
nv: nodal norms (vector) calculated from nodesurfnorm.m
with dimensions of (size(node,1),3), this can be
calcuated from nodesurfnorm.m
pt: points to be calculated, 3 columns for x,y and z respectively
outputs:
d2surf: a vector of length of p, the distances from p(i) to the surface
cn: a integer vector with the length of p, the indices of the closest surface node
getplanefrom3pt.m
[a,b,c,d]=getplanefrom3pt(plane)
calculate the plane equation coefficients for a plane (determined by 3 points), the plane equation is a*x+b*y+c*z+d=0
parameters:
plane: a 3x3 matrix, each row is a 3d point in form of (x,y,z)
this is used to define a plane
outputs:
a,b,c,d: the coefficients of the plane equation
linextriangle.m
[isinside,pt,coord]=linextriangle(p0,p1,plane)
calculate the intersection of a 3d line (passing two points) with a plane (determined by 3 points)
parameters:
p0: a 3d point in form of (x,y,z)
p1: another 3d point in form of (x,y,z), p0 and p1 determins the line
plane: a 3x3 matrix, each row is a 3d point in form of (x,y,z)
this is used to define a plane
outputs:
isinside: a boolean variable, 1 for the intersection is within the
3d triangle determined by the 3 points in plane; 0 is outside
pt: the coordinates of the intersection pint
coord: 1x3 vector, if isinside=1, coord will record the barycentric
coordinate for the intersection point within the triangle;
otherwise it will be all zeros.
nodesurfnorm.m
nv=nodesurfnorm(node,elem)
calculate a nodal norm for each vertix on a surface mesh (surface can only be triangular or cubic)
parameters:
node: node coordinate of the surface mesh (nn x 3)
elem: element list of the surface mesh (3 columns for
triangular mesh, 4 columns for cubic surface mesh)
pt: points to be projected, 3 columns for x,y and z respectively
outputs:
nv: nodal norms (vector) calculated from nodesurfnorm.m
with dimensions of (size(v,1),3)
trisurfnorm.m
ev=trisurfnorm(node,elem)calculate the surface norms for each element (surface can be either triangular or cubic)
parameters:
node: node coordinate of the surface mesh (nn x 3)
elem: element list of the surface mesh (3 columns for
triangular mesh, 4 columns for cubic surface mesh)
outputs:
ev: norm vector for each surface element
3. Acknowledgement
This toolbox was developed with the support from NIH grant titled "Dynamic Inverse Solutions for Multimodal Imaging" (R01EB006385)