This example makes use of the Cmgui C++ API to perform the same mesh alignment task achieved by the standard Cmgui mesh alignment example.
/* ***** 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_constant.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, "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;
}
Name Modified Size
mesh-alignment-c.cpp 17-Mar-2014 14k Makefile 17-Mar-2014 2.4k femur.exelem 17-Mar-2014 170k femur.exnode 17-Mar-2014 140k landmarks.exnode 17-Mar-2014 616
Name Modified Size
examples_a_optimisation_mesh-alignment-c.tar.gz 09-Mar-2016 120k
| Status | Tested | Real time (s) | |
| x86-linux | |||
| cmgui-gtk-debug | Failure | Sun Mar 6 00:05:44 2016 | 0 |
| last break | Fri Aug 15 03:01:00 2014 | 0 |
| x86-linux | |||
| 2 | cmgui-gtk-debug: | error exit status 2. |
Html last generated: Wed Mar 9 16:02:28 2016
Input last modified: Wed Mar 9 15:49:41 2016