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

This example demonstrates the use of the 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 data point projections. The fitted mesh.

Screenshot of example a/optimisation/smooth_volume_fitting


The comfile run by this example is as follows:

# This Source Code Form is subject to the terms of the Mozilla Public
# License, v. 2.0. If a copy of the MPL was not distributed with this
# file, You can obtain one at http://mozilla.org/MPL/2.0/.
#
# Cmgui example: Smooth Volume Fitting
#
# Fit a volume mesh to surface data points while penalising
# excessive strains throughout the volume to keep overall shape.
#
# Note this example only runs on latest Cmgui developer releases v3.0.1
# since it relies on new field types: mesh_integral_squares, is_exterior, is_on_face.
# For Cmgui v3.0 use the smooth_volume_fitting_v300.com, which integrates
# using Gauss point generation, and lists outside faces explicitly.
#
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
# add faces that are exterior AND on xi2=0 of parent element
gfx define field is_exterior is_exterior;
gfx define field is_on_face_xi2_0 is_on_face face xi2_0;
gfx define field is_outside and fields is_exterior is_on_face_xi2_0;
gfx create egroup outside;
gfx modify egroup outside add faces conditional is_outside;
gfx list faces conditional outside;
#
# 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_wireframe 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 data_coordinates
# on the outside of the mesh (outside.mesh2d, the 2-D elements in group 'outside')
gfx define field found_location find_mesh_location find_nearest mesh outside.mesh2d mesh_field coordinates source_field data_coordinates;
#
# define a field for storing these locations so not recalculated during optimisation
# (Dynamically found locations could be used, but will be much slower and possibly unstable) 
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;
# assign 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;
#
if (!$TESTING) {
  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.datapoints;
#
# we aim to smooth the fit by penalising the displacement gradient in a least-squares sense
# (displacement gradient = du/dX, where u=displacement, X = 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_gradient gradient field displacement coordinate reference_coordinates;
# This example uses a relatively high alpha value (scaling displacement gradient)
# to try to keep close to original shape. In more typical fitting problems, use
# as small a value of alpha as possible which achieves a good looking fit.
# Experiment by varying alpha up and down in whole or multiple orders of
# magnitude to get close to the desired result, then fine tune.
# This example uses different scale factors for each strain component to permit
# more strain in the z direction (useful for this problem, but not in general):
gfx define field alpha constant 4 4 4 4 4 4 4 4 0.4;
gfx define field scaled_displacement_gradient multiply fields displacement_gradient alpha;
#
# form least-squares smoothing objective as integral of squares of above field
# over the mesh. Note this is integrated with the reference_coordinates field
# otherwise very non-linear!
gfx define field volume_smoothing_objective mesh_integral_squares integrand_field scaled_displacement_gradient coordinate reference_coordinates mesh mesh3d gaussian_quadrature numbers_of_points "4";
#
# if testing, run the optimisation otherwise the user can manually choose to run the optimisation
# When running interactively, just run the 'gfx minimise' command below
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 10-Jul-2014 7.7k boxpoints.exnode 17-Mar-2014 49k figure8.exelem 17-Mar-2014 33k figure8.exnode 10-Jul-2014 17k smooth_volume_fitting_v300.com 10-Jul-2014 8.9k

Download the entire example:

Name                                                  Modified     Size

examples_a_optimisation_smooth_volume_fitting.tar.gz 09-Mar-2016 466k

Testing status by version:

StatusTestedReal time (s)
i686-linux
cmgui-wxFailureSun Mar 6 00:06:41 20161
last breakMon Aug 25 13:12:00 20141
cmgui-wx-debugFailureSun Mar 6 00:05:45 20161
last breakMon Aug 25 13:12:00 20141
cmgui-wx-debug-memorycheckFailureSun Mar 6 00:05:46 20161
last breakMon Aug 25 13:12:00 20141
cmgui-wx-debug-valgrindFailureSun Mar 6 00:31:51 201618
last breakTue Feb 24 00:05:00 201522
x86_64-linux
cmgui-wxFailureSun Mar 6 00:01:32 20161
last breakFri Aug 15 00:33:00 20141
cmgui-wx-debugFailureSun Mar 6 00:01:32 20161
last breakFri Aug 15 00:33:00 20140
cmgui-wx-debug-memorycheckFailureSun Mar 6 00:01:32 20161
last breakFri Aug 15 00:33:00 20141
cmgui-wx-debug-valgrindFailureSun Mar 6 00:03:56 201616
last breakSun Mar 6 00:03:00 201616

Testing status by file:


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

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


CMISS Help / Examples / a / optimisation / smooth_volume_fitting