## Abstract

We present a novel, high-resolution magnetic resonance technique, fine structure analysis (FSA) for the quantification and analysis of amorphous and quasi-amorphous biological structures. The one-dimensional technique is introduced mathematically and then applied to one simulated phantom, two physical phantoms and a set of *ex vivo* biological samples, scanned with interpoint spacings of 0.0038–0.195 mm and cross-sectional sizes of 3 × 3 or 5 × 5 mm. The simulated phantom and one of the physical phantoms consists of randomly arranged beads of known size in two and three dimensions, respectively. The second physical phantom was constructed by etching lines on Perspex. The *ex vivo* samples are human bone specimens. We show that for all three phantoms, the FSA technique is able to elucidate the average spacing of the structures present within each sample using structural spectroscopy, the smallest of which was 180 µm in size. We further show that in samples of trabecular bone, FSA is able to produce comparable results to micro-computed tomography, the current gold standard for measuring bone microstructure, but without the need for ionizing radiation. Many biological structures are too small to be captured by conventional, clinically deployed medical imaging techniques. FSA has the potential for use in the analysis of pathologies where such small-scale repeating structures are disrupted or their size, and spacing is otherwise altered.

## 1. Introduction

Many disease states are characterized by the appearance of or disruption to fine structures. For example, the trabeculae in the bones of osteoporotic subjects tend to be less densely packed than those of normal subjects [1]. As a second example, early liver disease, is often manifest by the appearance of localized fibrotic structures [2]. As a third example, there is evidence that highly organized locally parallel neural structures become increasingly disorganized as degenerative neural disease progresses [3]. Each of these cases corresponds to a change in local spatial frequency: the changes are localized spatially and they also correspond to textural changes that can be characterized in the spatial frequency domain. In these and similar cases, the fine structures of interest span the range 0.05–0.2 mm.

There is a need for technology that can be deployed clinically to enable regions of interest with such local spatial frequencies to be detected *in vivo*. In many pathologies, such as the examples given above, therapies may achieve a better outcome if an early diagnosis is made, and physicians are often reluctant to subject a patient with borderline risk factors to an invasive and painful biopsy or other tests. Moreover, there is a need for repeated monitoring to evaluate treatment efficacy, not least, because there may be several potential therapies for a particular pathology, a selection of which needs to be evaluated. Repeated probing of the fine structures argues against technologies that use radiation such as computed tomography (CT), dual energy X-ray absorbtiometry (DEXA) and nuclear medicine; rather, it suggests the use of ultrasound or magnetic resonance imaging (MRI). In each of the three cases given above, ultrasound is challenged by access, signal-to-noise and spatial resolution. The obvious technology of choice is MRI, which meets all of the criteria except for spatial resolution. Several investigations have been performed to improve the resolution of MRI for various applications. For example, in the case of trabecular bone imaging, much smaller voxel sizes have been achieved, typically of the order of 100 µm^{3} [4,5]. Such techniques have generally been limited to peripheral sites because of signal-to-noise considerations; however, recent developments have shown progress in studying more central sites, for example [6]. In addition, magnetic susceptibility artefacts can cause significant errors when examining the features of images close to interfaces between tissues [7].

Other MRI methods to detect fine structures have also been reported. For example, diffusion tensor imaging (DTI) measures the diffusion of molecules *in vivo*, and because the diffusion processes often interact with tissue boundaries, structures of interest can be mapped. DTI is a true three-dimensional technique that can provide an image set of a whole brain *in vivo* with a scan of approximately 4 min on a 3 T machine [8]. Single shot DTI sequences are robust against motion effects, but can be susceptible to imaging artefacts [9,10].

This paper introduces the technical aspects of an MRI-based technique called fine structure analysis (FSA) [11] and applies it to a set of phantoms and *ex vivo* samples that are either (approximately) amorphous or varies only in one direction. Conventionally, MRI acquires a two-dimensional image by filling the whole of *k*-space (or a substantial fraction of it) then applies the Fourier transform to generate the image. Instead, in FSA, a single line of *k*-space is sampled finely, yielding a one-dimensional signal over the range of the finely sampled axis chosen by the clinician. The fine sampling rate depends upon the application and equipment available. The structural spectrum, which is a representation of the spatial frequencies present locally in the one-dimensional signal, is then computed. In essence, FSA trades the reduction from a two-dimensional image to a one-dimensional spatial spectrum at each selected anatomic location for higher-resolution sampling along the corresponding line of *k*-space, hence in the spatial frequency spectrum. Restricting attention to one line of *k*-space can potentially yield information at significantly higher resolution than is available in conventional imaging, and, in certain cases, this does not imply a loss of information. Furthermore, as resolution of MRI improves the resolution advantage of FSA relative to imaging will be maintained. In this paper, we present results that have interpoint spacings of 171 µm and 152.5 µm in a scan time of 3 and 1.5 min, respectively, on clinical 3 T MRI platforms. Furthermore, because the spatial spectrum is sensitive to changes in the pitch of repeating structures rather than changes in individual features of the image, the issue reported in the literature with magnetic susceptibility artefacts is mitigated. This is because the artefact added to a periodic signal will produce a spurious periodic signal at a frequency at least twice that of the fundamental signal that will typically be of sufficiently high frequency to be outside the range of interest of the clinician.

We note there are several approaches in the literature that have similarities to the method we present. For example, point resolved spectroscopy (PRESS) was developed for magnetic resonance spectroscopy which also consists of intersecting 90° and 180° RF pulses, with an additional orthogonal 180° RF pulse to completely limit the excited volume to a single voxel [12]. Although its implementation is similar to that of the pulse sequence we present, the PRESS sequence acquires data that are very different to the spatial information of interest in FSA. A number of other pulse sequences have been developed for single voxel spectroscopy; to the extent that these are similar in implementation to PRESS, they are also similar to the sequence we present [13]. Line scan imaging sequences form an image out of a set of one-dimensional excitations, so these sequences are most similar to our approach. Line scan approaches have significant advantages for robustness against motion artefacts [14], and have also been applied to diffusion sequences [15,16]. We are not interested in forming a two-dimensional image, however, and we are typically interested in fewer lines than would be required to form a line scan image. The lines will be treated separately and analysed with spatial frequency analysis rather than used to form an image.

There are two steps to the FSA process: (i) the generation of a set of spatial frequency signals and (ii) the (application-specific) analysis of those signals. These two steps are explained in detail in §§3 and 4, respectively. We stress that these two steps are largely independent, in the sense that the signal processing methods developed for the spatial frequency representation described in §4 can be applied to the signals acquired by any means of acquisition, including those outlined in §3. Section 5 provides descriptions of how the experiments were performed and presents the results from these experiments. We conclude with a brief discussion and outlook.

## 2. Amorphous structures and biological tissues

In materials science, diffraction patterns can be used to characterize the typical physical structure in a sample. Broadly, two classes of structures can be characterized in this way: crystalline and amorphous. Figure 1 illustrates an amorphous structure and its associated Fourier transform. The diffraction pattern gives an estimate of the spectral density function, which is related to the Fourier transform of the image of the structure in position space. The radius of the central ring is inversely proportional to the nearest-neighbour spacing of elements of the structure. In this way, the windowed Fourier transform gives an aggregate measure of a quantifiable property of the structure.^{1} By definition, an amorphous structure is (statistically) uniform in all directions, and this corresponds to the circular symmetry of the diffraction pattern. Note in particular that for amorphous structures the diffraction pattern is highly redundant: a single line of data, taken through the centre of the diffraction pattern in any direction, suffices to capture all the information that it embodies and for the nearest-neighbour spacing to be deduced.

One may visualize this argument in the following way: suppose we were to take each pair of points in the vertical direction of the position space image in figure 1 and average them together. This has the effect of halving the Nyquist frequency in the vertical axis direction, but has no effect on the Nyquist frequency of the horizontal axis. Averaging over all points in the vertical axis of the image has the effect of reducing the Nyquist frequency in that direction to zero. The remaining one-dimensional signal contains all the information of interest, specifically it gives information about the distance of peaks in the Fourier coefficients from the origin of *k*-space. We may examine the horizontal axis coefficients and derive all the information required to characterize the structure: the amorphous spatial frequency information is preserved, even if only a very small strip of *k*-space is sampled.^{2} This is illustrated in figure 2.

The practical relevance of the above example is that numerous biological structures are, to a good approximation, either amorphous or quasi-amorphous at least over a local area or volume. Therefore, it is reasonable to suppose that such structures may be characterized adequately using a one-dimensional technique, such as FSA. Furthermore, the spacing of such structures is often of clinical interest: in trabecular bone for example, the strength of the bone is reduced as the spacing of the trabeculae becomes larger. In addition, FSA can be used to derive structural information about crystalline substances so long as the orientation in which the data are acquired is carefully chosen. This is illustrated in §5.2. We have developed a fast MRI method that provides structural information (such as dominant spacing) at high spatial resolutions, and we have applied it to structures that are either amorphous or essentially one-dimensional. We make no claims that the method will give clinically important information in the general case, rather in these two special but frequently occurring cases. In §3, we summarize two methods for computing the spatial frequency information and then illustrate its performance on an increasingly complex set of cases.

## 3. Implementation on magnetic resonance systems

A novel MR pulse sequence, illustrated in figure 3, that acquires a one-dimensional spatial frequency signal from a small ‘prism’ of material has been used to gather data presented in this work.

The pulse sequence consists of 90° and 180° pulses which, when applied as shown in the timing diagram, excites an elongated volume of material whose cross section is square or rectangular. This volume corresponds to the intersection of the orthogonal planar slabs. A readout gradient is then applied along the long axis of the prism, and the resulting echo is recorded. The (one-dimensional) Fourier transform applied to this echo yields the one-dimensional signal magnitude *p*(*x*) as a function of position *x* along the length of the prism. The recorded signal is a one-dimensional representation of the signal magnitude along the length of the prism, and the signal magnitude averages the information over the prism cross section. The length of the prism is frequency encoded, with a gradient sufficient to distinguish the requested number of spatial points over the user-defined length of the prism. The region outside the length of the prism is also encoded, but is filtered out by the scanner in the same way as the region outside a requested image volume would be. The reduction from three or two-dimensions to one-dimension reduces the signal-to-noise ratio of a single measurement significantly, but now each measurement is gathered in a single *T _{R}*. This allows one to gather many more repeated measurements in the same time as one would use to gather a two-dimensional spin echo image, all of which are statistically independent. In this work, the samples are static, so motion is not problematic, but the same sequence has been used to gather clinical data in the spine and liver, the latter of which is limited to a single breath hold. The choice of the size of prism cross section and the number of spatial points in the frequency encode direction are both application-dependent; typically, the signal-to-noise ratio is improved by having a prism as large as possible while still fitting in the region of interest, and as few points as possible while still resolving the structures of interest. It is important to note that any tissue variation in directions other than the frequency encode direction will not be detected by the technique.

Note that a key requirement of the method is elimination of the free induction decay (FID). In practice, the 180° pulse is imperfectly prescribed by the scanner, and any material that falls within this region will contribute an FID-like artefact to the acquired data. This artefact can be minimized by using spoilers to suppress unwanted signals, and using phase cycling that reduces the impact of FID on measured refocused echoes [17]. In practice, phase cycling is very effective for this pulse sequence, as each measurement takes a short time to acquire, so many independent measurements can be made during a scan. Additionally, we are interested in the spatial frequency representation of the data, rather than a spatial image: the recorded echoes are approximately Hermitian functions which means that we may choose to examine only measured positive spatial frequencies, because negative spatial frequencies are more likely to be contaminated by this artefact. However, these components, recorded last in time, are subject to a greater degree of *T*_{2} decay which can be mitigated by not recording a complete set of negative spatial frequencies, shortening the echo time. Below, we quantify this offset using a partial Fourier fraction; the ratio of the length of the echo that is acquired to the length of the full echo containing a complete set of positive and negative spatial frequencies. In order to reconstruct the spatial signal in one dimension, we use a basic homodyne method, because we are not using the reconstructed negative spatial frequencies to produce spectra. For a review of this topic, see [18].

It is possible to acquire multiple prisms in a single *T*_{R}; the only constraint is that subsequent prisms do not lie on the plane of either the 90° or 180° pulses of any of the other prisms being acquired. We achieve this by arranging the prisms in a ‘staircase’ arrangement, ensuring that the order of prism acquisitions is interleaved such that neighbouring prisms are acquired neither immediately before nor immediately after each other, so mitigating the effect the aforementioned imperfect pulse prescription would have on signal level in adjacent prisms (figure 4). It may be advantageous to acquire multiple interleaved prisms across a tissue of interest as this enables a more representative overall characterization of the structures to be performed or the variation of the structures across the tissue to be assessed with no increase in scan time. This pulse sequence has been implemented on a number of Varian and Bruker preclinical scanners and Siemens clinical scanners, and tested in a number of clinical and preclinical applications including the acquisition of data in [19]. For the remainder of this work, we will refer to the pulse sequence we have described in this section as the *custom intersecting sequence*, whenever a distinction between it and another sequence is required.

The pulse sequence described above is by no means the only way to acquire ‘prism’ data. For example, one could use a standard two-dimensional spin echo sequence with a large number of frequency encode points and a small number of phase encode points. The gathered data in this case are in two-dimensional *k*-space, but it is straightforward to transform the two-dimensional data to a set of one-dimensional images, with the number of one-dimensional profiles given by the number of phase encodes used. Many MRI scanners support this type of data acquisition without modification; indeed, we have gathered data of this type using Philips, GE and Siemens clinical systems and measurements made on these different platforms have a coefficient of variation of 5.3%. The similarity between the sequence described above and a standard two-dimensional spin echo means that data from this one-dimensional sequence should be no more subject to the effects of patient motion than conventional imaging, indeed, because a two-dimensional spin echo requires the acquisition of a full set of phase encodes, gathering data with the custom intersecting sequence will make the acquired data more robust to the effects of motion, because each repetition is independent and can be gathered in a single *T*_{R}. For example, involuntary limb motion that causes image blurring, and thus erroneous quantification of trabecular distances during a 10–15 min acquisition is considered virtually impossible to prevent, even with limb restraints [20]. Note that in addition to the advantages discussed above the intersecting pulse sequence has a number of other advantages over a two-dimensional spin echo acquisition of FSA data. For example, there would be much more flexibility in setting the number and position of required prism elements, uncertainty can be estimated directly and the excited volume of the prism is generally much more clearly identified, because this volume is defined by RF pulses rather than a RF pulse and phase encodes.

## 4. Interpretation of data from a one-dimensional excitation

Because one-dimensional spatial ‘images’ cannot be interpreted in the same way as two-dimensional images, we choose to analyse the spatial signals gathered using spectroscopy. Structural spectra are generated from the spatial signal; structural spectra contain information about the presence of structures of specific sizes within the sample.

Because, in a clinical situation, it is unlikely the prism will lie entirely within a region of amorphous tissue, we may limit the spatial range of the data used to generate the structural spectra: this is done using a windowed Fourier transform.^{3} The choice of window function has been a subject of much research in past years [21]. The Blackman–Harris [22] window was chosen, because the frequency response far away from the central peak is substantially below the expected noise level of the MR signal, and the loss-of-frequency resolution is relatively small. The size of the window may have a significant impact on the spectra produced. In practice, over a wide range of examples, we have found a window size that is approximately 10 times the size of the expected feature size provides a good compromise between spatial and frequency resolution. Of course, the Fourier transform of a signal is generally complex valued, so we compute the absolute value of the spectra for presentation. Phase measurements are of interest, not least for signal-to-noise analysis, hence the identification of statistically significant peaks in the spatial frequency data. However, this will be discussed in a subsequent study. The resulting spectra have several features, but the most important clinically are their peaks; the Fourier transform of a periodic function is a sum of Dirac delta functions and hence, when one considers a windowed Fourier transform these delta functions will be broadened by the convolution with the window function in *k*-space. The spatial frequency at which the peak is located provides the inverse of the period of the repeating structure. We emphasize the spatial frequency of the peak represents the period of the repeating structure, or the pitch, and not the size of either the elements making up the repeating structure nor the gaps between them which we now demonstrate: the function known as a pulse train is a series of pulses at regular intervals. Suppose the pulse duration is *A*, whereas the time between pulses is *B*, so the period is *T* = *A* + *B*, the frequency is *f*_{0} = 1/*T* and the angular frequency is . The amplitude of the pulses is *M*. The first few terms of the Fourier transform of the pulse train are
4.1;where
4.2;and asterisk denotes complex conjugation. The frequency locations of the delta function peaks depends only on the period of the repeating structure; the sum of the structure size *A* and the size of the gaps between the structures *B*. It does not depend on either *A* or *B* individually.

The coefficients of the Dirac delta functions depend on the amplitude of the periodic signal, i.e. the contrast. As long as there is sufficient contrast, so that a peak is visible above the chosen noise cut-off (described in the following paragraph) minor changes to the contrast will have little effect on the spectra and no effect on the wavelength of the peak.

As usual, if the noise has zero mean as is the case here, multiple measurements may be made and averaged together in order to improve the signal-to-noise ratio. In addition, the noise can be quantified using these multiple measurements, so that an estimate of uncertainty can be calculated and used to ascribe confidence intervals to the generated spectra. The measured echoes can be separated into a signal and noise component. Both are complex, and the noise is additive white Gaussian which is circularly symmetric in the complex plane. The variance in any direction is *σ*^{2} [23]. Given we are calculating the absolute value of a complex number, the expected distribution will be a Rayleigh distribution if only the noise component is non-zero, and a Rician distribution when there is a non-zero signal present. To compute the confidence intervals, one must invert the cumulative distribution function, which cannot be done for the general case of a Rician distribution. We note that as the signal-to-noise ratio increases, the Rician distribution increasingly approximates a Gaussian distribution with variance *σ*^{2}, and the peaks of interest on structural spectra are all in the high signal to noise regime. Therefore, we assume the calculated sample variance is approximately equal to the Gaussian variance *σ*^{2}.
4.3;where *S*_{i} are the individual measurements and is the sample mean of the set of measurements. Confidence intervals can be calculated using the inverse of the cumulative distribution function of the Rayleigh distribution
4.4;Evaluating this formula gives the magnitude *M* above which we have at least *C* confidence that a peak does not derive from noise. Conventionally, a *p*-value of 0.05 (or a confidence of 0.95) is taken to be the threshold for discarding the null hypothesis.

## 5. Experiments

### 5.1. Simulated circle phantom

We first present a mathematical simulation to illustrate the above. For each sample, a set of circles were placed sequentially such that the arrangement was as closely packed as possible. It is well known that if the circles were all of equal size, the arrangement would be crystalline. However, if some randomness is introduced into the sizes of the circles, then the structure becomes amorphous. The diameters of the circles used followed a Gaussian distribution with mean 0.5 mm and a standard deviation of 0.1 mm. A dot of diameter 0.1 mm was placed at the centre of each circle. An example arrangement of dots is shown in figure 5. To create the one-dimensional signal that would be analysed as described in §5, we averaged over the short axis of each configuration of circles. Forty different configurations of circles were generated randomly, with circle sizes drawn from the same distribution.

#### 5.1.1. Simulated circle phantom results

The spectrum derived from the example arrangement of dots shown in figure 5 is presented in figure 6. A strong peak indicating repeating structure as described in §4 at the average pitch of the dots is evident. Spectra were generated for each of the 40 circle configurations and the highest peak was found in between the wavelength cut-offs of 0 and 1.5 mm. The wavelengths of these peaks are shown in a histogram in figure 7. The mean and standard deviation of the highest peak wavelength found was 0.47 mm and 0.03 mm, respectively, hence, the mean measured pitch was within one standard deviation of the actual mean pitch of the dot structure.

### 5.2. Laser-etched line phantom

Next, we describe an MR scan of a physical phantom comprising laser-etched lines on a Perspex sheet. This experiment will show FSA can be used to deduce a pitch in a non-amorphous (essentially one-dimensional) phantom, as well as to provide some comparison between the two implementations of FSA outlined in §3. The pitch of the lines scanned in this experiment was measured at 497.7 ± 6.2 µm with 2*σ* confidence. The whole phantom was submerged in a solution of copper sulfate and sodium chloride to mimic the MR parameters of biological tissue (*T*_{1} ≈ 200 ms, *T*_{2} ≈ 95 ms) [24]. If the prism is positioned with the long axis orthogonal to the lines the technique can be used to deduce the line spacing; FSA is not sufficient to fully characterize the three-dimensional structure by gathering data in one orientation only. The relevance of a phantom of this type is validation—there is a very clear repeating structure that should lead to a very strong spectral peak at a size indicating the repeat distance of the structure. Data were gathered at the University of Swansea Institute of Life Sciences using a 3 T Siemens Skyra system (see table 1 for the scanning parameters used). The average spacing in the sample found by the scan is presented in §5.2.1 (figure 8). Furthermore, we compare the results derived from scanning the phantom against the custom intersecting sequence with data gathered using a more conventional two-dimensional spin echo imaging sequence with a small number of phase encodes and a larger number of frequency encodes. Data were gathered during the same scanning session so the geometry of the measuring coils and the phantom would not change. The parameters used in the two-dimensional spin echo data acquisition are shown in table 2. The equivalent of 32 prisms were gathered in the phase encode direction, whereas there were 512 points in the frequency encode direction. We acknowledge the differences in the scanning parameters; the greatest effect on the signal-to-noise ratio will come from the different voxel sizes and scan time, which will improve the signal of the two-dimensional spin echo data relative to the custom intersecting sequence, and the bandwidth, which will be detrimental to the signal-to-noise ratio of the two-dimensional spin echo data relative to the custom intersecting sequence. With these two effects, we estimate the signal-to-noise ratio of the two-dimensional spin echo sequence is underestimated by approximately 8%. A more thorough and robust investigation to compare various methods for generating FSA data is left for future work.

#### 5.2.1. Laser-etched line phantom results

Results from this experiment are shown in figure 8. We have set the upper wavelength cut-off of the plot to 1.25 mm as this is well above the expected size of the structure (0.5 mm). There is a very prominent peak close to 500 µm on both spectra, indicating the FSA method has detected a repeating structure in the sample with a pitch close to 500 µm using both data acquisition methods, as expected.

A noise estimate was derived for both data acquisitions and was used to normalize the spectra so that the noise level lies at 1 on the spectral intensity axis, hence, the spectra can be compared directly. Figure 8 shows the custom intersecting sequence has slightly improved the signal-to-noise ratio at the spatial wavelength of interest compared with the spectrum derived from the two-dimensional spin echo data. Furthermore, the two-dimensional spin echo data took 42% longer to acquire than the intersecting sequence data. Despite the caveats mentioned in §5.2, we believe that these results show the custom intersecting sequence is at least as good as the two-dimensional spin echo sequence, but a more thorough comparison between the custom intersecting pulse sequence and existing approaches is left for future work.

### 5.3. Microspheres phantom

Next, we describe the scanning of an amorphous phantom manufactured from microspheres. Four such phantoms were produced from spheres with diameters in the ranges 180–212 µm and 355–425 µm, which were polyethylene microspheres supplied by Cospheric LLC, and diameters in the ranges 500–600 µm and 1180–1400 µm which were glass microspheres supplied by Corpuscular Inc. The two former bead sizes were scanned using a 7 T Varian preclinical MR system at Newcastle University, whereas the latter two bead sizes were scanned on a clinical 3 T Siemens Skyra system at the University of Swansea.

#### 5.3.1. Polyethylene microspheres

The two phantoms made with polyethylene were produced by filling two 1.5 ml microcentrifuge tubes with polydimethylsiloxane (PDMS) liquid and adding the beads, so that they sank under gravity. Specifications for the microspheres state that more than 90% of particles are within the quoted size range, and greater than 90% of particles have an aspect ratio of between 1 and 1.1, i.e. orthogonal diameters differ by a maximum of 10% for 90% of particles. Constructing the phantom in this way ensured the spheres assumed a random close packing arrangement. PDMS was used as it gives an MRI signal and has a significantly lower surface tension than water, reducing the likelihood of pockets or bubbles being formed in the phantom [25]. Prisms were placed along the long axis of the microcentrifuge tubes to ensure the finely sampled axis was as long as possible. Scanning parameters are given in table 3.

#### 5.3.2. Glass microspheres

The two phantoms that comprised glass microspheres were made in a similar way to the polyethylene microspheres phantom described above. Prisms were placed through the centre of the container to maximize ROI length. Scanning parameters are shown in table 4.

#### 5.3.3. Microspheres phantom results

Figure 9 shows the spectra derived from the spatial scans of the microspheres phantom. We set the upper wavelength cut-offs of the plots to be well above the expected pitch of all four sizes of microspheres scanned. It is evident in each case that there are peaks located close to the average bead sizes within each reported range. The peak derived from the scan of the smaller size of glass microsphere is marginally below the minimum size specified by the manufacturer; however, measurement of the bead sizes photographically indicated the large beads were 1.46 ± 0.223 mm in diameter, and the small beads were 0.498 ± 0.093 mm, both with 2*σ* confidence. The spatial wavelengths of the highest peaks on the spectra coincide with these measurements within uncertainties, indicating the repeating pattern of beads in each sample has been detected.

### 5.4. Cadaver bone specimens

In addition to the simulation and physical phantom scans above, we have carried out a comparison between FSA spectra and data from micro-CT (µCT) scans of human cadaver vertebrae. Osteoporosis commonly affects the spine; it is for example one of the locations scanned under the DEXA protocol, therefore, it is valuable to test FSA in this location. Four lumbar (L2–5) and six thoracic (T7–12) vertebrae were extracted from the human cadaver and preserved by freezing. No fixative agent was used, so the tissue contrast would remain as close as possible to that of a living subject. µCT data were acquired in collaboration with the Orthopaedic Bioengineering Research Laboratory (ORBL) at Colorado State University (CSU), using a Scanco µCT 80 scanner. The µCT datasets had isotropic voxels which were 37 µm on each side. MR data were acquired at the University of California, Santa Barbara, using a Siemens TIM Trio 3 T scanner (see table 5 for the MR parameters used).

To compare the two modalities, a structural spectrum was computed for each of the µCT datasets, as described below. Both these one-dimensional spatial signals could then be analysed with the same techniques and compared with evaluate the efficacy of FSA in the analysis of biological structures.

MR analysis of trabecular bone relies on the signal from bone marrow, and the marrow of cadaver bone specimens can contain bubbles of gas caused by air ingress and initial decomposition. Bubbles in the samples would make the trabeculae appear thicker on average in MRI than they would appear in µCT, so it was necessary to extract the gas bubbles before data acquisition. Each sample was thawed and held underwater prior to and during vacuum pumping to remove bubbles. µCT can be used to directly measure bone mineral content, so vacuum pumping of the samples was not strictly required, but was done in advance of both µCT and MR data acquisition for consistency.

Two FSA prisms were acquired for each vertebrae and were positioned mediolaterally in the vertebral body, one central and one anterior with respect to the spinal cord. An example of the placement of the prisms is shown in figure 10. Regions of interest for the prisms were set so as to limit the tissue sampled to trabecular bone. The thoracic vertebrae can be quite small, so we require that the ROI is long enough to accommodate at least one filter window. Out of the 20 samples analysed, all but the anteriorly positioned element of T7 were of sufficient length.

A one-dimensional signal from each of the µCT image sets was derived as follows: data from the same location as the volume of the MR prism were extracted from the full image sets, and rotated such that the axes of the image slices were orthogonal to the prism axes. Only the long axis of the prism, as shown in figure 3 is finely sampled, so the other axes were averaged over. This leaves a one-dimensional spatial signal that can be Fourier transformed and processed in the same way as an echo derived from the MR prism acquisition would be.

Data were analysed with a 19 mm window, and spectra generated were limited to structural wavelengths of between the Nyquist wavelength for the MR data acquisition (0.31 mm) and 2.5 mm as the average structural size in the trabeculae of human vertebrae is expected to be in this range [1]. Spectra were normalized such that
5.1;and compared with the normalized inner product
5.2;where *S*_{C} and *S*_{M} are µCT and MR spectra, respectively.

#### 5.4.1. Cadaver bone specimens results

The main sources of differences between modalities are errors in positioning and the different contrast mechanism, not least because in MRI, signal is derived from the bone marrow, whereas in CT, imaging signal comes from the bone mineral itself. There is no appreciable signal from the bone mineral in MR and the bone marrow in CT. However, because the locations of peaks depends only on the period of the repeating structure, not the spacing or size of the individual elements of the structure as described in §3. The spectra show strong similarity, with a minimum normalized dot product of 0.71. Because MR spectra are generally more noisy than the spectra derived from CT data, the high peaks give a better indication of the structures present in the sample and, furthermore, the peaks on this spectra are indicative of repeating structures which for this application are the most clinically relevant features of the spectrum. For each peak in the CT-derived spectrum, we found the peak in the MR-derived spectrum that was closest, as measured by Euclidean distance in the structural frequency/magnitude plane. The list of MR peak frequencies was compared with the list of CT peak frequencies using Pearson's correlation coefficient that indicated very strong correlation, whose minimum value was 0.992. Histograms of both of these results are shown in figure 11, along with example spectra plots showing both modalities and identified peaks.

## 6. Conclusion and outlook

We have shown that using conventional MRI scanners, and so without ionizing radiation, one can acquire information about amorphous structures, in particular biological structures, which is sufficient to characterize its physical properties using a one-dimensional imaging sequence. Moreover, we have shown the technique presented can characterize repeating structures as small as 180 µm. Furthermore, if one has prior knowledge of a quasi-amorphous or crystalline material or tissue, one can use this same one-dimensional MRI pulse sequence to gather data that will characterize the structures of interest.

This technique could be applied to any pathology where the physical size of structures within organs changes, for example in the osteoporotic degeneration of bone. We have shown above that the one-dimensional MR data acquired for FSA correlates strongly with the equivalent object constructed from high-resolution CT data. Initial studies have been undertaken in a preclinical [19] setting to establish how changes to the structural spectrum reflect changes to the trabecular structure of bone and also how information from these structural spectra can be related to existing measures of bone density. In bone density measurement, the gold standard is μCT and where that technology can be deployed it will produce results that are more informative than the technique presented here. However, because the resolution available for a full three-dimensional scan of *in vivo* bone tissue is often insufficient for characterization of bone microarchitecture for reasons of radiation exposure, the technique we present here may be of interest to clinicians who cannot achieve high enough resolution with conventional imaging to resolve amorphous or quasi-amorphous structures of interest. Further study is required to evaluate the behaviour of the technique *in vivo*.

As remarked in §3, most conventional MRI techniques are applicable to the pulse sequence described in this paper. An exciting potential application would be to acquire data from the brains of people experiencing neural degenerative disorders, perhaps using a pulse sequence incorporating inversion recovery to investigate physical changes during disease progression and better understand this important class of diseases. We expect this technique to be applicable to the problem of understanding brain structure and diseases causing cognitive impairment, because the cortex of the brain consists of highly organized, locally parallel neuronal fibres. These fibres are arranged in bundles which have a pitch of approximately 100–130 µm, and are disrupted as the disease progresses [3]. Such structures are amorphous in the plane perpendicular to the fibres, so placing the prism in this plane should be enable one to characterize the structure present, as long as sufficient contrast can be achieved.

## Authors' contributions

J.R. drafted the manuscript, J.R. and L.F. contributed to the development of the technique and analysis of the data. T.J. invented the technique and supervised the initial data acquisition. D.C. contributed to the development of the technique. J.H. and M.B. supervised the project. All authors consented to the submission of the manuscript.

## Competing interests

We declare that J.R., L.F., J.H. and M.B. are formerly employees and T.J., D.C., J.H. and M.B. are shareholders of Acuitas Medical Ltd, a company that holds intellectual property on the MR pulse sequence implementation and analysis of data acquired by FSA described in the paper.

## Funding

We received no funding for this study.

## Acknowledgements

The authors are grateful to Dominick McIntyre of Cambridge Research Institute for implementing the pulse sequence on the Varian scanner platform and Peter Jezzard and Chris Rodgers of the University of Oxford for implementing the pulse sequence on the Siemens scanner platform, and particularly to Peter Jezzard for very helpful comments on initial drafts. We are indebted to Ross Maxwell, Ian Wilson and Gilberto Almeida of Newcastle University for acquiring the polyethylene bead phantom data, Richard Hugtenberg and Paola Griffiths of Swansea University for acquiring the line phantom and glass bead phantom data and Christian Puttlitz, Snehal Shetye and Kirk McGlivray of Colorado State University and Mario Mendoza of UCSB for acquiring the cadaver bone data. We thank Samantha A. Telfer for providing helpful notes on the analysis of the cadaver vertebrae and comments on initial drafts and Amanda Cox for helpful discussions and comments on initial drafts.

## Footnotes

↵1 Image taken from http://www.caltech.edu/content/caltech-led-team-creates-damage-tolerant-metallic-glass.

↵2 We are grateful to one of the referees for providing this remark.

↵3 In the context of temporal signals, this technique is often called short-time Fourier transform.

- Received July 26, 2016.
- Accepted September 12, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.