FE geometric face-fitting to laser scanned bone data

Example 21i: FE geometric face-fitting to laser scanned bone data

This example covers:

creation of an initial mesh comprised of linear Lagrange elements using nodes obtained using our Faro Arm.
acquisition of data points on the surface of a cast model of an irregular human bone (the sternum) using the Faro-mounted Model Maker 4 Laser Scanner.
fitting of the initial mesh to this surface data using face-fitting for volume elements.

This example links to the following examples which also contain relevant information:

Example 21d: FE geometric fitting to laser scenned bone data covers creation of an initial mesh comprised of linear Lagrange elements using nodes obtained using our Faro Arm.
Example 21g: Face fitting(geometric) in volume elements is the original example on face-fitting from the person who wrote most of the code (Kumar Mithraratne).

The process described below is that used by the NERF-funded musculoskeletal modelling project currently being undertaken in the Bioengineering Institute for production of volume meshes for the irregular bones of the human body.

Created by: Steven Thrupp 27 Aug 2002

Creating an initial mesh: using the Faro Arm
The initial linear Lagrange mesh used in this process was created using our Faro Arm. The process of using the Faro Arm in constructing an initial mesh is described in example 21d.

Obtaining Data for Fitting: Model Maker 4
While use of the HLS (Hand Held Laser Scanner) is described in example 21d, this example gives a brief description of our other laser scanner, Model Maker 4 (and a pointer to more extensive documentation for its use and calibration). Model Maker 4 is a Faro Arm-mounted laser scanner. As such it has a number of advantages over the HLS:

It uses the Faro-Arm to provide its coordinate system therefore is not subject to inaccuracies related to interference with the magnetic tracking system which were a problem with the HLS.
Has a higher degree of accuracy associated with it.
Able to change the sampling rate in the x and y direction when scanning.
You have the ability to interact more directly with the resulting surface mesh therefore allowing repair or smoothing of mesh in certain regions and also limiting the area where data may be acquired when scanning (thereby avoiding overlap of data).
Perhaps most importantly the Faro-mounted laser scanner enables the user to obtain reference points which can subsequently be used to register scans of an object that has been moved between scans. This allows the user to move an object into an optimal position for scanning it and then move it later when another position is required to adequately scan the surface geometry.

Figure 1: Faro-mounted laser scanner
In order to use the Laser scanner, or if using it for the first time after a period of disuse, there are two processes the user should go through first. First, you should ensure that the Faro Arm itself is calibrated and second, ensure that the laser scanner is correctly aligned with the Faro-Arm. This is a slightly time-consuming process however it is essential to ensure that these two things are correct otherwise one will run into difficulty down the line when trying to align your data within the Model Maker package.

Faro Arm Calibration
Calibrating the Faro arm is straightforward. It is described in example 21d however is described here again on the Windows 2000 OS. Program Menu >> Faro >> Caliper 3D. Within Caliper 3D go Settings >> Probes. Within This dialogue boss select either Point Probe or Ball Probe. If Ball probe you will need to specify the diameter. (If point probe no diameter is necessary.) Then look in the Calibration box within the current dialogue box and go 1"Ball. Depressing this button will provide you with a new dialogue box. From this point on you are required to obtain 27 points by clicking on the 1"Ball magnetically connected to the base of the Faro Arm. The points should be obtained in a reasonably uniform fashion (ie, 5 points around the horizontal equator of the ball, then 5 around the vertical, followed by e.g. around the top of the ball etc.). Points are obtained by clicking a front button on the Faro Arm when in position and then clicking a back button in order to accept the the point you have just obtained. The software will count off the number of points you have obtained. When you have obtained these 27 points the software will calculate an RMS for you which will indicate the calibration status.

Sensor Alignment
Once you are confident with the calibration of the Faro Arm the next step is to align the Laser scanner with the Faro Arm. This is done within the Model Maker software package.

FE geometric face-fitting to bone data
Once your data has been acquired you need to align it with your mesh using CMGUI. This process is described in example 21d. Once this is done you need to convert your frontend files to backend files. A script is given in the Notes section below that can be used to convert 8-noded volume elements from the .exelem to .ipelem format (trilinearex2ipelem.pl). Once your files are in this backend format you are ready to use cm to fit the initial linear mesh to the data using tri-cubic Hermite basis functions. The command file and necessary files for this process are given in the example below.

Notes
The following perl scripts are used in various aspects of the above process:

#!/usr/local/perl5.6/bin/perl

# ........................................................................

# Name: asci2ex.pl

# Description: Reads the asci file format exported by the faro-mounted laser
# scanner (Model Maker) which contains coordinates of data points acquired by the
# scanner and converts it to an .exdata file.

# Updates:

# Usage: asci2ex.pl input_file[.asc]

# ........................................................................


# get input arguments
$input_argument = $ARGV[0];
$input_file = $input_argument . ".asc";
$output_file = $input_argument . ".exdata";
$node_value = "1";

# Open .asc file.
open (input,$input_file) or die "Error opening ";

# Open new exdata output file.
open (output, ">$output_file") or die "Error opening ";

# Print .exdata header to output

print output " Group name: laser_scanned_data\n";
print output " #Fields=1\n";
print output " 1) coordinates, coordinate, rectangular cartesian, #Components=3\n";
print output "  x.  Value index=1, #Derivatives=0, #Versions=1\n";
print output "  y.  Value index=2, #Derivatives=0, #Versions=1\n";
print output "  z.  Value index=3, #Derivatives=0, #Versions=1\n";

# Read Input and assign first value on each line as x coordinate,
# second as y coordinate and third as z coordinate. Print these and Node value.

  while () {
      chomp;
      ($x_coord, $y_coord, $z_coord) = split;
      print output " Node: $node_value\n";
      print output " $x_coord\n";
      print output " $y_coord\n";
      print output " $z_coord\n";
      $node_value++; 
    }

# Close the files.

close input;
close output;

#! /bin/perl
#
# ........................................................................

# Name: trilinearex2ipelem.pl

# Description: Reads a frontend .exelem file containing 8-noded elements and converts that file to the backend. Very simple e.g. does not handle versions.

# Updates:

# Usage: trilinearex2ipelem.pl exelemfile ipelemfile

# ........................................................................

use strict;

if (scalar @ARGV < 2)
  {
	 print "Usage $0 exelemfile ipelemfile\n";
	 die;
  }

my $filenameIn = $ARGV[0];
my $filenameOut = $ARGV[1];
my $data;
my $number_location;
my $element_number;
my $number_of_nodes;
my $number_of_elements;
my $collapsed_basis;
my $one;
my $two;
my $three;
my $four;
my $five;
my $six;
my $seven;
my $eight;

open (INPUT_FILE,"<$filenameIn") || (die "Could not open file $filenameIn"); 
open (OUTPUT_FILE, ">$filenameOut") || (die "Could not open output file $filenameOut");
$number_of_elements=0;

print OUTPUT_FILE <<END_HEADER;
 CMISS Version 1.21 ipelem File Version 2
 Heading:

END_HEADER
print OUTPUT_FILE " The number of elements is [1]: ";
$number_location = tell OUTPUT_FILE;
print OUTPUT_FILE "XXXXX\n\n";

while (defined($data = <INPUT_FILE>))
  {
	 if ($data =~ m/#Nodes=\s*(\d+)/)
		{
		  $number_of_nodes = $1;
		}
	 if (($data =~ m/Element:\s*([0-9]+).*/) && ($1 > 0))
		{
		  #Only want top level elements.

		  $element_number = $1;
		  #print "found element $1\n";
		  if ($number_of_nodes == 8)
			 {
				$number_of_elements++;
				print OUTPUT_FILE " Element number [1]:      $element_number\n";
				
				#Read until we get the nodes
				while ((defined ($data = <INPUT_FILE>)) && !($data =~ m/Nodes:/))
				  {
				  }
				#Read the node numbers
				if ((defined ($data = <INPUT_FILE>)) && 
				  ($data =~ m/\s*([0-9]*)\s*([0-9]*)\s*([0-9]*)\s*([0-9]*)\s*([0-9]*)\s*([0-9]*)\s*([0-9]*)\s*([0-9]*)/))
				  {
					 $collapsed_basis = 0;
					 $one = $1;
					 $two = $2;
					 $three = $3;
					 $four = $4;
					 $five = $5;
					 $six = $6;
					 $seven = $7;
					 $eight = $8;
						  print OUTPUT_FILE <<END_ELEMENT;
 The number of geometric Xj-coordinates is [3]: 3
 The basis function type for geometric variable 1 is [1]:  1
 The basis function type for geometric variable 2 is [1]:  1
 The basis function type for geometric variable 3 is [1]:  1
 Enter the 8 global numbers for basis 1:  $one   $two   $three   $four   $five   $six   $seven   $eight

END_ELEMENT
              }
				else
				  {
					 die "Node numbers not found for element $element_number";
				  }
			 }
		  else
			 {
				die "Fatal error: $number_of_nodes nodes for the top level element is not supported.\n";
			 }
		}
  }
close INPUT_FILE;

seek (OUTPUT_FILE, $number_location, 0);
printf (OUTPUT_FILE "%5d", $number_of_elements);
close OUTPUT_FILE;

Results can be view with view.com.

The comfile run by this example is as follows:

$MULTIPLEFITS=0

FEM def para;r;sternum;example
FEM def coor;r;sternum;example
FEM def base;r;sternum;example
#
# Reads in an initial mesh and calculates linear approximations
#  for the derivatives
#
FEM def node;r;sternum;example
FEM def elem;r;sternum;example
FEM up node deriv 1 linear
FEM up node deriv 2 linear
FEM up node deriv 3 linear
FEM update scale_factor normalise
#
# The initial mesh
#
FEM export node;sternum_init as sternum_init offset 1000
FEM export elem;sternum_init as sternum_init offset_elem 1000
#
# Calculate the intial xi projections
#
FEM def data;r;sternum;example
FEM def xi;c closest_face search 3 external
FEM li data err 
FEM export data;sternum_data as sternum_data error
#
# Setup the fitting problem
#
FEM up field from geometry
FEM def fit;r;sternum;example geometry

Additional testing commands:

#
# First fit
#
FEM fit
FEM update node fit
FEM update scale_factor normalise 
FEM def xi;c closest_face search 3 external
FEM li data err 
FEM export data;sternum_data_1 as sternum_data_1 error
FEM export node;fitted_sternum_1 as fitted_sternum_1 offset 2000
FEM export elem;fitted_sternum_1 as fitted_sternum_1 offset_elem 2000

if ($MULTIPLEFITS==1)
{

#
# Second fit
#
FEM fit
FEM update node fit
FEM update scale_factor normalise
FEM def xi;c closest_face external old
FEM li data err 
FEM export data;sternum_data_2 as sternum_data_2 error
FEM export node;fitted_sternum_2 as fitted_sternum_2 offset 3000
FEM export elem;fitted_sternum_2 as fitted_sternum_2 offset_elem 3000
#
# Third fit
#
FEM fit
FEM update node fit
FEM update scale_factor normalise
FEM def xi;c closest_face external old
FEM li data err 
FEM export data;sternum_data_3 as sternum_data_3 error
FEM export node;fitted_sternum_3 as fitted_sternum_3 offset 4000
FEM export elem;fitted_sternum_3 as fitted_sternum_3 offset_elem 4000
#
# Fourth fit
#
FEM fit
FEM update node fit
FEM update scale_factor normalise
FEM def xi;c closest_face external old
FEM li data err 
FEM export data;sternum_data_4 as sternum_data_4 error
FEM export node;fitted_sternum_4 as fitted_sternum_4 offset 5000
FEM export elem;fitted_sternum_4 as fitted_sternum_4 offset_elem 5000
#
# Fifth fit
#
FEM fit
FEM update node fit
FEM update scale_factor normalise
FEM def xi;c closest_face external old
FEM li data err 
FEM export data;sternum_data_5 as sternum_data_5 error
FEM export node;fitted_sternum_5 as fitted_sternum_5 offset 6000
FEM export elem;fitted_sternum_5 as fitted_sternum_5 offset_elem 6000
#
# Sixth fit
#
FEM fit
FEM update node fit
FEM update scale_factor normalise
FEM def xi;c closest_face external old
FEM li data err 
FEM export data;sternum_data_6 as sternum_data_6 error
FEM export node;fitted_sternum_6 as fitted_sternum_6 offset 7000
FEM export elem;fitted_sternum_6 as fitted_sternum_6 offset_elem 7000

FEM def no;w;fitted_sternum_6
FEM def elem;w;fitted_sternum_6

quit

}

Files used by this example are:

Name             Modified     Size

example_21i.com  03-May-2005  803 
sternum.ipbase   29-Aug-2002  3.3k
sternum.ipcoor   29-Aug-2002  570 
sternum.ipdata   29-Aug-2002  141k
sternum.ipelem   29-Aug-2002   13k
sternum.ipfit    13-Apr-2007  1.8k
sternum.ipnode   29-Aug-2002   93k
sternum.ippara   12-Nov-2002  5.9k
test_output.com  03-May-2005  2.0k

Download the entire example:

Name                      Modified     Size

examples_2_21_21i.tar.gz  27-Aug-2010  575k

Testing status by version:

	Status	Tested	Real time (s)
hpc_cm64_irix	Success	Thu Apr 1 13:32:37 2004	328
hpc_cm_irix	Success	Fri Jul 29 01:38:59 2005	180
hpc_cmo64_irix	Success	Thu Apr 1 13:27:09 2004	76
hpc_cmo_irix	Success	Sun Jul 31 01:39:46 2005	42
rs6000-aix
cm	Success	Wed Mar 4 01:15:34 2009	12
cm-debug	Success	Mon Mar 2 01:32:10 2009	123
cm64	Success	Wed Mar 4 01:15:36 2009	12
cm64-debug	Success	Tue Mar 3 01:38:41 2009	138
x86_64-linux
cm	Success	Sun Mar 6 00:01:16 2016	4
cm-debug	Success	Sat Mar 5 00:04:29 2016	19

Testing status by file:

cmiss_test.log

Success	hpc_cm64_irix:	cmiss_test.log.retain.
Success	hpc_cm_irix:	cmiss_test.log.retain.
Success	hpc_cmo64_irix:	cmiss_test.log.retain.
Success	hpc_cmo_irix:	cmiss_test.log.retain.
rs6000-aix
Success	cm:	cmiss_test.log.retain.
Success	cm-debug:	cmiss_test.log.retain.
Success	cm64:	cmiss_test.log.retain.
Success	cm64-debug:	cmiss_test.log.retain.
x86_64-linux
Success	cm:	cmiss_test.log.retain.
Success	cm-debug:	cmiss_test.log.retain.

cmiss_test.out

Success	hpc_cm64_irix:	ndiff test: no significant differences with generic answer.
Success	hpc_cm_irix:	ndiff test: no significant differences with generic answer.
Success	hpc_cmo64_irix:	ndiff test: no significant differences with generic answer.
Success	hpc_cmo_irix:	ndiff test: no significant differences with generic answer.
rs6000-aix
Success	cm:	ndiff test: no significant differences with generic answer.
Success	cm-debug:	ndiff test: no significant differences with generic answer.
Success	cm64:	ndiff test: no significant differences with generic answer.
Success	cm64-debug:	ndiff test: no significant differences with generic answer.
x86_64-linux
Success	cm:	ndiff test: no significant differences with generic answer.
Success	cm-debug:	ndiff test: no significant differences with generic answer.

Graphical output from this problem is given here.

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

Input last modified: Thu Aug 26 13:17:08 2010

CMISS Help / Examples / 2 / 21 / 21i