Computational Complexity

In order to understand the computational difficulty of the problem it is helpful to consider the amount of calculations required. For instance, take a typical field of view (FOV) for the brain as 200mm cubed, and let each voxel be a cube of length d mm. Therefore, there are (200/d)³ voxels in the volume.

To interpolate this volume requires, generally, one evaluation of the interpolation function (acting on the floating image) for each voxel location in the reference image. For trilinear interpolation there are at least 8 floating point operations required for each interpolation evaluation. Therefore, the total cost of interpolation is at least $8 \times (200/d)^3$ operations.

In addition, to calculate the cost function typically requires at least one operation (to estimate variance) per voxel, making the number of operations required per voxel at least 9. Therefore, on a machine that could achieve 500 MFLOPS, the time required for one cost function evaluation, with voxel size of 1mm would be at least $9 \times 200^3 / (5 \times 10^8) = 0.144$ seconds. However, in practice there are some additional calculations required as well as memory access overheads, and so our implementation (on a Pentium3 500MHz) takes approximately 1 second for this evaluation.

Now consider implementing a search through the transformation space by constructing a grid with N different values for each parameter, so that there are a total of N^D elements (where D is the Degree Of Freedom -- 12 for affine). The total number of calculations required to evaluate the cost function at each point is $7.2 \times 10^7 (N^D / d^3)$ which would ideally take $0.144 \times (N^D / d^3)$ seconds at 500 MFLOPS. Table 1 shows the times required for a range of moderate parameters. From this table it can be seen that evaluations for 12 DOF using 1mm voxels are prohibitive, taking over 4000 years for a grid with N=10 and over 21 hours for a simple neighbourhood N=3. However, using 8mm cubed voxels improves the execution time although the biggest improvement comes from restricting the dimensionality of the transformation space (that is, the DOF).

 

Table 1: Some typical computational requirements for implementing a simple grid-based search strategy (see text for explanation of terms). The times are based on an idealised 500 MFLOP performance with no memory overhead. In practice our implementation was approximately 7 times slower than this.

N	D	d	No. of Operations	Time (sec)
10	12	1	$7.2 \times 10^{19}$ $\quad$	$1.44 \times 10^{11}$
10	12	2	$9.0 \times 10^{18}$ $\quad$	$1.8 \times 10^{10}$
10	12	4	$1.13 \times 10^{18}$ $\quad$	$2.3 \times 10^{9}$
10	12	8	$1.4 \times 10^{17}$ $\quad$	$2.8 \times 10^{8}$
10	6	1	$7.2 \times 10^{13}$ $\quad$	$1.44 \times 10^{5}$
10	6	2	$9.0 \times 10^{12}$ $\quad$	18 000
10	6	4	$1.13 \times 10^{12}$ $\quad$	2 250
10	6	8	$1.4 \times 10^{15}$ $\quad$	281
3	12	1	$3.8 \times 10^{13}$ $\quad$	76 528
3	12	2	$4.8 \times 10^{12}$ $\quad$	9 566
3	12	4	$6.0 \times 10^{11}$ $\quad$	1 196
3	12	8	$7.5 \times 10^{10}$ $\quad$	149

 

For a simple, but usable, search routine it would be necessary to search with a grid of at least N=100, in which case even for 6 DOF and 8mm voxels the time taken would exceed 7 hours. This would be unacceptable for most users, especially given that our implementation would be likely to take over 50 hours. However, in combination with a local optimisation method, a coarse search method using 8mm voxels can be used to create a useful optimisation method, which is pursued in section 4.3.