Determination of the concentration distribution of sonically dispersive elements

Abstract

A medical apparatus (200, 300, 400, 500) determines the concentration distribution of sonically dispersive elements (606, 2001) within a subject (306, 604, 1004), wherein the medical apparatus comprises: a memory (212) for storing machine executable instructions (224, 226, 228, 230, 232, 318) and a processor (206) for executing the machine executable instructions. Execution of the instructions cause the processor to: receive (100) shear wave data (214) descriptive of the propagation of shear waves (310, 608, 1118) within the subject for at least two frequencies; determine (102) a mechanical property (316, 618, 620) of the subject using the shear wave data at each of the at least two frequencies; determine (104) a power law relationship (218, 702) between the at least two frequencies and the mechanical property; and determine (106) the concentration distribution of the sonically dispersive elements within the subject using the power law relationship and calibration data (222, 704, 800).

Claims

1. A medical apparatus for determining a concentration distribution of sonically dispersive elements within a subject, the medical apparatus comprising: a medical imaging system; a vibration system; a memory for storing machine executable instructions; a processor for executing the machine executable instructions, wherein execution of the instructions cause the processor to: control the vibration system to cause shear waves in the subject; control the medical imaging system to acquire shear wave data using the medical imaging system, wherein the shear wave data is descriptive of a propagation of the shear waves within the subject for at least two frequencies; determine a mechanical property of the subject using the shear wave data; determine a power law relationship between the at least two frequencies and the mechanical property; determine calibration data by modeling scattering of shear waves by the sonically dispersive elements as a function of the shear wave frequency; determine the concentration distribution of the sonically dispersive elements within the subject using the power law relationship and the calibration data; and displaying the concentration distribution of the sonically dispersive elements on a display, wherein the sonically dispersive elements are blood vessels within the subject, and execution of the instructions cause the processor to determine the concentration distribution of the blood vessels within the subject using the power law relationship and the calibration data wherein the calibration data comprises an average size of blood vessels within the subject.

2. The medical apparatus of claim 1, wherein the medical imaging system is an ultrasound system, wherein the ultrasound system is arranged to acquire ultrasound data, wherein the ultrasound system is arranged to determine the shear wave data by tracking speckle patterns in the ultrasound data.

3. The medical apparatus of claim 1, wherein the medical imaging system is a magnetic resonance imaging system, wherein the magnetic resonance imaging system is arranged to acquire magnetic resonance elastography data, wherein the magnetic resonance imaging system is arranged to determine the shear wave data using the magnetic resonance elastography data.

4. The medical apparatus of claim 1, wherein the vibration system comprises one of an ultrasound transducer or a mechanical actuator.

5. The medical apparatus of claim 1, wherein the vibration system comprises a high intensity focused ultrasound system, wherein the high intensity focused ultrasound system is operable for inducing the shear waves using sonic radiation force.

6. The medical apparatus of claim 1, wherein the vibration system is operable to generate shear waves with a frequency of 10 Hz to 1000 Hz.

7. The medical apparatus of claim 1, wherein the shear wave data, the mechanical property, the power law relationship, and the concentration distribution of the sonically dispersive elements have a two-dimensional spatial dependence.

8. The medical apparatus of claim 1, wherein execution of the machine executable instructions further cause the processor to perform an operation selected from the group consisting of storing the concentration distribution of the sonically dispersive elements in the memory, sending the concentration distribution of the sonically dispersive elements to a computer system via a computer network, and combinations thereof.

9. The medical apparatus of claim 1, wherein execution of the machine executable instructions further cause the processor to generate the calibration data by modeling the scattering of shear waves by the sonically dispersive elements as a function of shear wave frequency and sonically dispersive element size.

10. The medical apparatus of claim 1, wherein the mechanical property selected from the group consisting of elasticity, viscosity, propagation, attenuation, and the dispersion relation.

11. The medical apparatus of claim 1, wherein the shear wave data, the mechanical property, the power law relationship, and the concentration distribution of the sonically dispersive elements have a three-dimensional spatial dependence.

12. A non-transitory computer readable medium comprising machine executable instructions for execution by a processor controlling a medical apparatus, the medical apparatus arranged to determine a density distribution of blood vessels within a subject, wherein execution of the instructions cause the processor to: receive shear wave data, wherein the shear wave data describes a propagation of shear waves within the subject for at least two frequencies; determine a mechanical property of the subject using the shear wave data at each of the at least two frequencies; determine a power law relationship between the at least two frequencies and the mechanical property; and determine the density distribution of blood vessels within the subject using the power law relationship and calibration data comprising an average or typical size of blood vessels within the subject; wherein the density distribution of blood vessels within the subject is used to determine information of the blood vessels to diagnosis the patient with a pathology.

13. The non-transitory computer readable medium of claim 12, wherein the shear wave data, the mechanical property, the power law relationship, and the density distribution of blood vessels in within the subject have a three-dimensional spatial dependence.

14. The non-transitory computer readable medium of claim 12, wherein execution of the machine executable instructions further cause the processor to perform at least one operation selected from the group consisting of: storing the density distribution of blood vessels within the subject in a memory, displaying the density distribution of blood vessels within the subject on a display, and sending the density distribution of blood vessels within the subject to a computer system via a computer network.

15. A medical apparatus for determining a density distribution of blood vessels within a subject, the medical apparatus comprising: a medical imaging system; a vibration system; and the non-transitory computer readable medium of claim 12, wherein execution of the instructions further cause the processor to: control the vibration system to cause the shear waves within the subject using the vibration system; and control the medical imaging system to acquire the shear wave data.

16. A method of determining the concentration distribution of sonically dispersive elements within a subject, wherein the method comprises the steps of: receiving shear wave data, wherein the shear wave data describes the propagation of shear waves within the subject for at least two frequencies; determining a mechanical property of the subject using the shear wave data at each of the at least two frequencies, the mechanical property including one of elasticity, propagation, and the dispersion relation; determining a power law relationship between the at least two frequencies and the mechanical property; and determining the concentration distribution of the sonically dispersive elements within the subject using the power law relationship and calibration data; wherein the concentration distribution of the sonically dispersive elements within the subject is used to determine information of the sonically dispersive elements to diagnosis the patient with a pathology; and wherein the sonically dispersive elements are blood vessels within the subject, and determining the concentration distribution of the blood vessels within the subject using the power law relationship and the calibration data wherein the calibration data comprises an average size of blood vessels within the subject.

17. The method of claim 16, wherein the method further comprises the steps of: measuring a calibration power law relationship as a function of the concentration distribution of dispersive elements; and determining the calibration data empirically using the calibration power law relationship.

18. The method of claim 16, wherein the shear wave data, the mechanical property, the power law relationship, and the concentration distribution of the sonically dispersive elements have a three-dimensional spatial dependence.

19. The method of claim 16, further including at least one of: storing the concentration distribution of the sonically dispersive elements in a memory, displaying the concentration distribution of the sonically dispersive elements on a display, and sending the concentration distribution of the sonically dispersive elements to a computer system via a computer network.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) In the following preferred embodiments of the invention will be described, by way of example only, and with reference to the drawings in which:

(2) FIG. 1 shows a flow chart which illustrates a method according to an embodiment of the invention.

(3) FIG. 2 illustrates an example of a medical apparatus;

(4) FIG. 3 illustrates a further example of a medical apparatus;

(5) FIG. 4 illustrates a further example of a medical apparatus;

(6) FIG. 5 illustrates a further example of a medical apparatus;

(7) FIGS. 6 and 7 illustrate a series of images which are used to illustrate a portion of the method of determining the concentration distribution of the sonically dispersive elements;

(8) FIG. 8 illustrates multiple solutions to the calculation of the particulate concentration;

(9) FIG. 9 shows a flow diagram which illustrates how the correct solution can be determined;

(10) FIG. 10 shows a light microscopy image of a colloidal gel specimen with particulates embedded;

(11) FIG. 11 illustrates an experimental setup;

(12) FIG. 12 shows a plot of experimental data;

(13) FIG. 13 shows a further plot of experimental data;

(14) FIG. 14 shows Box-counting results for different densities of a fixed particle size;

(15) FIG. 15 shows an example of the fractal dimension df as a function of concentration;

(16) FIG. 16 shows the characteristic length (in units of pixels=m here) as a function of

(17) ${}^{-\frac{1}{D}}$
with D=2; and

(18) FIG. 17 shows the schematic depiction of the two contributions for the total lag-time distribution a(t).

DETAILED DESCRIPTION OF THE EMBODIMENTS

(19) Like numbered elements in these figures are either equivalent elements or perform the same function. Elements which have been discussed previously will not necessarily be discussed in later figures if the function is equivalent.

(20) FIG. 1 shows a flow chart which illustrates a method according to an embodiment of the invention. First in step 100 shear wave data is received which is descriptive of the propagation of shear waves within the subject for at least two frequencies. Next in step 102 a mechanical property of the subject is determined using the shear wave data at each of the at least two frequencies. Then in step 104 a power law relationship between the at least two frequencies and the mechanical property is determined. Then in step 106 the concentration distribution of the sonically dispersive element within the subject is determined using the power law relationship and calibration data. For instance a lookup table could be created which contains entries for various power law relationships values as a function of the concentration distribution of the sonically dispersive elements.

(21) It should be noted that in some instances it may be beneficial to use a priori knowledge of the size or approximate size of the sonically dispersive elements within the subject. In some cases the size of the sonically dispersive elements may be well known. For instance if the sonically dispersive elements are blood vessels the average or typical size of blood vessels within the subject may be known beforehand and may be useful also in conjunction with the calibration data. For instance the calibration data could be for a particular type or distribution of dispersive elements.

(22) FIG. 2 illustrates an example of a medical apparatus 200. The medical apparatus 200 comprises a computer 202. The computer 202 has a hardware interface 204 connected to a processor 206. The processor 206 is also connected to a user interface 208 and computer storage 210 and computer memory 212. Within the computer storage 210 is stored shear wave data 214. The shear wave data 214 is used to drive a mechanical property 216 which is also stored in the computer storage 210. The computer storage 210 also contains a power law relationship 218 derived or calculated from the mechanical property 216. The computer storage 210 is also shown as containing a concentration of sonically dispersive elements 220 which was calculated using the power law relationship 218 by comparing it to calibration data 222. The calibration data 222 is also shown as being stored in the computer storage 210. The shear wave data 214, the mechanical property 216, the power law relationship 218, and the concentration of the sonically dispersive elements 220 or concentration distributrion of the sonically dispersive elements may have a spatial dependence.

(23) The computer memory 212 is shown as containing a control module 224. The control module 224 enables the processor 206 to control the operation and function of the medical apparatus 200. In the case of additional components being added to the medical apparatus 200 such as a mechanical actuator or the system for generating shear waves or a medical imaging system, the processor 206 may be enabled by the control module 224 to control them via the hardware interface 204. The computer memory 212 is further shown as containing a shear wave data processing module 226. The shear wave data processing module 226 contains computer-executable code which enables the processor 206 to determine the mechanical property 216 from the shear wave data 214.

(24) The computer memory 212 is further shown as containing a power law determination module 228. The power law determination module 228 contains computer-executable code which enables the processor 206 to determine the power law relationship 218 from the mechanical property 216. The computer memory 212 is further shown as containing a concentration determination module 230. The concentration determination module 230 enables the processor 206 to determine the concentration of sonically dispersive elements 220 or the concentration distribution of sonically dispersive elements using the power law relationship 218 and the calibration data 222.

(25) Finally the computer memory 212 is shown as containing a calibration data generation module 232. The calibration data generation module 232 is an optional module which in some embodiments would be used to theoretically calculate the calibration data 222. In other embodiments the calibration data generation module may use empirical measurements to derive or calculate the calibration data 222.

(26) FIG. 3 shows a further example of a medical imaging system 300. The medical imaging system in FIG. 3 is similar to that shown in FIG. 2 except there is additionally a medical imaging system 302 and a mechanical actuator 310 are also shown as being included. The medical imaging system 302 is intended to be representative and may be any medical imaging system which is able to detect shear waves traveling through a subject 306. The medical imaging system in particular may be representative of a magnetic resonance imaging system or an ultrasound system. There is a subject 306 shown as reposing on a subject support 308 partially within an imaging zone 304.

(27) There is a mechanical actuator 310 in contact with the subject 306 which is generating shear waves 314. In some instances there may be a catheter or object inserted into an orifice to locally generate shear waves also. The mechanical actuator 310 is connected to a mechanical actuator controller 312 which supplies electrical power or other actuation for moving the mechanical actuator 310. In the case of a magnetic resonance imaging system the mechanical actuator controller 312 may for instance provide pneumatic power to the mechanical actuator 310 or may move a non-magnetic rod. The medical imaging system 302 and the mechanical actuator controller 312 are shown as being connected to a hardware interface 204. This enables the processor 206 to control the operation and function of the various components of the medical imaging system 300.

(28) The computer storage 210 is shown as containing medical image data 316 that was acquired using the medical imaging system 302. The computer memory 212 is shown as additionally containing an image processor module 318. The imaging processing module 318 enables the processor 206 to generate the shear wave data 214 from the medical image data 316.

(29) FIG. 4 shows a further example of a medical apparatus 400. In this example the medical imaging system is a magnetic resonance imaging system 402. The magnetic resonance imaging system comprises a magnet 404. The magnet 404 is a cylindrical type superconducting magnet with a bore 406 through the center of it.

(30) The magnet has a liquid helium cooled cryostat with superconducting coils. It is also possible to use permanent or resistive magnets. The use of different types of magnets is also possible for instance it is also possible to use both a split cylindrical magnet and a so called open magnet. A split cylindrical magnet is similar to a standard cylindrical magnet, except that the cryostat has been split into two sections to allow access to the iso-plane of the magnet, such magnets may for instance be used in conjunction with charged particle beam therapy. An open magnet has two magnet sections, one above the other with a space in-between that is large enough to receive a subject: the arrangement of the two sections area similar to that of a Helmholtz coil. Open magnets are popular, because the subject is less confined. Inside the cryostat of the cylindrical magnet there is a collection of superconducting coils. Within the bore 406 of the cylindrical magnet there is an imaging zone 408 where the magnetic field is strong and uniform enough to perform magnetic resonance imaging. The mechanical actuator 310 is shown as being in the bore of the magnet 406.

(31) Within the bore 406 of the magnet there is also a set of magnetic field gradient coils 410 which are used for acquisition of magnetic resonance data to spatially encode magnetic spins within the imaging zone 408 of the magnet 404. The magnetic field gradient coils are connected to a magnetic field gradient coil power supply 412. The magnetic field gradient coils 410 are intended to be representative. Typically magnetic field gradient coils contain three separate sets of coils for spatially encoding in three orthogonal spatial directions. A magnetic field gradient power supply 412 supplies current to the magnetic field gradient coils 410. The current supplied to the magnetic field coils is controlled as a function of time and may be ramped or pulsed.

(32) Adjacent to the imaging zone 408 is a radio-frequency coil 414 for manipulating the orientations of magnetic spins within the imaging zone 408 and for receiving radio transmissions from spins also within the imaging zone. The radio-frequency coil may contain multiple coil elements. The radio-frequency coil may also be referred to as a channel or an antenna. The radio-frequency coil 414 is connected to a radio frequency transceiver 416. The radio-frequency coil 414 and radio frequency transceiver 416 may be replaced by separate transmit and receive coils and a separate transmitter and receiver. It is understood that the radio-frequency coil 414 and the radio-frequency transceiver 416 are representative. The radio-frequency coil 414 is intended to also represent a dedicated transmit antenna and a dedicated receive antenna. Likewise the transceiver 416 may also represent a separate transmitter and receivers.

(33) The mechanical actuator controller 312, the transceiver 416 and the magnetic field gradient coil power supply 412 are shown as being connected to the hardware interface 204 of the computer 202.

(34) In this example the medial image data is magnetic resonance data. The computer storage 210 is further shown as containing a pulse sequence 420. The pulse sequence 420 is a set of commands or information which may be used to derive a set of commands for controlling the magnetic resonance imaging system 402 to acquire the magnetic resonance data 316. For instance the control module 224 could use the pulse sequence 420 to acquire the magnetic resonance data 316.

(35) FIG. 5 shows a medical instrument 500 that is similar to the embodiment shown in FIG. 4. However, in FIG. 5 a high-intensity focused ultrasound system 506 is used instead of the mechanical actuator. The high-intensity focused ultrasound system 506 focuses the ultrasound to a point 522. By switching the ultrasound on or off or modulating it shear waves can be generated within the subject 306. The shear waves 314 can be shown as radiating outwards from the focal point 522.

(36) A subject 306 is shown as reposing on a subject support 308. The medical apparatus 200 comprises a high-intensity focused ultrasound system 506. The high-intensity focused ultrasound system comprises 506 a fluid-filled chamber 508. Within the fluid-filled chamber 508 is an ultrasound transducer 510. Although it is not shown in this figure the ultrasound transducer 510 comprises multiple ultrasound transducer elements each capable of generating an individual beam of ultrasound. This may be used to steer the location of a focal point 522 electronically by controlling the phase and/or amplitude of alternating electrical current supplied to each of or groups of the ultrasound transducer elements. Point 522 represents the adjustable focus of the medical apparatus 500.

(37) The ultrasound transducer 510 is connected to a mechanism 512 which allows the ultrasound transducer 510 to be repositioned mechanically. The mechanism 512 is connected to a mechanical actuator 514 which is adapted for actuating the mechanism 512. The mechanical actuator 512 also represents a power supply for supplying electrical power to the ultrasound transducer 510. In some embodiments the power supply may control the phase and/or amplitude of electrical power to individual ultrasound transducer elements.

(38) The ultrasound transducer 510 generates ultrasound which is shown as following the path 516. The ultrasound 516 goes through the fluid-filled chamber 508 and through an ultrasound window 518. In this embodiment the ultrasound then passes through a gel pad 520. The gel pad 520 is not necessarily present in all embodiments but in this embodiment there is a recess in the subject support 308 for receiving a gel pad 520. The gel pad 520 helps couple ultrasonic power between the transducer 510 and the subject 306. After passing through the gel pad 520 the ultrasound 516 is focused to a sonication volume 522 or target zone.

(39) The sonication volume 522 may be moved through a combination of mechanically positioning the ultrasonic transducer 510 and electronically steering the position of the sonication volume 522. By modulating or pulsing the intensity of ultrasound focused at the focal point 522 shear waves 314 can be induced in the subject.

(40) FIG. 6 shows a series of images which are used to illustrate a portion of the method of determining the concentration distribution of the sonically dispersive elements. First image 600 shows a transducer for mechanical actuator 602 in contact with a matrix 604 filled with a variety of particulates 606. The ultrasound transducer 602 is able to induce a shear wave which is partially dispersed by the particulates 606. Next image 608 shows a shear wave image 608 or shear wave data that is descriptive of the transport of the shear wave through the matrix 604. A shear wave imaging system permits to acquire or register shear wave propagation.

(41) Image 610 represents reconstructed images used for determining the mechanical properties. Dedicated software allows the reconstruction of mechanical properties from the theory of wave propagation. This may include but is not limited to the elasticity, viscosity, propagation, attenuation and the dispersion relation of the waves. Next image 613 shows two plots of the logarithm of the frequency 614 versus the logarithm of the mechanical property 616. A local volume 612 is indicated in the image 610. The local volume 612 is examined and the mechanical property at that particular point for this frequency is determined. This is plotted value 618 on the first plot. This is then repeated at multiple frequencies and the multiple plot values are shown as 620. The experiment is repeated at different frequencies in order to obtain a frequency dependence of the mechanical property.

(42) Next, FIG. 7 shows a portion of image 613 again. The multiple values 620 are plotted and a power law fit indicated by the dashed line is performed. The value of alpha is the slope of this line 702 and represents the power law fit. Extraction of the power law exponent alpha 702 characterizes the frequency dependence of the mechanical property. The image 704 illustrates one theoretical model which may be used for interpreting the particulate density in terms of alpha 702. This plot shows the spatial scale versus a filling factor calculated for a theoretical representation of particulates dispersed within a matrix. A box counting algorithm is used to characterize the volume or area being examined. On the spatial scale there is a value 710 which characterizes a transition between a normal fractal filling regime 712 and a normal Euclidian filling regime 714. The distribution of obstacles which may be for instance blood vessels or particles in space can be analyzed as a function of the spatial scale that is whether we look at the micrometer, the millimeter or the centimeter scale. It can be shown for example in image 704 that the filling space where obstacles changes between a classical geometrical filling and a fractal filling space for a given spatial scale 710. The anomalous part is linked to the exponent or the slope of alpha. This part is defined by the micro-architecture of the material, hence it is possible to deduce the details of the microarchitecture from the value of alpha 702, This is only feasible when the wavelength of the shear wave is sufficiently small to sense the fractal filling regime 712.

(43) FIG. 8 illustrates that there may be multiple solutions to the particulate density. Image 613 is used to represent calculating alpha again. In plot 800 the concentration distribution versus the value of alpha is shown. There may be two solutions 802 for a particular measured value of alpha. This value was derived using theoretical derivation of the dispersion properties of the shear wave:

(44) $\mathrm{Disp}={}_F^{\mathrm{df}-d}\left(df-d,\left(\frac{1}{{}_F}i\right)t_0\right),$
where Disp is measure if the measurement is within the fractal or Euclidian regime, df is the fractal dimension, d is the dimension, is the frequency being investigated, and t.sub.0 is a characteristic time which is equivalent to the radius of the particulates. When Disp is much greater than one then there is a fractal effect, and when Disp is much less than one there no fractal effect. A brief outline of the theory of this above equation is contained in the following appendix. In particular see the portion of Eq. 9 labeled general weight in the appendix.

(45) For a given measurement the concentration or the distribution of concentration may have more than one solution. To get the concentration from the experimental data the use of a knowledge of the macroscopic pathology may be used to determine the radius and a priori information may be used to distinguish between the two possible solutions. For instance it may be known what sort of blood vessels or particulates are inside of the subject. This would allow elimination of one of the possible solutions.

(46) FIG. 9 shows a method of determining the correct solution. In step 900 data acquisition and extraction of the value of .sub.mp is determined, wherein .sub.mp is the power law relationship for the mechanical property mp. Next in step 902 an estimation of the obstacle size and concentration or concentration distribution according to an observed pathology is determined. Next in step 904 the injection of rO.sub.E, the size range of the particulates, and C.sub.E, the concentration range, in the fractal effect formula (Disp, shown above) is performed. Choosing rO.sub.E and C.sub.E imparts some .sub.a priori knowledge of the microstructure and enables determining the solution.

(47) If the value of Disp is much less than 1 then branch 906 is selected. In this case there is no fractal effect 908 and it is not possible to extract the concentration 910. In this case the frequency is changed 914 and the method returns to step 900. In case the value Disp is much greater than 1 912 then there is a fractal effect 924. If C.sub.E is much less than C.sub.min (the lower minimum concentration solution) then the low concentration 918 is selected. If the concentration is at the minimum 920 then there is only one solution and the solution is known. If C.sub.E is greater than C.sub.min then the high concentration is selected 926. The use of equation Disp is not necessary. Experiments or numerical simulations could be performed to determine the relation between .sub.mp and the concentration.

(48) Magnetic Resonance Elastography (MRE) is a technique capable of noninvasively assessing the mechanical properties of tissues. The assessment of these properties is done indirectly via the measurement of low frequency mechanical shear waves traversing the tissue. It can be hypothesized that the presence of micro-obstaclessimilar to effects leading to the apparent diffusion coefficientchanges the dispersion relation of propagating shear waves and hence might influence at the macroscopic scale the apparent mechanical properties of the medium. In diffusion weighted imaging (DWI), disordered media can lead to two effects: reduction of the typical diffusion length leading to the apparent diffusion coefficient and/or a mean-square displacement which is not anymore proportional to time but to a fractional power of time not equal to one (so-called anomalous diffusion). In DWI, micro-structural information is lost due to the massive averaging that occurs within the imaging voxel and can only be revealed when exploring the tissue using different b-values. Similarly here, where the propagation of a mechanical wave enters into the diffusive regime due to multiple scattering effects, the frequency dependence of the mechanical properties could allow the assessment of the sub-voxel microarchitecture. In this study we investigate the propagation of shear waves in calibrated phantoms containing accurately controlled size distributions of scattering particles and demonstrate for the first time that shear waves are able to reveal at the macroscopic scale the hidden microarchitecture properties of the material.

(49) To test this experimental, gel phantoms were fabricated using an agarose solution at 15 g/L (BRL, Type 5510UB) prepared in a water bath at 80 C. In order to create well defined scattering particle size distributions, colloidal suspensions of polystyrene microspheres with precisely known diameter (1 m, 5 m, 10 m, 30 m and 150 m diameter, Sigma-Aldrich) and concentrations were added to the gel before solidification. This is shown in FIG. 10.

(50) FIG. 10 shows a light microscopy image 1000 of a colloidal gel specimen 1004 with particulates embedded 1002. The image was taken at a magnification of 50 using a Leica microscope. Clearly, particles of different sizes can be identified. The thereby measured diameter distribution per volume corresponds to the expected theoretical value hence validating the desired microarchitectural properties of the gel.

(51) The aim was to maintain for all prepared gels a concentration of 8% of spheres (similar to the volume fraction of blood vessels in tissue). The polystyrene microspheres have an extremely elevated shear modulus (MPa) and hence can serve as microscopic scatterers in the soft gel (kPa). Different sample were prepared: gels without spheres serving as reference, gels with only one type of spheres (so-called monosize gel) and gels with particle size distributions which followed a power law and hence possessed fractal properties. Different exponents of power-law particle size distributions (#d.sup.Y, with d the particle diameter) were fabricated (=2, 1, 0). A -value of zero indicates a flat distribution meaning that as many small as large particles are present. MRE was performed on a horizontal 7 T imaging scanner (Pharmascan, Bruker, Erlangen, Germany). Mechanical vibrations were generated by a toothpick placed in the center of the sample to induce a circular propagation. An electromagnetic shaker located outside the MR scanner was used to transmit mechanical vibrations via a flexible carbon fiber rod to the toothpick. This is shown in FIG. 11.

(52) FIG. 11 shows a schematic description of the experimental setup used. The gel is filled into an insert which is mounted onto the MRE setup. FIG. 11 shows an example of the experimental setup. There is a electromagnetic shaker 1100 which is connected to a carbon rod 1102. The electromagnetic shaker 1100 causes the carbon rod 1102 to move in the direction indicated by the arrows 1104. The carbon rod 1102 is connected to a cradle 1106. The cradle 1106 translates the motion of the rod 1104 into a different motion indicated by the arrows 1108. 1108 is transverse to 1104. A toothpick 1110 is mounted in the cradle 1106. The toothpick 1110 is inserted into a container 1112 that is filled with a gel 1114 and sealed with parafilm 1116. The toothpick 1110 vibrates up and down inducing shear waves 1118 in the gel 1114.

(53) Samples placed around the toothpick 1110 were always at the same height via a home-made support. A surface receiver coil was placed around the sample at the level of the gel to assure optimal signal-to-noise. For each phantom a steady-state MRE sequence was applied with a mechanical excitation frequency in the range of 150 to 300 Hz and the following sequence parameter: 8 dynamics, 7 contiguous transverse slices with slice thickness of 0.4 mm, field of view=25 mm25 mm, matrix size=256256, TE/TR=2717/427353 ms and acquisition time in the range of 6 to 10 min depending on the excitation frequency and on the number of motion encoding gradient periods. The MRE sequence was acquired for the three spatial direction of motion in order to obtain volumetric images of the 3D propagating mechanical wave inside the phantom. In order to take into account a potential temporal evolution of the gel during the entire acquisition time (up to 300 mins!), the first experiment was repeated at the end of the acquisition time. This allowed correcting for potentially drying effect. Data was reconstructed with an isotropic reconstruction technique.

(54) In examining the experimental results, the complex-shear modulus (G*) of each phantom increased by a maximum of 10% between the beginning and the end of the multifrequency-MRE experiment due to aging effects. As presented in FIG. 2, results show that the macroscopic shear modulus is frequency-dependant for the four investigated samples and follows a power law with |G*()|=.Math..sup.z0. The power coefficient z.sub.0 of a gel with the 10 m-monosize distribution of microspheres is almost unchanged as compared to z.sub.0 of the reference gel, shown in FIGS. 12 and 13.

(55) FIG. 12 shows a plot of experimental data. FIG. 12 is a plot of the frequency in Hertz 1200 versus the normalized complex shear modulus 1202. The + marked points 1204 are measurements for gel with 10 m. The line 1206 is a power law fit to the data 1204. The points marked with an x 1208 are taken for the reference gel. The line 1210 is a power law fit to 1208.

(56) FIG. 13 shows more experimental data. The points labeled 1304 or a + correspond to a fractal of a gamma=0. The line 1306 is a power law fit to the data 1304. The data marked with an x 1308 is the reference gel. The line 1310 is a power law fit to the data 1308.

(57) However, in the presence of a fractal distribution of microspheres, z.sub.0 increases significantly compared to the reference gel by a factor of 2.2. All other fractal gels demonstrated equally a significant increase in z.sub.0.

(58) The experimental tests demonstrate that the frequency-dependence of mechanical shear wave diffusion can allow probing sub-voxel distributions of scattering structures and as a consequence overcome the spatial resolution limitation relying intrinsically on the MR imaging sensitivity. These experimental results have been theoretically and numerically via FEM simulations confirmed (not shown). However, in this study mechanical properties of the gel were critically relying on the fabrication process and only relative slopes of different gels have been compared. The solidification process of the colloidal gels must be improved and additional imaging modalities should be involved such as CT-scans in order to image the microspheres distribution in phantoms after solidification of the gel that probably induces microspheres aggregation into fractal flocs. Moreover, the studied gels consisted of very simplified biphasic structural arrangements with particles being about 1000 times stiffer than the background gel. Biological tissue represents a far more complex arrangement with variations not only in size, but also in stiffness contrast and length distribution. Phantoms with microspheres exhibiting multi-size distributions and multiple elasticity properties would be better to simulate real tissue. The here observed effect might play an important role in understanding the influence of microscopic tissue components on mechanical properties as measured by elastography techniques. It opens the perspective of detecting and describing micro-inclusions, such as small metastases or neo-vascularisation, from elastography data, which are not directly detectable by MRE.

(59) While the invention has been illustrated and described in detail in the drawings and foregoing description, such illustration and description are to be considered illustrative or exemplary and not restrictive; the invention is not limited to the disclosed embodiments.

(60) Other variations to the disclosed embodiments can be understood and effected by those skilled in the art in practicing the claimed invention, from a study of the drawings, the disclosure, and the appended claims. In the claims, the word comprising does not exclude other elements or steps, and the indefinite article a or an does not exclude a plurality. A single processor or other unit may fulfill the functions of several items recited in the claims. The mere fact that certain measures are recited in mutually different dependent claims does not indicate that a combination of these measured cannot be used to advantage. A computer program may be stored/distributed on a suitable medium, such as an optical storage medium or a solid-state medium supplied together with or as part of other hardware, but may also be distributed in other forms, such as via the Internet or other wired or wireless telecommunication systems. Any reference signs in the claims should not be construed as limiting the scope.

APPENDIX

(61) Geometrical Characterization of the Material Via Box-Counting: Pair-Correlation Function

(62) We will investigate how to characterize a simple homogeneous elastic medium (no viscosity) which is filled with very stiff particles of a fixed size. For that purpose we will use the box-counting method as shown in FIG. 14. FIG. 14 shows Box-counting results for different densities and a fixed particle diameter of 10 m.

(63) This diameter corresponds to 1=log(10) on the x-axis.

(64) We can identify two distinct regions for this type of composite material which separate at the characteristic length such that
N(r)r.sup.df<(1)
N(r)r.sup.D<(2)
where we have introduced df as fractal dimension. This can be considered here as its definition and as such df represents the power-law exponent within a certain bandwidth.

(65) For low densities we observe df1 while for higher densities df approaches 2 for those 2D simulations (see FIG. 15). FIG. 15 shows the fractal dimension df (i.e. slope of non-euclidean part in FIG. 14) as a function of concentration.

(66) This dependence is not derived from first principles. It is merely an experimental result. However, since the material does not form any complex aggregates, the characteristic length must follow under these conditions the simple geometric relationship

(67) $\begin{array}{cc}{}^{-\frac{1}{D}}& \left(3\right)\end{array}$
with the density of particles submersed in the homogeneous background. This relationship is shown in FIG. 16.

(68) In order to evaluate the probability of finding at distance r a particle (i.e. we are searching the so-called pair-correlation function), we follow the approach of Teixeira, see Teixeira, J. Small-angle scattering by fractal systems. J. Appl. Cryst 21, 781-785 (1988). The number of particles within the radius r from the origin can be written as
N(r)=.sub.0.sup.rdrP(r).Math.(2r).sup.D-1.(4)

(69) Differentiation of Eqs. 2 and 4 leads to the following expression for the probability density function:
P(r)r.sup.dfD(5)

(70) Apparently, in case of an euclidean distribution with dfD, we find P(r)r.sup.0 which leads

(71) to N(r)r.sup.0 as expected. If we consider the background of the material as constant and isotropic,

(72) this probability density represents the lag-time distribution for a material since it describes how likely it is to find at distance r another obstacle. That is the key idea in order to connect this fractal part with the ODA theory. FIG. 16 shows the characteristic length (in units of pixels=m here) as a function of

(73) ${}^{-\frac{1}{D}}$
with D=2 here.

(74) The results of FIGS. 14 to 16 hence provide the following parameterizations
df 1/15(p5)+1(6)
155[m].Math..sup.1/2(7)
with the concentration in %.

(75) With this knowledge we can now construct a lag-time distribution which allows to analytically solve the necessary equations in ODA to calculate the dispersion properties of the propagation . Hence, the lag-time distribution a(t) is composed of two terms: one term describing the fractal part and one the classical euclidean part. The limit of validity of the fractal part is given by the characteristic length which is called in the temporal domain .sub.F=/c.sub.0 with c.sub.0 the speed of the wave in the background material. In order to render the equations analytical we use exponential functions for the suppression. The Euclidean part is accordingly suppressed for small distances by 1-e.sup.t/F. In order to prevent lag-times of infinite value, the euclidean part is furthermore suppressed by an exponential function with the characteristic time constant .sub.D with .sub.F<.sub.D This yields the following lag-time distribution (see FIG. 17):
a(t)=t.sup.dfD1e.sup.t/.sup.F+t.sup.1(1e.sup.t/.sup.P)e.sup.t/.sup.D,(8)
where an addition 1/r has been introduced since we want to use the probability density as developed for the dimensionality D for the ODA theory which operates in 1D!

(76) FIG. 17 shows the schematic depiction of the two contributions for the total lag-time distribution a(t). The finite particle size limits the analysis to t>t.sub.0=r.sub.0/c.sub.0.

(77) The translation from lag-time distribution to dispersion relation for necessitates to calculate the Fourier sinus transform of a(l), i.e. we need to calculate the characteristic equation (see Gradshteyn, I. S. & Ryzhik, I. M. Table of Integrals, Series, and Products

(78) (Academic Press, Burlington, Mass., 2007), 7th edn. p.498/eq.2):

(79) $\begin{array}{cc}2.{}_{}^{}{x}^{-1}{e}^{-x}\sin x\mathrm{dx}=\frac{2}{}{\left(+i\right)}^{-}\left[,\left(+i\right)u\right]-\frac{2}{}{\left(-i\right)}^{-}\left[,\left(-i\right)\right]\left[\mathrm{Re}>.Math.\mathrm{Im}.Math.\right]& \mathrm{ET}|318\left(9\right)\\ {}_{t_0}^{}\mathrm{dt}{t}^{\left(\mathrm{df}-D\right)-1}{e}^{-t/{}_F}\sin \left(t\right)=\frac{i}{2}{\left(\frac{1}{{}_F}+i\right)}^{D-\mathrm{df}}\left(\mathrm{df}-D,\left(\frac{1}{{}_F}+i\right)t_0\right)-\frac{i}{2}{\left(\frac{1}{{}_F}-i\right)}^{D-\mathrm{df}}\left(\mathrm{df}-D,\left(\frac{1}{{}_F}-i\right)t_0\right)={}_F^{\left(\mathrm{df}-D\right)}\left(\mathrm{df}-D,\left(\frac{1}{{}_F}+i\right)t_0\right).Math.\frac{i}{2}\left[{\left(1+i{}_F\right)}^{D-\mathrm{df}}-{\left(1-i{}_F\right)}^{D-\mathrm{df}}\right]=\underset{\underset{\mathrm{general}\mathrm{weight}}{}}{(-){}_F^{\left(\mathrm{df}-D\right)}\left(\mathrm{df}-D,\left(\frac{1}{{}_F}+i\right)t_0\right)}.Math.{\left(1+{\left({}_F\right)}^2\right)}^{\frac{D-\mathrm{df}}{2}}\sin \left(\left(D-\mathrm{df}\right)a\tan \left({}_F\right)\right)& \left(9\right)\end{array}$

(80) Apparently, for df.fwdarw.D the Fourier sinus integral yields zero. Thus, the multiple reflections from the Euclidean part of the distribution do not contribute to . The different terms of Eq. 8 yield hence the following expression for the propagation of the wave:

(81) $\begin{array}{cc}\left(\right)=\underset{\underset{\mathrm{direct}\mathrm{beam}}{}}{}+\underset{\underset{\mathrm{reflected}\mathrm{beam}}{}}{{\left(1+{\left({}_F\right)}^2\right)}^{\frac{D-\mathrm{df}}{2}}\sin \left(\left(D-\mathrm{df}\right)a\tan \left({}_F\right)\right)},& \left(10\right)\end{array}$
with a scale factor for the direct beam (which is of the order of .sub.F) and

(82) $\begin{array}{cc}=\left({}_D*{}_F\right)/\left({}_D+{}_F\right)& \left(11\right)\\ {}_F=N{t}_0=\frac{155\left[\mathrm{.Math.}m\right]}{c_0}.Math.{}^{-\frac{1}{2}}& \left(12\right)\\ {}_D=M{t}_0>{}_F& \left(13\right)\\ \mathrm{df}=\frac{1}{15}\left(-5\right)+1& \left(14\right)\end{array}$

LIST OF REFERENCE NUMERALS

(83) 200 medical apparatus 202 computer 204 hardware interface 206 processor 208 user interface 210 computer storage 212 computer memory 214 shear wave data 216 mechanical property 218 power law relationship 220 concentration distribution of sonically dispersive elements 222 calibration data 224 control module 226 shear wave data processing module 228 power law determination module 230 concentration determination module 232 calibration data generation module 300 medical apparatus 302 medical imaging system 304 imaging zone 306 subject 308 subject support 310 mechanical actuator 312 mechanical actuator controller 314 shear waves 316 medical image data 318 image processing module 400 medical apparatus 402 magnetic resonance imaging system 404 magnet 406 bore of magnet 408 imaging zone 410 magnetic field gradient coils 412 magnetic field gradient coils power supply 414 radio-frequency coil 416 transceiver 420 pulse sequence 500 medical apparatus 506 high intensity focused ultrasound system 508 fluid filled chamber 510 ultrasound transducer 512 mechanism 514 mechanical actuator/power supply 516 path of ultrasound 518 ultrasound window 520 gel pad 522 focal point 600 excitation step 602 ultrasound transducer 604 matrix 606 particulate 608 shear wave imaging 610 reconstructed images 612 local value 613 plot of mechanical parameter vs. frequency 614 log of frequency 616 log of mechanical property 618 value 620 multiple values 700 power law fit 702 alpha 704 plot 706 spatial scale 708 filling factor 710 concentration inflection 712 abnormal fractal filling 714 normal Euclidean filling 800 plot 802 two solutions 1000 image 1002 particulates 1004 gel 1100 electromagnetic shaker 1102 carbon rod 1104 mechanical motion of rod 1106 cradle 1108 motion of cradle 1110 toothpick 1112 container 1114 gel 1116 cover 1118 shear waves 1200 frequency Hz 1202 Normalized complex shear modulus 1204 gel with 10 m microspheres 1206 power law fit to 1204 1208 reference gel 1210 power law fit to 1208 1304 fractal with gamma=0 1306 power law fit to 1304 1308 reference gel 1310 power law fit to 1308