Vertebra cohort studies

The first step is to load an individual's vertebra 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 vertebra, 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 L1 vertebrae. We start with a seven degree-of-freedom similarity transformation (Registration->Similarity). wxRegSurf calculates a rough, initial registration by aligning the centres of gravity and assuming that both vertebrae are oriented identically in the CT coordinate system. 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) 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.
wxRegSurf also offers the option of stipulating fussy matches. This involves looking beyond the closest HR vertex and favouring instead HR vertices with similar surface normals and curvatures. In the example on the left, we have selected 15% fussy matches using the Fussy % slider. This forces fussy matches for the least well aligned 15% of the LR vertices. Show matches displays the fussy matches in magenta and the standard matches in white.

By definition, fussy matches extend over a longer range than regular matches. They can therefore introduce substantial errors into the alignment if the matches are wrong. Nevertheless, some degree of fussy matching is necessary when registering vertebrae, since otherwise the narrow processes tend to be incorrectly matched and aligned.

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 seven degree-of-freedom similarity transformation has brought the surfaces closer together, but we need to apply some nonrigid adjustments to finish the job. Fussy matching is especially important at this stage. Looking closely at the matches, it is apparent that the tips of the LR transverse processes are correctly matched to the tips of the HR transverse processes by fussy matches (magenta). Standard matching would pair the LR tips with the closest vertices on the HR processes, which are definitely not the HR tips. There would then be nothing to "pull" the LR processes outwards.

There are three options for the nonrigid registration: 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. wxRegSurf's low resolution FFDs do not allow sufficient local control for vertebra alignment, so we prefer LADs.

The SSM option is available only when a valid .ssm file is found in the same directory as the canonical .ply|.wrl file. It fits the first n modes of a statistical shape model trained on similarity+LAD registrations of many vertebrae (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, check the Auto sequence checkbox and select Registration->Similarity+LAD (preferable) or Registration->Similarity+SSM (if in a rush), before pressing Start. The automatic sequencing gradually reduces the search range and the proportion of fussy matches as the algorithm converges, thus minimizing the impact of incorrect fussy matches that span large distances. 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 the transverse processes. The initial search range in the automatic sequence is conservative, since incorrect fussy matches over long distances are particularly problematic at the nonrigid stage, when they can result in surface folding. But sometimes this restricted search range is insufficient to match up the tips of the transverse processes. If this happens, press the Reset button to reinitialize the LAD parameters, reset the Fussy % slider to 15, increase the search range until the tips are correctly matched, then run one or two Similarity+LAD iterations to bring the tips closer together. Next, check the Auto sequence checkbox and run a further 20 or so 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 (Stradwin thickness), _density.bin (Stradwin 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 Stradwin v4.6. The white areas are where the deconvolution was not successful.

wxRegSurf maps each LR data value with its matched HR value. At this stage, we set Fussy % to zero, ensuring that each LR vertex is matched to the nearest HR vertex (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 sixteen 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, 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.


# 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
#                                                                                                                     -1=auto          -=as->sm            -1 = no                 0 = no       0 = no
#                                                                                                                                                          0 = yes                 1-6 = angle  1 = yes
#                                                                                                                                                          >0 = smth               - = LR only  2 = to rims

# similarity registration
canonical_L1_vertebra_8329.wrl  individual.wrl  n/a      n/a                       0        0         0        100    15       200     1         1         -1            0         0            1

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

# review
canonical_L1_vertebra_8329.wrl  individual.wrl  n/a      n/a                       0        0         0        100    0        200     3         0         -1            0         1            1

# map
canonical_L1_vertebra_8329.wrl  individual.wrl  n/a      individual_thickness.bin  0         0        3        100    0        0       3         0         16            0         0            1

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 a similarity 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). The search range is 10mm (Search range slider set to 100) and there are 15% fussy matches. 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 with automatic sequencing of the search range and the proportion of fussy matches (this is indicated by the -1 in the fussy 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). Adding a minus sign in front of the review flag saves a snapshot of the LR surface only, which is useful when checking for surface folding caused by incorrect fussy matches. 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 a 10mm search range, maps the smoothed HR data across to the LR surface, performs a further sixteen smoothing cycles on the LR surface and then saves the smoothed and mapped data to a .bin file.