Example 21d: FE geometric fitting to laser scanned bone data

This example covers:

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

Created by: Steven Thrupp 6 Feb 2001

Modified: Steven Thrupp 22 Aug 2002

Creating an initial mesh: using the Faro Arm

This process is quite straightforward. First you physically draw a mesh on your object and then rigidly mount the bone and obtain the intersecting points on the mesh as nodes using the Faro Arm. Then connect the nodes using CMGUI and this results in a linear initial mesh. Currently the Faro Arm is connected solely to ESU30 and so users will need to obtain a login on this computer in order to use the Faro Arm.

The Faro Arm is straightforward to use. In the directory from which you start CMGUI you will need a file called faro.init. The content of this file is given below.
FaroInitialisationFile
serial_port 2
baud_rate 9600
_D1
        
This file enables the Faro Arm to be initialized when CMGUI is started so that you can obtain points directly via a CMGUI session. You need to ensure that only one CMGUI session at a time is using the Faro Arm.

Once the Faro Arm is initialised when you start CMGUI open up the 3d digitizer window using either the CMGUI drop-down menu or using GFX create 3d_digitizer.

The default settings in the 3d digitizer window are to obtain points as position and as node.

Create a node group (
gfx cre ngroup femur) and then switch on the input device (click over `Input Device' in the Data Grabber window). You will now be able to obtain node points to add to the group you have selected in `Group'.

To obtain points as data (i.e. .exdata), create a data group (
gfx cre dgroup femur) and instead of `Data as Node', select `Data as Point'.

For some applications the user may also wish to place the Faro Arm on Stream Mode. Do this by selecting the `Stream' box at the bottom of the Data Grabber window and hit `record'. (Note that a point seems to be automatically acquired at the position the Faro Arm is placed when the record button is depressed. This is a minor bug and your .exnode or .exdata file should be edited to remove this point).


Data Grabber/3d_digitizer Window in CMGUI

Once you have undertaken the above steps you are ready to acquire data. Note that any of the four buttons on the Faro Arm may be depressed to obtain a point. Also note that the green diode signifies the Arm is switched on. The red diodes are error messages and indicate that one or more of the swivel-points on the Arm is fully extended. When this occurs the Arm is unable to obtain points and the Arm will need to be rearranged so that no swivel-point is fully extended.

Calibrating the Faro Arm

Periodically the Faro Arm needs to be calibrated and also requires calibration when one of the heads is changed. In order to calibrate the Faro Arm the user needs to become user Faro on ESU30 and use SoftWindows ('swin'). Once SoftWindows is started open up the Faro Arm program on the desktop. Go to the settings menu and from there you can either choose the Point Probe and or the Ball Probe (by choosing Ball Probe you will need to specify the diameter of the probe in the diameter field). In the same box you will see a calibrate sub-box. Choose 1"Ball and follow the directions to click 27 points on the one inch ball physically located at the base of the Faro Arm. Finally, another setting which you may need to check from time to time is the Baud Rate which should be set to 9600. This is found also in the Settings Menu under `Communications'.

Obtaining Data for Fitting: Using the Hand Held Laser Scanner

Once you have created your initial mesh using nodes obtained using the Faro Arm and used CMGUI to create your linear elements you now require data to fit your mesh to. For the musculoskeletal modelling project this was done using the 3Space® FastSCAN™ Hand Held Laser Scanner. However at the time of revision of this page the HLS has been superseded by the Faro-Mounted ModelMaker4 Scanner. Use of this later scanner will be covered elsewhere (and will eventually be linked to from this example).

The HLS is held by Andrew Pullan he should be consulted over its use. It requires an NT operating system.

First you need to connect the scanner. The scanner box is connected to the computer's parallel port. The HLS wand, receiver and FastSCAN magnetic tracking device should be connected to the scanner box prior to the box being switched on (otherwise there seems to be interference which results in poor quality data acquisition). When the hardware is connected the scanning software can be opened from the computer's main menu (look for Polhemus FastSCAN in the menu).

Before starting the scanner you need to mount the object. The HLS uses a magnetic tracking device hence any metal object within a metre of the device will distort the magnetic field. All metal objects therefore need to be kept outside of this distance. Use the wand to obtain points. Set the transmitter as reference in the `scanner' menu. Try to minimise the number of sweeps you need to make. It is important to mount the object such that each face of it is easily accessible by the two cameras of the wand. Once you have obtained a surface you are happy with you should export the surface as an `.obj' file. Select export in the `file' menu and export withthe settings `tracker reference', `cloud of points' and format as unix text. You are able to set export parameters. Generally it is wise to combine facets, set the resolution you wish, and limit objects since this last function will remove extraneous points that are invariably incidentally picked up when you are scanning.

(A brief note on resolution. Optimally the export surface ought to coincide with the raw data. However often this results in an unreasonably large exported file (containing, for example over 100,000, individual points). I tended to use a profile smoothing set to low (and tended to use an export resolution of 1 mm and filter radius of 2.5 mm. See the documentation on the scanner for information on the resolution options for the scanner (including user examples).

Once you Transfer your exported .objfile to your SGI, use the executable obj2ex fileName.obj fileName.exdata. A copy of obj2ex can be found in /usr/people/thrupp/.bin/. You will probably need to manually edit the resultant .exdata file to remove any obviously extraneous points.


FE geometric fitting to bone data

Once you are happy with your data it is necessary to transform it so that it aligns with your mesh. This can be done in CMGUI using the scene editor which contains a `transform' option for each group. When you have used this to move the object within your scene you need to use the command gfx edit graph groupName applywhen you are satisfied that your data and mesh are closely aligned and then write out the transformed data.

This example has covered the process from creating an initial mesh by drawing the mesh on an object, then obtaining those mesh points using the Faro Arm, followed by using CMGUI to create and initial linear mesh by connecting the points, followed by obtaining and converting data using the HLS. Then CMGUI is again used to align the data and mesh. Once all of this is done, you need to convert your frontend files to backend files and use
cm to fit the initial mesh to the data using cubic Hermite basis functions. The command file and necessary files for this process are given in the example below.

Notes

Initial Mesh Fitted Mesh

Results can be view with view.com.

This example borrows from Leo Cheng's example on finite element geometric fitting to epicardial data. (Although his example uses simplex elements rather than this example which uses collapsed node elements with versions at the apex and base.)


The comfile run by this example is as follows:

$MULTIPLEFITS=0

{
  eval { require Time::HiRes };
  if( $@ )
    {
      print "warning: $@";
    }
  else
    {
      import Time::HiRes qw(time);
    }
}

sub cpu_time()
  {
    # Add user/system parent/child cpu_times.
    my $result = 0;
    map { $result += $_ } times;
    $result;
  }

sub print_timing($$$)
  {
    my ($message,$wall,$cpu) = @_;
    printf "%s: %.3fs wall, %.3fs cpu.\n", $message, $wall, $cpu;
  }

sub print_memory_use()
  {
    # Different operating systems have different data and ps uses different
    # arguments so try a series of commands in an attempt to collect the
    # information.
    system qw(ps -o vsz -p), $$;
    system qw(ps -l -p), $$;
    system qw(ps v), $$;
  }

FEM def para;r;femur;example
FEM def coor;r;femur;example
FEM def base;r;femur;example

#
# Reads in an initial mesh and calculates linear approximations
#  for the derivatives
#
FEM def node;r;femur;example
FEM def elem;r;femur;example
FEM up node deriv 1 linear node 2..401
FEM up node deriv 2 linear node 2..401
FEM update scale_factor normalise

#
# The initial mesh
#
FEM export node;femur_init as femur_init offset 10000
FEM export elem;femur_init as femur_init offset_elem 10000


#
# Calculate the intial xi projections
#
FEM def data;r;femur;example
FEM def xi;c close
#FEM def xi;w;femur 
FEM li data err 
FEM export data;femur_data as femur_data error


#
# Setup the fitting problem
#
FEM up field from geometry
FEM def fit;r;femur;example geometry

($wall_setup,$cpu_setup) = (time, cpu_time);
print_timing 'setup', $wall_setup - $^T, $cpu_setup;

#
# First fit
#
FEM fit
FEM update node fit
FEM update scale_factor normalise 
FEM def xi;c close old
FEM li data err 
FEM export data;femur_data_1 as femur_data_1 error
FEM export node;fitted_femur_1 as fitted_femur_1 offset 20000
FEM export elem;fitted_femur_1 as fitted_femur_1 offset_elem 20000

if( $TESTING ) {
  fem list data;fitted1 errors;
  fem list node;fitted1
}

if ($MULTIPLEFITS==1)
{

#
# Second fit
#
FEM def fit;r;femur_2;example geometry


FEM fit
FEM update node fit
FEM update scale_factor normalise
FEM def xi;c close old
FEM li data err 
FEM export data;femur_data_2 as femur_data_2 error
FEM export node;fitted_femur_2 as fitted_femur_2 offset 30000
FEM export elem;fitted_femur_2 as fitted_femur_2 offset_elem 30000

#
# Third fit
#
FEM fit
FEM update node fit
FEM update scale_factor normalise
FEM def xi;c close old
FEM li data err 
FEM export data;femur_data_3 as femur_data_3 error
FEM export node;fitted_femur_3 as fitted_femur_3 offset 40000
FEM export elem;fitted_femur_3 as fitted_femur_3 offset_elem 40000

#
# Fourth fit
#
FEM fit
FEM update node fit
FEM update scale_factor normalise
FEM def xi;c close old
FEM li data err 
FEM export data;femur_data_4 as femur_data_4 error
FEM export node;fitted_femur_4 as fitted_femur_4 offset 50000
FEM export elem;fitted_femur_4 as fitted_femur_4 offset_elem 50000
#
# Fifth fit
#
FEM fit
FEM update node fit
FEM update scale_factor normalise
FEM def xi;c close old
FEM li data err 
FEM export data;femur_data_5 as femur_data_5 error
FEM export node;fitted_femur_5 as fitted_femur_5 offset 60000
FEM export elem;fitted_femur_5 as fitted_femur_5 offset_elem 60000
#
# Sixth fit
#
FEM fit
FEM update node fit
FEM update scale_factor normalise
FEM def xi;c close old
FEM li data err 
FEM export data;femur_data_6 as femur_data_6 error
FEM export node;fitted_femur_6 as fitted_femur_6 offset 70000
FEM export elem;fitted_femur_6 as fitted_femur_6 offset_elem 70000

FEM def no;w;fitted_femur_6
FEM def elem;w;fitted_femur_6

}

($wall_solve,$cpu_solve) = (time, cpu_time);
print_timing 'solve', $wall_solve - $wall_setup, $cpu_solve - $cpu_setup;
print_memory_use;

Files used by this example are:

Name             Modified     Size

example_21d.com 29-Sep-2005 3.6k femur.ipbase 23-Aug-2002 1.5k femur.ipcoor 23-Aug-2002 688 femur.ipdata 23-Aug-2002 206k femur.ipelem 23-Aug-2002 147k femur.ipfit 13-Apr-2007 3.8k femur.ipnode 23-Aug-2002 393k femur.ippara 28-Sep-2005 5.9k femur_2.ipfit 13-Apr-2007 15k

Download the entire example:

Name                      Modified     Size

examples_2_21_21d.tar.gz 27-Aug-2010 815k

Testing status by version:

StatusTestedReal time (s)
hpc_cm64_irixSuccessThu Apr 1 14:19:18 2004358
hpc_cm_irixSuccessFri Jul 29 01:54:39 2005309
hpc_cmo64_irixSuccessThu Apr 1 14:16:51 2004201
hpc_cmo_irixSuccessSun Jul 31 01:48:17 2005166
rs6000-aix
cmSuccessWed Mar 4 01:19:50 200921
cm-debugSuccessMon Mar 2 01:40:41 2009176
cm64SuccessWed Mar 4 01:20:06 200923
cm64-debugSuccessTue Mar 3 01:46:57 2009192
x86_64-linux
cmSuccessSun Mar 6 00:01:22 20166
cm-debugSuccessSat Mar 5 00:05:16 201627

Testing status by file:


Graphical output from this problem is given here.


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

Input last modified: Thu Aug 26 13:09:40 2010


CMISS Help / Examples / 2 / 21 / 21d