Example a_backup/optimisation/smooth_volume_fitting: Least-squares volume fitting with smoothing.

This example demonstrates the use of the non-linear least-squares optimisation method to fit outside surfaces of a three-dimensional tricubic Hermite mesh (with multiple derivative versions) to a set of data points while penalising deformation of the volume so it broadly maintains its shape. Uses the least-squares method as it is more efficient and well-suited to this type of problem.
Initial mesh showing the original mesh and the data point projections. The fitted mesh.

Screenshot of example a_backup/optimisation/smooth_volume_fitting


The comfile run by this example is as follows:

# ***** BEGIN LICENSE BLOCK *****
# Version: MPL 1.1/GPL 2.0/LGPL 2.1
#
# The contents of this file are subject to the Mozilla Public License Version
# 1.1 (the "License"); you may not use this file except in compliance with
# the License. You may obtain a copy of the License at
# http://www.mozilla.org/MPL/
#
# Software distributed under the License is distributed on an "AS IS" basis,
# WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
# for the specific language governing rights and limitations under the
# License.
#
# The Original Code is CMISS smooth volume fitting example.
#
# The Initial Developer of the Original Code is
# Auckland Uniservices Ltd, Auckland, New Zealand.
# Portions created by the Initial Developer are Copyright (C) 2011
# the Initial Developer. All Rights Reserved.
#
# Alternatively, the contents of this file may be used under the terms of
# either the GNU General Public License Version 2 or later (the "GPL"), or
# the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
# in which case the provisions of the GPL or the LGPL are applicable instead
# of those above. If you wish to allow use of your version of this file only
# under the terms of either the GPL or the LGPL, and not to allow others to
# use your version of this file under the terms of the MPL, indicate your
# decision by deleting the provisions above and replace them with the notice
# and other provisions required by the GPL or the LGPL. If you do not delete
# the provisions above, a recipient may use your version of this file under
# the terms of any one of the MPL, the GPL or the LGPL.
#
# ***** END LICENSE BLOCK *****

gfx define tessellation default minimum_divisions "1" refinement_factors "12";
# Read the tricubic Hermite 'figure 8' mesh with 'coordinates' field
gfx read node $example/figure8.exnode;
gfx read element $example/figure8.exelem;
# Rename 'coordinates' to 'reference_coordinates' to keep at initial value
gfx modify field coordinates rename reference_coordinates;
# Re-read field coordinates:
gfx read node $example/figure8.exnode;
gfx read element $example/figure8.exelem;
# create a group containing the outside faces for finding data point projections on
gfx create egroup outside;
gfx modify egroup outside add faces 3,8,13,18,23,27;
# read some random data points on a box around the mesh to fit to
gfx read data $example/boxpoints.exnode;

# visualise the figure 8 mesh and data points
gfx modify g_element "/" general clear;
gfx modify g_element "/" lines coordinate coordinates tessellation default LOCAL native_discretization NONE select_on material default selected_material default_selected;
gfx modify g_element "/" surfaces coordinate coordinates exterior tessellation default LOCAL native_discretization NONE select_on material silver selected_material default_selected render_shaded invisible;
gfx modify g_element "/" node_points coordinate coordinates LOCAL glyph sphere general size "0.1*0.1*0.1" centre 0,0,0 font default select_on material default selected_material default_selected;
gfx modify g_element outside general clear;
gfx modify g_element outside lines coordinate coordinates tessellation default LOCAL native_discretization NONE select_on material default selected_material default_selected invisible;
gfx modify g_element outside surfaces coordinate coordinates tessellation default LOCAL native_discretization NONE select_on material muscle selected_material default_selected render_shaded;
gfx modify g_element boxpoints general clear;
gfx modify g_element boxpoints data_points coordinate data_coordinates LOCAL glyph cross general size "0.1*0.1*0.1" centre 0,0,0 font default select_on material green selected_material default_selected;

# create a field which dynamically finds the nearest location to the scaled_data_coordinates on the 1-D mesh 
gfx define field found_location find_mesh_location find_nearest mesh outside.cmiss_mesh_2d mesh_field coordinates source_field data_coordinates;

# define a field for storing these locations so not recalculated during optimisation
# (You can use the dynamically found locations but it is slower and probably unstable without smoothing) 
gfx define field stored_location finite_element element_xi;
# define storage for the the stored_location field at the boxpoints:
gfx modify data group boxpoints define stored_location;
# set stored_location field from found_location at the data points
gfx evaluate dgroup boxpoints destination stored_location source found_location;

# define a field giving the coordinates at the stored_location on the mesh:
gfx define field projected_coordinates embedded element_xi stored_location field coordinates;
# define a field which calculates the error vector between the data point and its element project:
gfx define field error_vector add fields projected_coordinates data_coordinates scale_factors 1 -1;
# get the magnitude for visualisation
gfx define field error magnitude field error_vector;

# visualise error bars on the boxpoints:
gfx modify g_element boxpoints data_points coordinate data_coordinates LOCAL glyph line general size "0*0.1*0.1" centre 0,0,0 font default orientation error_vector scale_factors "1*0*0" select_on material green data error spectrum default selected_material default_selected;
gfx modify spectrum default autorange

gfx create window 1
gfx modify window 1 view perspective eye_point 7.97495 -10.9082 12.2888 interest_point 0 0 -0.00470632 up_vector -0.363096 0.567435 0.739039 view_angle 39.2727;

# define surface fit objective function as sum of squared error_vector over the data points
# (least-squares method works directly with individual terms from sum-of-squares field types like this)
gfx define field outside_surface_fit_objective nodeset_sum_squares field error_vector nodeset boxpoints.cmiss_data;

# create 4*4*4 Gauss points on the mesh to integrate smoothing terms with
gfx define field gauss_location finite_element number_of_components 1 field element_xi;
gfx define field gauss_weight finite_element number_of_components 1 field real;
gfx define field gauss_coordinates embedded element_xi gauss_location field coordinates;
gfx create dgroup gausspoints region "/";
gfx create gauss_points region "/" mesh cmiss_mesh_3d order 4 gauss_location_field gauss_location gauss_point_nodeset gausspoints.cmiss_data gauss_weight_field gauss_weight first_identifier 1001;

# visualise gauss points
gfx modify g_element gausspoints general clear;
gfx modify g_element gausspoints data_points coordinate gauss_coordinates LOCAL glyph cross general size "0.05*0.05*0.05" centre 0,0,0 font default select_on material tissue selected_material default_selected;

# obtain the differential reference coordinates change with chart coordinate 'xi'
gfx define field dX_dxi gradient field reference_coordinates coordinate xi;
gfx define field dV determinant field dX_dxi;

# we aim to smooth the fit by penalising du/dX in a least-squares sense
# (du/dX = change of displacement w.r.t. reference coordinates)
# Note this is just one possible smoothing approach; other ways include
# penalising second derivatives d2u/dX2 to limit curvature.
gfx define field displacement add fields coordinates reference_coordinates scale_factors 1 -1;
gfx define field displacement_xy composite displacement.1 displacement.2
gfx define field displacement_gradient gradient field displacement_xy coordinate reference_coordinates;

# use 'embedded' fields to obtain the above quantities at the gauss points
gfx define field gauss_displacement_gradient embedded element_xi gauss_location field displacement_gradient;
gfx define field gauss_dV embedded element_xi gauss_location field dV;
# multiply to get integration terms
gfx define field gauss_weight_dV multiply_components fields gauss_weight gauss_dV;
# 'alpha' is a scale factor for the smoothing. Try a range of values varying by whole orders of magnitude
# to achieve a solution which looks sufficiently like the original mesh but close to the data points.
gfx define field alpha constant 20;
gfx define field weight multiply_components fields gauss_weight_dV alpha;
gfx define field scaled_gauss_displacement_gradient multiply_components fields gauss_displacement_gradient weight;

# form least-squares smoothing objective as sum of squares of above field over gauss points
gfx define field volume_smoothing_objective nodeset_sum_squares field scaled_gauss_displacement_gradient nodeset gausspoints.cmiss_data;

# if testing, run the optimisation otherwise the user can manually choose to run the optimisation
if ($TESTING) {
	# Setting alpha = 0 or removing '& volume_smoothing_objective' from the following
        # Gives a solution with no change to interior DOFs. Note omitting smoothing causes
	# oscillations in the solution where few data points are present to constrain it.
	# 1 iteration should be enough because the problem is essentially linear
	gfx minimise LEAST_SQUARES_QUASI_NEWTON region "/" objective_fields outside_surface_fit_objective & volume_smoothing_objective independent_fields coordinates maximum_iterations 1;
	gfx write nodes field coordinates non_recursive fitted_figure8.exnode;
}

# Extra ideas:
# 1. Try updating the stored_location field and solve again
# gfx evaluate dgroup boxpoints destination stored_location source found_location;
# 2. Try fitting with dynamically found mesh locations (warning: much slower & non-linear so may need to iterate)
# gfx define field projected_coordinates embedded element_xi found_location field coordinates;


Files used by this example are:

Name                       Modified     Size

smooth_volume_fitting.com 20-Apr-2012 9.3k boxpoints.exnode 20-Apr-2012 49k figure8.exelem 20-Apr-2012 33k figure8.exnode 20-Apr-2012 17k

Download the entire example:

Name                                                         Modified     Size

examples_a_backup_optimisation_smooth_volume_fitting.tar.gz 12-Aug-2014 474k

Testing status by version:

StatusTestedReal time (s)
i686-linux
cmgui-wxSuccessSun Mar 6 00:36:55 201621
cmgui-wx-debugSuccessSun Mar 6 00:51:24 201634
cmgui-wx-debug-memorycheckSuccessSun Mar 6 00:55:53 201638
cmgui-wx-debug-valgrindFailureSun Mar 6 03:20:57 2016868
last breakSun Mar 6 03:06:00 2016868
last successTue Feb 10 03:27:00 2015840
x86_64-linux
cmgui-wxFailureSun Mar 6 00:01:37 20160
last breakSun Mar 6 00:01:00 20160
last successWed Jun 3 00:20:00 201523
cmgui-wx-debugFailureSun Mar 6 00:01:37 20160
last breakSun Mar 6 00:01:00 20160
last successWed Jun 3 00:25:00 201530
cmgui-wx-debug-memorycheckFailureSun Mar 6 00:01:38 20160
last breakSun Mar 6 00:01:00 20160
last successWed Jun 3 00:28:00 201532
cmgui-wx-debug-valgrindFailureSun Mar 6 00:03:13 20169
last breakSun Mar 6 00:03:00 20169
last successWed Jun 3 02:05:00 2015644

Testing status by file:


Html last generated: Sun Mar 6 05:50:59 2016

Input last modified: Fri Apr 20 16:00:35 2012


CMISS Help / Examples / a_backup / optimisation / smooth_volume_fitting