Femur cohort studies

The first step is to load an individual's femur surface and also the canonical surface onto which we wish to map the cortical data (File->Load surfaces, select the two .ply|.wrl files by clicking on the first and ctrl-clicking on the second, or just select the individual's .ply|.wrl file and let wxRegSurf look for a canonical .ply|.wrl in the same directory). It is important that the canonical surface's filename includes the string canonical: this tells wxRegSurf that this is the "low resolution" (LR, red) surface that needs registering to the "high resolution" (HR, green) surface. It is also important that the canonical surface's filename includes the string femur, since there are some minor differences in the way wxRegSurf's algorithms are tuned for different bone types. Although wxRegSurf offers a number of different options for the registration transformations, we shall describe just one combination that we have found to work well with proximal femurs. We start with a nine degree-of-freedom aspect ratio registration (Registration->Aspect), which allows a rigid body transformation and independent scaling in the x, y and z directions. wxRegSurf calculates a rough, initial alignment by analysing the surfaces' shapes and topologies. If this fails, directly manipulate the LR surface (Mouse->LR surface, left/middle mouse button to rotate, right mouse button to translate, scrollwheel to scale, single click to reflect) until the surfaces are roughly aligned.
By default, wxRegSurf matches each vertex on the LR surface with the closest vertex on the HR surface, as long as the distance between the vertices is less than the Search range. To see these matches, click Show matches. wxRegSurf's registration algorithms aim to minimize the summed lengths of these matches. In this example, we see that there is a lot of work to be done. The individual's femur is a lot larger than the canonical femur and also different in shape. As well as being larger, the individual scan also extends much further down the femoral shaft. This is important: the canonical femur stops at the lesser trochanter and the individual scan must extend at least this far if we are to map its cortical data across.
Next, we set the Iterations slider to something suitably high and press Start. At each iteration, wxRegSurf finds the transformation that minimizes the summed lengths of the matches and then recomputes the matches. After a short while this process converges to a stable alignment and we can either interrupt it (Stop) or wait until the iteration count is reached.
The nine degree-of-freedom aspect ratio transformation has brought the surfaces closer together, but we need to apply some nonrigid adjustments to finish the job. We could perform the nonrigid stage from this starting point (Registration->Aspect+FFD or Registration->Aspect+LAD), but then the individual femur's shape would be encoded in both the initial aspect ratio transformation and the subsequent nonrigid transformation. Often we are interested in the individual's shape as well as the cortical parameters, and for this reason we suggest selecting Registration->Similarity and then pressing the Aspect->Similarity button. This constructs a seven degree-of-freedom similarity transformation from the current aspect ratio transformation, preserving the z scaling but adjusting the x and y scales to match. Compared with a direct similarity registration, the Aspect->Similarity route tends to do a better job of aligning the lesser trochanters.
We are now in a position to apply a nonrigid transformation which will, incidentally, fully encode the individual's shape (but not size, since this depends also on the similarity transformation). There are three options: Similarity+FFD, Similarity+LAD and Similarity+SSM, where FFD indicates a B-spline free form deformation, LAD indicates a locally affine deformation and SSM indicates a statistical shape model. We did all our early work with FFDs before switching to LADs, which we found to perform better at the shaft rim where there is relatively little information to guide the alignment.

The SSM option is available only when a valid statistical shape model has been loaded or calculated. We provide a suitable .ssm file containing a model trained on registrations of many femurs (the modes can be visualized using the controls in the Utilities tab). SSM registration is faster but less accurate than direct LAD registration.

To perform the nonrigid registration, we set the Iterations slider to around 20 and select Registration->Similarity+LAD (preferable) or Registration->Similarity+SSM (if in a rush), before pressing Start. When complete, File->Save registration saves the alignment in .reg and .breg files.

An unsuccessful registration may require some manual intervention. The most common failure mode involves misalignment of the lesser trochanters, with the canonical surface riding too high up the individual. If this happens, press the Reset button to reinitialize the LAD parameters, select Mouse->LR surface to directly manipulate the LR surface, use the scrollwheel to enlarge the LR surface until the lesser trochanters are approximately level, then run 20 or so Similarity+LAD iterations.

The next stage of the process is to load the cortical data for the HR surface (File->Load surface data, choose between the sw thickness file, sw density file and surface data file filters, then select the appropriate .bin file). wxRegSurf will quietly load any corresponding error file it finds in the same directory: these contain precision estimates at each vertex and are used to weight any data smoothing we might perform. The filenames for the cortical data must terminate in either _thickness.bin (Stradview thickness), _density.bin (Stradview density) or _dat.bin (general surface data).

To see the data we have loaded, we select Surface shading->Surface data. In this example, we see a cortical thickness estimate obtained by deconvolution in Stradview v4.6. The white areas are where the deconvolution was not successful.

wxRegSurf maps each LR data value with its matched HR value (Show matches). This sort of nearest neighbour mapping might at first seem naive, but it is in fact perfectly adequate if used carefully in conjunction with suitable smoothing. Data smoothing is, in any case, a prerequisite for cohort analysis techniques like statistical parametric mapping (SPM). Every time we press the Smooth button, wxRegSurf replaces each data value with a weighted average of immediately neighbouring values, with the weights set according to the precision estimates found in the error file. If there is no error file, uniform weighting is used.

How much smoothing is necessary? For SPM, this depends on the resolution of the CT data and the expected effect size. Suppose we have determined that we need a 10mm full-width-half-maximum (FWHM) kernel. We need to calculate how many smoothing cycles are required to achieve this. The number will depend on the average edge length in the mesh. wxRegSurf helps by displaying an estimate of the current FWHM after each smoothing cycle. The image on the left shows the result of three smoothing cycles.

The next step is to map the smoothed HR data across to the LR surface. This is achieved by pressing the Map button. To visualize the mapped data clearly, without the obscuring clutter of the HR data, we uncheck the Show HR checkbox.

Finally, we perform four further smoothing cycles on the LR mesh, achieving (in this instance) approximately 10mm FWHM overall. Note that wxRegSurf smoothes any data that is currently displayed, so if we only wish to smooth the LR data we must ensure that the Show HR checkbox is unchecked. The smoothed and mapped data can be saved by selecting File->Save LR data. The saved files provide the input for subsequent SPM analysis.

Splitting the smoothing between the HR and LR meshes is something we do routinely. Sufficient smoothing should be performed on the HR mesh to (a) remove all the invalid (white) data patches and (b) ensure that the subsequent nearest neighbour mapping stage adequately captures all the HR data in the vicinity of each LR vertex (wxRegSurf will warn if there is a danger of undersampling). Performing the rest of the smoothing on the LR mesh removes any noise caused by spatial imprecision in the nearest neighbour mapping (the matches are not necessarily orthogonal to either surface) and minimizes any potential bias caused by different average edge lengths in the individual meshes.

We might wish to repeat this process for tens or hundreds of bones. wxRegSurf provides a mechanism for sequential processing using batch (.run) files. An example batch file is shown below. We could cut-and-paste these lines into a file called, for example, MyBatchFile.run. A batch file can be run from the user interface (File->Load and run batch file) or by supplying its name as a command line argument when wxRegSurf is invoked.

/path/to/surface/and/data/files/including/trailing/slash/

# LR surface              HR surface      LR data  HR data                   reflect  smth LR  smth HR  range  fussy    n iter  reg type  save reg  save LR data  save map  review       rims         save HR  match
#                                                                                                              -1=auto          -=as->sm            -1 = no                 0 = no       0 = no       surface  colours
#                                                                                                                                                   0 = yes                 1-6 = angle  1 = yes
#                                                                                                                                                   >0 = smth               - = LR only  2 = to rims

# aspect ratio registration
canonical_femur_5580.wrl  individual.wrl  n/a      n/a                       1        0        0        999    0        200     7         1         -1            0         0            1            0        0

# similarity+LAD registration
canonical_femur_5580.wrl  individual.wrl  n/a      n/a                       1        0        0        999    0        20      -3        1         -1            0         0            1            0        0

# review
canonical_femur_5580.wrl  individual.wrl  n/a      n/a                       1        0        0        999    0        200     3         0         -1            0         1            1            0        0

# map
canonical_femur_5580.wrl  individual.wrl  n/a      individual_thickness.bin  1        0        3        999    0        0       3         0         4             0         0            1            0        0

A batch file may include any number of blank or comment lines. The first non-comment line contains the full path to the working directory, where wxRegSurf will look for .ply|.wrl and .bin files. Each subsequent line contains a set of instructions for wxRegSurf to load files, process them and then save the results. Draft batch files can be produced using wxRegSurf's automatic batch file writer, then edited manually as required.

In the example above, the first line loads two surface files (but no thickness or density files) and registers them with an aspect ratio transformation (reg type is 0 for rigid, 1 for similarity, 2 for similarity+FFD and so on, mirroring the options in the Registration choice control). Since the canonical surface is a right femur and the individual scan is of a left femur, we need to reflect the canonical surface before starting: hence the reflect flag being set to 1. The search range is effectively unconstrained (999). 200 registration iterations are performed before the results are saved in a .reg file.

The second line performs 20 iterations of a similarity+LAD registration, after first mapping the aspect ratio transformation to the similarity transformation (this is indicated by the minus sign in front of the 3 in the reg type field). The results are saved in .reg and .breg files. We could fill the batch file with hundreds of similar lines, registering different bones to the canonical surface, and leave wxRegSurf to run unattended.

Following registration, we must check that all bones were aligned satisfactorily. This is what the third line does in the example above. It performs 200 iterations in "review" mode, whereby the surfaces are loaded and spun around for the user to inspect. Lower numbers of iterations reduce the inspection time, higher numbers increase it. Alternatively, we could set the iteration count to zero, in which case wxRegSurf instead saves a snapshot of the aligned bones (.png) for offline inspection. Changing the review flag from 1 to some other number saves a different view (1=front, 2=left, 3=back, 4=right, 5=top, 6=bottom). Any misregistered surfaces must be corrected manually, with the updated results saved to the .reg and .breg files.

The final line of the above example loads the HR cortical data, performs three smoothing cycles on the HR data, applies the similarity+LAD transformation, matches the LR and HR vertices over an effectively unconstrained search range, maps the smoothed HR data across to the LR surface, performs a further four smoothing cycles on the LR surface and then saves the smoothed and mapped data to a .bin file.