Example a/optimisation/mesh-alignment: Aligning a mesh using landmarks.

This example uses specified landmarks to align a geometric mesh with a data cloud for further fitting. In this example, a six component constant field is defined and used to define the rigid body rotation and translation of the geometric mesh. The six components are the rotation about and the translation along each of the X, Y, and Z axes.
Initial mesh showing the original landmark locations (cubes) and the corresponding locations in the data cloud (spheres). Mesh aligned using QN method with the initial mesh location shown with the lines.


The comfile run by this example is as follows:

if (! defined $example){
	$example = ".";
}

# read in the initial mesh
gfx read node "$example/femur.exnode";
gfx read element "$example/femur.exelem";
# and the data points (not used in this example apart from visualisation)
gfx read nodes "$example/data.exdata";
# and the identified landmark nodes from the mesh that we want to fit to. These should be the nodes from geometric mesh but 
# with the coordinates you want to align the mesh with  (i.e., the corresponding coordinate in the data cloud).
gfx read nodes "$example/landmarks.exnode";

$i = 1;
sub defineRotationFields {
	my $axis = shift;
	my $n1,$n2,$n3,$n4;
	$n1 = "theta" . $axis;
	gfx define field $n1 component rotations.$i;
	$i++;
	$n2 = "cos" . $axis;
    gfx define field $n2 cos field $n1;
    $n3 = "sin" . $axis;
    gfx define field $n3 sin field $n1;
    $n4 = "minusSin" . $axis;
    gfx define field $n4 scale field $n3 scale_factor -1;
}

# This is the field defining the parameters to be fit. The components of this field will be used to define:
# 1) Rotation about the X axis;
# 2) Rotation about the Y axis;
# 3) Rotation about the Z axis;
# 4) Translation along the X axis;
# 5) Translation along the Y axis; and
# 6) Translation along the Z axis;
gfx define field rigidBodyMotion constant 0 0 0 0 0 0;

gfx define field rotations component rigidBodyMotion.1 rigidBodyMotion.2 rigidBodyMotion.3;
&defineRotationFields("X");
gfx define field rotationMatrixX composite 1 0 0 0 cosX minusSinX 0 sinX cosX; 
&defineRotationFields("Y");
gfx define field rotationMatrixY composite cosY 0 sinY 0 1 0 minusSinY 0 cosY; 
&defineRotationFields("Z");
gfx define field rotationMatrixZ composite cosZ minusSinZ 0 sinZ cosZ 0 0 0 1; 

gfx define field coordinatesRotatedX matrix_multiply number_of_rows 3 fields rotationMatrixX coordinates;
gfx define field coordinatesRotatedY matrix_multiply number_of_rows 3 fields rotationMatrixY coordinates;
gfx define field coordinatesRotatedZ matrix_multiply number_of_rows 3 fields rotationMatrixZ coordinates;

gfx define field bob1 matrix_multiply number_of_rows 3 fields rotationMatrixX rotationMatrixY;
gfx define field bob2 matrix_multiply number_of_rows 3 fields bob1 rotationMatrixZ;
gfx define field rotatedCoordinates matrix_multiply number_of_rows 3 fields bob2 coordinates;

gfx define field translations component rigidBodyMotion.4 rigidBodyMotion.5 rigidBodyMotion.6;
gfx define field translatedCoordinates add fields rotatedCoordinates translations;

gfx define field landmark_diff add fields translatedCoordinates landmark_data_location scale_factors 1 -1;
gfx define field error magnitude field landmark_diff;
gfx create ngroup landmarks add_ranges 3160027,3160053,3160063,3160212,3160285,3160402;
gfx define field objective_value nodeset_sum field error nodeset landmarks.nodes;

sub updateNodalCoordinates {
	my $file = shift;
	# FIXME: this only does one thing, i.e., value. need to call multiple times to get derivatives evaluated also...
	#        need to define nodal value fields for each derivative and version and evaluate them individually and reassemble,
	#        or update the evaluate code to work properly :)
	gfx evaluate source translatedCoordinates dest coordinates ngroup femur_left;
	if (defined $file) {
		gfx write node group femur_left fields coordinates "aligned.exnode";
	}
}

gfx modify g_element femur_left general clear;
gfx modify g_element femur_left lines coordinate coordinates select_on material default selected_material default_selected;
gfx modify g_element femur_left lines coordinate translatedCoordinates select_on material gold selected_material default_selected;

# draw the data points
gfx modify g_element data_points general clear;
gfx modify g_element data_points node_points coordinate data_coordinates glyph point general size "1*1*1" centre 0,0,0 font default select_on material green selected_material default_selected;

# draw the landmark locations
gfx modify g_element landmarks general clear;
gfx modify g_element landmarks node_points coordinate translatedCoordinates glyph cross general size "25*25*25" centre 0,0,0 font default select_on material green selected_material default_selected;
gfx modify g_element landmarks node_points coordinate landmark_data_location glyph sphere_hires general size "10*10*10" centre 0,0,0 font default select_on material green selected_material default_selected;


gfx create window 1 double_buffer;
gfx modify window 1 image scene default infinite_viewer_lighting two_sided_lighting;
gfx modify window 1 image add_light default;
gfx modify window 1 image add_light default_ambient;
gfx modify window 1 layout simple ortho_axes z -y eye_spacing 0.25 width 512 height 512;
gfx modify window 1 set current_pane 1;
gfx modify window 1 background colour 0 0 0 texture none;
gfx modify window 1 view parallel eye_point -346.262 -1150.62 74.6099 interest_point 90.5635 127.151 28.2347 up_vector 0.0531056 0.0180817 0.998425 view_angle 30.2232 near_clipping_plane 13.5118 far_clipping_plane 4828.64 relative_viewport ndc_placement -1 1 2 2 viewport_coordinates 0 0 1 1;
gfx modify window 1 set transform_tool current_pane 1 std_view_angle 40 normal_lines no_antialias depth_of_field 0.0 fast_transparency blend_normal;

if ($TESTING) {

   # use the quasi-Newton Opt++ method to minimise the objective function field.
   # [this method minimises a scalar valued field by altering the provided independent field(s)]
   # and demonstrates the use of the 'hide_output' option on the minimise command.
   gfx minimise QUASI_NEWTON region "/" objective_fields objective_value independent_fields rigidBodyMotion hide_output;
   
   # list out the optimised field for testing
   gfx list field rigidBodyMotion;

   # reset and use least squares minmisation
   gfx define field rigidBodyMotion constant 0 0 0 0 0 0;
   
	# use sum of squared diffs for least squares fit
	gfx define field objective_value_lsq nodeset_sum_squares field landmark_diff nodeset landmarks.nodes;
	
	# fit the data using quasi-Newton least-squares
	gfx minimise LEAST_SQUARES_QUASI_NEWTON region "/" objective_fields objective_value_lsq independent_fields rigidBodyMotion hide_output;

   # list out the optimised field for testing
   gfx list field rigidBodyMotion;

}


Files used by this example are:

Name                Modified     Size

mesh-alignment.com 09-Mar-2016 6.1k data.exdata 17-Mar-2014 71k femur.exelem 17-Mar-2014 170k femur.exnode 17-Mar-2014 140k landmarks.exnode 17-Mar-2014 616

Download the entire example:

Name                                           Modified     Size

examples_a_optimisation_mesh-alignment.tar.gz 09-Mar-2016 267k

Testing status by version:

StatusTestedReal time (s)
i686-linux
cmgui-wxFailureSun Mar 6 00:10:51 20162
last breakTue Feb 24 03:16:00 20152
cmgui-wx-debugFailureSun Mar 6 00:11:16 20164
last breakTue Feb 24 03:15:00 20152
cmgui-wx-debug-memorycheckFailureSun Mar 6 00:11:58 20164
last breakTue Feb 24 03:16:00 20152
cmgui-wx-debug-valgrindFailureSun Mar 6 01:10:50 201646
last breakSun Mar 6 01:10:00 201646
x86_64-linux
cmgui-wxFailureSun Mar 6 00:01:31 20160
last breakSun Mar 6 00:01:00 20160
cmgui-wx-debugFailureSun Mar 6 00:01:32 20160
last breakSun Mar 6 00:01:00 20160
cmgui-wx-debug-memorycheckFailureSun Mar 6 00:01:32 20161
last breakSun Mar 6 00:01:00 20161
cmgui-wx-debug-valgrindFailureSun Mar 6 00:02:56 20169
last breakSun Mar 6 00:02:00 20169

Testing status by file:


Html last generated: Wed Mar 9 16:02:27 2016

Input last modified: Wed Mar 9 15:49:41 2016


CMISS Help / Examples / a / optimisation / mesh-alignment