Example a_backup/optimisation/mesh-alignment-c: Mesh alignment with the Cmgui C++ API.

This example makes use of the Cmgui C++ API to perform the same mesh alignment task achieved by the standard Cmgui mesh alignment example.


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 cmgui meshing_api 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.
 *
 * Contributor(s):
 *
 * 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 ***** */

extern "C"
{
	#include "api/cmiss_context.h"
	#include "api/cmiss_field.h"
	#include "api/cmiss_field_arithmetic_operators.h"
	#include "api/cmiss_field_composite.h"
	#include "api/cmiss_field_finite_element.h"
	#include "api/cmiss_field_group.h"
	#include "api/cmiss_field_module.h"
	#include "api/cmiss_field_matrix_operators.h"
	#include "api/cmiss_field_nodeset_operators.h"
	#include "api/cmiss_field_subobject_group.h"
	#include "api/cmiss_field_trigonometry.h"
	#include "api/cmiss_field_vector_operators.h"
	#include "api/cmiss_element.h"
	#include "api/cmiss_graphic.h"
	#include "api/cmiss_graphics_module.h"
	#include "api/cmiss_node.h"
	#include "api/cmiss_optimisation.h"
	#include "api/cmiss_region.h"
	#include "api/cmiss_rendition.h"
	#include "api/cmiss_status.h"
	#include "api/cmiss_tessellation.h"
}
#include <iostream>
#include <vector>

static Cmiss_field_id rotateCoordinates(Cmiss_field_id coordinates, Cmiss_field_id rigidBodyMotion);
static Cmiss_field_id translateCoordinates(Cmiss_field_id coordinates, Cmiss_field_id rigidBodyMotion);
static Cmiss_field_id createObjectiveField(Cmiss_field_id coordinates, Cmiss_field_id rigidBodyMotion);
//static Cmiss_field_id createDifferenceField(Cmiss_field_id field1, Cmiss_field_id field2);
static int evaluateConstantField(Cmiss_field_module_id fieldModule, Cmiss_field_id field, int numberOfValues, double* values);

int main()
{
	int return_code = 0;
	Cmiss_context_id context = Cmiss_context_create("example");

	if (context && (CMISS_OK == Cmiss_context_enable_user_interface(context, NULL)))
	{
		// get the root region
		Cmiss_region_id rootRegion = Cmiss_context_get_default_region(context);
		// load in the mesh and data points
		if ((CMISS_OK == Cmiss_region_read_file(rootRegion, "femur.exnode"))
				&& (CMISS_OK == Cmiss_region_read_file(rootRegion, "femur.exelem"))
				&& (CMISS_OK == Cmiss_region_read_file(rootRegion, "landmarks.exnode")))
		{
			std::cout << "Mesh loaded!" << std::endl;
			Cmiss_field_module_id fieldModule = Cmiss_region_get_field_module(rootRegion);
			// the coordinates of the mesh
			Cmiss_field_id coordinates = Cmiss_field_module_find_field_by_name(fieldModule, "coordinates");
			/*
			 * Define a field that will be used as the independent field for the
			 * optimisation. 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;
			 */
			const double zero[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
			Cmiss_field_id rigidBodyMotion = Cmiss_field_module_create_constant(fieldModule,6,zero);

			// create the field to be used as the objective function
			Cmiss_field_id objectiveField = createObjectiveField(coordinates, rigidBodyMotion);
			double value;
			if (evaluateConstantField(fieldModule, objectiveField, 1, &value))
			{
				std::cout << "Initial objective field value = " << value << std::endl;
			}
			else
			{
				std::cerr << "Oops, objective function can't be evaluated?" << std::endl;
			}

#if 0
			if ((CMISS_OK == Cmiss_context_execute_command(context,
					"gfx write nodes group femur_left final_mesh.exnode")))
			{
				std::cout << "Exported rotated coordinates field?" << std::endl;
			}
			else
			{
				std::cout << "Doh!" << std::endl;
			}
#endif

			// Create and set up the optimisation object
			Cmiss_optimisation_id optimisation = Cmiss_field_module_create_optimisation(fieldModule);
			Cmiss_optimisation_set_method(optimisation, CMISS_OPTIMISATION_METHOD_QUASI_NEWTON);
			Cmiss_optimisation_set_attribute_integer(optimisation, CMISS_OPTIMISATION_ATTRIBUTE_MAXIMUM_ITERATIONS, 100);
			Cmiss_optimisation_set_attribute_integer(optimisation,
					CMISS_OPTIMISATION_ATTRIBUTE_MAXIMUM_FUNCTION_EVALUATIONS, 1000);
			Cmiss_optimisation_set_attribute_real(optimisation, CMISS_OPTIMISATION_ATTRIBUTE_FUNCTION_TOLERANCE, 1.0e-8);
			Cmiss_optimisation_add_objective_field(optimisation, objectiveField);
			Cmiss_optimisation_add_independent_field(optimisation, rigidBodyMotion);

			// perform the optimisation
			Cmiss_optimisation_optimise(optimisation);

			if (evaluateConstantField(fieldModule, objectiveField, 1, &value))
			{
				std::cout << "Optimised objective field value = " << value << std::endl;
			}
			else
			{
				std::cerr << "Oops, objective function can't be evaluated?" << std::endl;
			}

			// clean up
			Cmiss_field_destroy(&objectiveField);
			Cmiss_field_destroy(&rigidBodyMotion);
			Cmiss_field_destroy(&coordinates);
			Cmiss_optimisation_destroy(&optimisation);
			Cmiss_field_module_destroy(&fieldModule);
		}
		Cmiss_region_destroy(&rootRegion);
		Cmiss_context_destroy(&context);
	}

	return return_code;
}

static Cmiss_field_id createObjectiveField(Cmiss_field_id coordinates, Cmiss_field_id rigidBodyMotion)
{
	Cmiss_field_module_id fieldModule = Cmiss_field_get_field_module(coordinates);
	// define the rotated coordinates field for the input mesh
	Cmiss_field_id rotatedCoordinates = rotateCoordinates(coordinates, rigidBodyMotion);
	Cmiss_field_set_name(rotatedCoordinates, "rotatedCoordinates");
	// and translate the rotated coordinates
	Cmiss_field_id translatedCoordinates = translateCoordinates(rotatedCoordinates, rigidBodyMotion);
	Cmiss_field_set_name(translatedCoordinates, "translatedCoordinates");

	// compute the difference between the translated coordinates and the target landmark locations
	Cmiss_field_id landmarkCoordinates = Cmiss_field_module_find_field_by_name(fieldModule, "landmark_data_location");
	Cmiss_field_id landmarkDiff = Cmiss_field_module_create_subtract(fieldModule, translatedCoordinates, landmarkCoordinates);
	Cmiss_field_id error = Cmiss_field_module_create_magnitude(fieldModule, landmarkDiff);

	if (error == NULL) std::cout << "Doh!! big time..." << std::endl;

	// create a group field with the landmark nodes only
	Cmiss_nodeset_id master_nodeset = Cmiss_field_module_find_nodeset_by_name(fieldModule, "cmiss_nodes");
	Cmiss_field_id landmarkGroupField = Cmiss_field_module_create_node_group(fieldModule, master_nodeset);
	Cmiss_field_node_group_id landmarkNodes = Cmiss_field_cast_node_group(landmarkGroupField);
	Cmiss_nodeset_group_id landmark_nodeset_group = Cmiss_field_node_group_get_nodeset(landmarkNodes);
	Cmiss_field_node_group_destroy(&landmarkNodes);
	Cmiss_field_destroy(&landmarkGroupField);
	int landmarkNodeIDs[] = {3160027, 3160053, 3160063, 3160212, 3160285, 3160402};
	for (int i = 0; i < 6; i++)
	{
		Cmiss_node_id node = Cmiss_nodeset_find_node_by_identifier(master_nodeset, landmarkNodeIDs[i]);
		Cmiss_nodeset_group_add_node(landmark_nodeset_group, node);
		Cmiss_node_destroy(&node);
	}
	Cmiss_nodeset_destroy(&master_nodeset);

	Cmiss_field_id objectiveField = Cmiss_field_module_create_nodeset_sum(fieldModule, error,
		Cmiss_nodeset_group_base_cast(landmark_nodeset_group));
	Cmiss_nodeset_group_destroy(&landmark_nodeset_group);

	if (objectiveField == NULL) std::cout << "Doh!! even bigger time..." << std::endl;

	Cmiss_field_destroy(&error);

	Cmiss_field_destroy(&rotatedCoordinates);
	Cmiss_field_destroy(&translatedCoordinates);
	Cmiss_field_destroy(&landmarkCoordinates);
	Cmiss_field_destroy(&landmarkDiff);
	return objectiveField;
}


static Cmiss_field_id rotateCoordinates(Cmiss_field_id coordinates, Cmiss_field_id rigidBodyMotion)
{
	int i;
	Cmiss_field_id matrixFields[9];
	Cmiss_field_module_id fieldModule = Cmiss_field_get_field_module(coordinates);
	double mOne[] = {-1};
	double One[] = {1};
	double Zero[] = {0};
	Cmiss_field_id minusOne = Cmiss_field_module_create_constant(fieldModule, 1, mOne);
	Cmiss_field_id one = Cmiss_field_module_create_constant(fieldModule, 1, One);
	Cmiss_field_id zero = Cmiss_field_module_create_constant(fieldModule, 1, Zero);
	Cmiss_field_id thetaX = Cmiss_field_module_create_component(fieldModule, rigidBodyMotion, 0);
	Cmiss_field_id cosX = Cmiss_field_module_create_cos(fieldModule, thetaX);
	Cmiss_field_id sinX = Cmiss_field_module_create_sin(fieldModule, thetaX);
	Cmiss_field_id minusSinX = Cmiss_field_module_create_multiply(fieldModule, sinX, minusOne);
	Cmiss_field_id thetaY = Cmiss_field_module_create_component(fieldModule, rigidBodyMotion, 1);
	Cmiss_field_id cosY = Cmiss_field_module_create_cos(fieldModule, thetaY);
	Cmiss_field_id sinY = Cmiss_field_module_create_sin(fieldModule, thetaY);
	Cmiss_field_id minusSinY = Cmiss_field_module_create_multiply(fieldModule, sinY, minusOne);
	Cmiss_field_id thetaZ = Cmiss_field_module_create_component(fieldModule, rigidBodyMotion, 2);
	Cmiss_field_id cosZ = Cmiss_field_module_create_cos(fieldModule, thetaZ);
	Cmiss_field_id sinZ = Cmiss_field_module_create_sin(fieldModule, thetaZ);
	Cmiss_field_id minusSinZ = Cmiss_field_module_create_multiply(fieldModule, sinZ, minusOne);
	// create the X-axis rotation matrix
	matrixFields[0] = Cmiss_field_module_create_identity(fieldModule, one);
	matrixFields[1] = Cmiss_field_module_create_identity(fieldModule, zero);
	matrixFields[2] = Cmiss_field_module_create_identity(fieldModule, zero);
	matrixFields[3] = Cmiss_field_module_create_identity(fieldModule, zero);
	matrixFields[4] = cosX;
	matrixFields[5] = minusSinX;
	matrixFields[6] = Cmiss_field_module_create_identity(fieldModule, zero);
	matrixFields[7] = sinX;
	matrixFields[8] = Cmiss_field_module_create_identity(fieldModule, cosX);
	Cmiss_field_id rotationMatrixX = Cmiss_field_module_create_concatenate(fieldModule, 9, matrixFields);
	for (i=0;i<9;++i) Cmiss_field_destroy(&(matrixFields[i]));
	// create the Y-axis rotation matrix
	matrixFields[0] = cosY;
	matrixFields[1] = Cmiss_field_module_create_identity(fieldModule, zero);
	matrixFields[2] = sinY;
	matrixFields[3] = Cmiss_field_module_create_identity(fieldModule, zero);
	matrixFields[4] = Cmiss_field_module_create_identity(fieldModule, one);
	matrixFields[5] = Cmiss_field_module_create_identity(fieldModule, zero);
	matrixFields[6] = minusSinY;
	matrixFields[7] = Cmiss_field_module_create_identity(fieldModule, zero);
	matrixFields[8] = Cmiss_field_module_create_identity(fieldModule, cosY);
	Cmiss_field_id rotationMatrixY = Cmiss_field_module_create_concatenate(fieldModule, 9, matrixFields);
	for (i=0;i<9;++i) Cmiss_field_destroy(&(matrixFields[i]));
	// create the Z-axis rotation matrix
	matrixFields[0] = cosZ;
	matrixFields[1] = minusSinZ;
	matrixFields[2] = Cmiss_field_module_create_identity(fieldModule, zero);
	matrixFields[3] = sinZ;
	matrixFields[4] = Cmiss_field_module_create_identity(fieldModule, cosZ);
	matrixFields[5] = Cmiss_field_module_create_identity(fieldModule, zero);
	matrixFields[6] = Cmiss_field_module_create_identity(fieldModule, zero);
	matrixFields[7] = Cmiss_field_module_create_identity(fieldModule, zero);
	matrixFields[8] = Cmiss_field_module_create_identity(fieldModule, one);
	Cmiss_field_id rotationMatrixZ = Cmiss_field_module_create_concatenate(fieldModule, 9, matrixFields);
	for (i=0;i<9;++i) Cmiss_field_destroy(&(matrixFields[i]));

	Cmiss_field_id rot1 = Cmiss_field_module_create_matrix_multiply(fieldModule, 3,
			rotationMatrixX, rotationMatrixY);
	Cmiss_field_id rot2 = Cmiss_field_module_create_matrix_multiply(fieldModule, 3,
			rot1, rotationMatrixZ);
	Cmiss_field_id rotatedCoordinates = Cmiss_field_module_create_matrix_multiply(fieldModule, 3,
			rot2, coordinates);

	Cmiss_field_destroy(&rot1);
	Cmiss_field_destroy(&rot2);
	Cmiss_field_destroy(&rotationMatrixX);
	Cmiss_field_destroy(&rotationMatrixY);
	Cmiss_field_destroy(&rotationMatrixZ);
	Cmiss_field_destroy(&minusOne);
	Cmiss_field_destroy(&one);
	Cmiss_field_destroy(&zero);
	Cmiss_field_destroy(&thetaX);
	Cmiss_field_destroy(&thetaY);
	Cmiss_field_destroy(&thetaZ);

	return rotatedCoordinates;
}

static Cmiss_field_id translateCoordinates(Cmiss_field_id coordinates, Cmiss_field_id rigidBodyMotion)
{
	int i;
	Cmiss_field_id translations[3];
	Cmiss_field_module_id fieldModule = Cmiss_field_get_field_module(coordinates);
	translations[0] = Cmiss_field_module_create_component(fieldModule, rigidBodyMotion, 3);
	translations[1] = Cmiss_field_module_create_component(fieldModule, rigidBodyMotion, 4);
	translations[2] = Cmiss_field_module_create_component(fieldModule, rigidBodyMotion, 5);
	Cmiss_field_id translationsField = Cmiss_field_module_create_concatenate(fieldModule, 3, translations);
	Cmiss_field_id translatedCoordinates = Cmiss_field_module_create_add(fieldModule, coordinates, translationsField);
	for (i=0;i<3;++i) Cmiss_field_destroy(&(translations[i]));
	return translatedCoordinates;
}

static int evaluateConstantField(Cmiss_field_module_id fieldModule, Cmiss_field_id field, int numberOfValues, double* values)
{
	Cmiss_field_cache_id cache = Cmiss_field_module_create_cache(fieldModule);
	int returnCode = Cmiss_field_evaluate_real(field, cache, numberOfValues,
			values);
	Cmiss_field_cache_destroy(&cache);
	return returnCode;
}


Files used by this example are:

Name                  Modified     Size

mesh-alignment-c.cpp 23-Aug-2012 14k Makefile 20-Apr-2012 2.4k femur.exelem 20-Apr-2012 170k femur.exnode 20-Apr-2012 140k landmarks.exnode 20-Apr-2012 616 mesh-alignment-c.o 23-Aug-2012 9.8k

Download the entire example:

Name                                                    Modified     Size

examples_a_backup_optimisation_mesh-alignment-c.tar.gz 12-Aug-2014 124k

Testing status by version:

StatusTestedReal time (s)
x86-linux
cmgui-gtk-debugFailureSun Mar 6 00:06:35 20160
last breakFri Jan 11 00:12:00 20132
last successThu Jan 10 00:11:00 20134

Testing status by file:


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

Input last modified: Thu Aug 23 10:09:19 2012


CMISS Help / Examples / a_backup / optimisation / mesh-alignment-c