# A broadband Ferromagnetic Resonance dipper probe for magnetic damping
measurements from 4.2 K to 300 K

Shikun He<EMAIL_ADDRESS>Division of Physics and Applied Physics, School of
Physical and Mathematical Sciences, Nanyang Technological University,
Singapore 637371  Christos Panagopoulos<EMAIL_ADDRESS>Division of
Physics and Applied Physics, School of Physical and Mathematical Sciences,
Nanyang Technological University, Singapore 637371

###### Abstract

A dipper probe for broadband Ferromagnetic Resonance (FMR) operating from 4.2
K to room temperature is described. The apparatus is based on a 2-port
transmitted microwave signal measurement with a grounded coplanar waveguide.
The waveguide generates a microwave field and records the sample response. A
3-stage dipper design is adopted for fast and stable temperature control. The
temperature variation due to FMR is in the milli-Kelvin range at liquid helium
temperature. We also designed a novel FMR probe head with a spring-loaded
sample holder. Improved signal-to-noise ratio and stability compared to a
common FMR head are achieved. Using a superconducting vector magnet we
demonstrate Gilbert damping measurements on two thin film samples using a
vector network analyzer with frequency up to 26 GHz: 1) A Permalloy film of 5
nm thickness and 2) a CoFeB film of 1.5 nm thickness. Experiments were
performed with the applied magnetic field parallel and perpendicular to the
film plane.

##  I introduction

In recent years, the switching of a nanomagnet by spin transfer torque (STT)
using a spin polarized current has been realized and intensively
studied.(Slonczewski, 1996; Katine  _et al._ , 2000; Mangin  _et al._ , 2006)
This provides avenues to new types of magnetic memory and devices, reviving
the interest on magnetization dynamics in ultrathin films.(Locatelli, Cros,
and Grollier, 2014; Brataas, Kent, and Ohno, 2012) High frequency techniques
play an important role in this research direction. Among them, Ferromagnetic
Resonance (FMR) is a powerful tool. Most FMR measurements have been performed
using commercially available systems, such as electron paramagnetic resonance
(EPR) or electron spin resonance (ESR).(Farle, 1998) These techniques take
advantage of the high Q-factor of a microwave cavity, where the field
modulation approach allows for the utilization of a lock-in
amplifier.(Heinrich, 2006) The high signal-to-noise ratio enables the
measurement of even sub-nanometer thick magnetic films. However, the operating
frequency of a metal cavity is defined by its geometry and thus is fixed. To
determine the damping of magnetization precession, which is in principle
anisotropic, several cavities are required to study the relation between the
linewidth and microwave frequency at a given magnetization direction.(Heinrich
_et al._ , 2003; Lenz  _et al._ , 2006; Heinrich  _et al._ , 2011)The apparent
disadvantage is that changing cavities can be tedious and prolong the
measurement time.

Recently, an alternative FMR spectrometer has attracted considerable
attention.(Kalarickal  _et al._ , 2006; Harward  _et al._ , 2011; Wei  _et
al._ , 2014; Montoya  _et al._ , 2014; Głowiński  _et al._ , 2014; Bailleul,
2013; Bilzer  _et al._ , 2007) The technique is based on a state of the art
vector network analyzer (VNA) and a coplanar waveguide (CPW). Both VNA and CPW
can operate in a wide frequency range hence this technique is also referred to
as broadband FMR or VNA-FMR. The broadband FMR technique offers several
advantages. First, it is rather straightforward to measure FMR over a wide
frequency range. Second, one may fix the applied magnetic field and acquire
spectra with sweeping frequency in a matter of minutes.(Bilzer  _et al._ ,
2007) Furthermore, a CPW fabricated on a chip using standard photolithography
enables FMR measurements on patterned films as well as on a single
device.(Neudecker  _et al._ , 2006) In brief, it is a versatile tool suitable
for the characterization of magnetic anisotropy, investigation of
magnetization dynamics and the study of high frequency response of materials
requiring a fixed field essential to avoid any phase changes caused by
sweeping the applied field.

Although homebuilt VNA-FMRs are designed mainly for room temperature
measurements, a setup with variable temperature capability is of great
interest both for fundamental studies and applications. Denysenkov et al.
designed a probe with variable sample temperature, namely, 4-420 K,
(Denysenkov and Grishin, 2003) however, the spectrometer only operates in
reflection mode. In a more recent effort, Harward et al. developed a system
operating at frequencies up to 70 GHz.(Harward  _et al._ , 2011) However the
lower bound temperature of the apparatus is limited to 27 K. Here we present a
2-port broadband FMR apparatus based on a superconducting magnet. A 3-stage
dipper probe has been developed which allows us to work in the temperature
range 4.2 - 300 K. Taking advantage of a superconducting vector magnet,
measurements can be performed with the magnetic moment saturated either
parallel or perpendicular to the film plane. We also designed a spring-loaded
sample holder for fast and reliable sample mounting, quick temperature
response and improved stability. This setup allows for swift changes of the
FMR probe heads and requires little effort for the measurement of devices. To
demonstrate the capability of this FMR apparatus we measured the temperature
dependence of magnetization dynamics of thin film samples of Permalloy (Py)
and CoFeB in different applied magnetic field configurations.

##  II apparatus

###  II.1 Cryostat and superconducting magnet system

Figure 1: View of the FMR dipper probe. Top panel: The schematic of the entire
design with a straight type FMR head. All RF connectors are 2.4 mm. The vacuum
cap mounted on the 4 K stage, using In seal, and the radiation shield mounted
on the second stage are not shown for clarity. Bottom panel: photograph of the
components inside the vacuum cap.

Our customized cryogenic system was developed by Janis Research Company Inc.
and includes a superconducting vector magnet manufactured by Cryomagnetics
Inc. A vertical field up to 9 T is generated by a superconducting solenoid.
The field homogeneity is $\pm$0.1% over a 10 mm region. A horizontal split
pair superconducting magnet provides a field up to 4 T with uniformity
$\pm$0.5% over a 10 mm region. The vector magnet is controlled by a Model
4G-Dual power supply. Although the power supply gives field readings according
to the initial calibration, to avoid the influence of remnant field we employ
an additional Hall sensor. The cryostat has a 50 mm vertical bore to
accommodate variable temperature inserts and dipper probes. Our dipper probe
described below is configured for this cryostat, however, the principle can be
applied also to other commercially available superconducting magnets and
cryostats.

###  II.2 Dipper probe

Fig. 1 shows a schematic of our dipper probe assembly and a photograph of the
components inside the vacuum cap. The dipper probe is 1.2 m long and is
mounted to the cryogenic system via a KF50 flange. The sliding seal allows a
slow insertion of the dipper probe directly into the liquid helium space.
Supporting arms (not shown) lock the probe and minimize vibration, with the
sample aligned to the field center. The connector box on top has vacuum tight
Lemo and Amphenol connectors for 18 DC signal feedthrough. Two 2.4 mm RF
connection ports allow for frequencies up to 50 GHz . A vacuum pump port can
be shut by a Swagelok valve. We adopted a three stage design as shown in the
photograph of Fig. 1. The 4 K stage and the vacuum cap immersed in the He bath
provide cooling power for the probe. The intermediate second stage acts as an
isolator of heat flow and as thermal sink for the RF cables, providing
improved temperature control. Furthermore, it allows one to change probe heads
conveniently as we discuss later. A separate temperature sensor on the second
stage is used for monitoring purpose. The third stage, namely, the FMR probe
head with the spring loaded sample holder, is attached to the lower end of the
intermediate stage using stainless steel rods.

A pair of 0.086” stainless steel Semi-Rigid RF cables run from the top of the
connector box to the non-magnetic bulkhead connector (KEYCOM Corp.) mounted on
the second stage. BeCu non-magnetic Semi-Rigid cables (GGB Industries, Inc.)
are used for the connection between the second stage and the probe head. The
cables are carefully bent to minimize losses. The rods connecting the stages
are locked by set screws. Loosening the set screw allows the rod length to be
adjusted to match the length of the RF cables. Reflection coefficient (S11)
and transmission coefficient (S21) can be recorded simultaneously with this
2-port design. The leads for the temperature sensors, heater, Hall sensor and
for optional transport measurements are wrapped around Cu heat-sinks at the 4
K stage before being soldered to the connection pins.

###  II.3 Probe head with spring-loaded sample holder

The key part of the dipper probe, namely, the FMR probe head is schematically
shown in Fig. 2. The assembly is placed in a radiation shield tube with an
inner diameter of 32 mm. To maximize thermal conduction between parts,
homebuilt components are machined from Au plated Cu. The 1” long customized
grounded coplanar waveguide (GCPW) has a nominal impedance of 50 Ohm. The
straight-line shape GCPW was made on duroid® R6010 (Rogers) board, with a
thickness of 254$\,\mu\text{m}$ and dielectric constant 10.2. The width of the
center conductor is 117 $\mu\text{m}$ and the gap between the latter and the
ground planes is 76 $\mu\text{m}$. For the connection, first the GCPW is
soldered to the probe head, and subsequently the center pin of the flange
connector (Southwest Microwave) is soldered to the center conductor of the
GCPW. The response of the dipper with the straight-line shape GCPW installed
is shown in Fig. 3. The relatively large insertion loss (-16.9 dB at 26 GHZ)
is due to a total cable length of more than 3 m and multiple connectors. The
high frequency current flowing in the CPW generates a magnetic field of the
same frequency. This RF field drives the precession of the magnetic material
placed on top of the signal line, and gives rise to a change in the system’s
impedance, which in turn alters the transmitted and reflected signals.

Figure 2: Schematic of the spring-loaded FMR probe head with straight shape
grounded coplanar waveguide (GCPW). 1 Au plated Cu housing; 2 straight shape
GCPW; 3 flange connector; 4 strain gauge thin film heater; 5 CernoxTM
temperature sensor; 6 Hall-sensor housing; 7 housing for 4-pin Dip socket or
pingo pin; 8 sample; 9 sample holder; 10 Cu spring; 11 spring housing; 12
sample holder handle nut.

A spring-loaded sample holder depicted in Fig. 2 by items 9 to 12 is designed
to mount the sample. The procedure for loading a sample is as follows: 1) Pull
up the handle nut and apply a thin layer of grease (Apiezon N type) to the
sample holder; 2) Place the sample at the center of the sample holder; 3)
Mount the spring-loaded sample holder to the FMR head; 4) Release the handle
nut gradually so that the spring pushes the sample towards the waveguide. The
mounting-hole of the spring-housing is slightly larger than the outer diameter
of the spring. This allows the sample holder to match the surface of the GCPW
self-adaptively. With the spring-loaded FMR head design, the sample mounting
is simple and leaves no residue from the commonly used tapes. It maximizes the
signal by minimizing the gap between waveguide and sample, and enhances the
stability. Furthermore, it is suitable for variable temperature measurements
due to the enhanced thermal coupling between the sample, cold head and sensors
( items 9 to 12 in Fig. 2.).

The temperature sensor is mounted at the backside of the probe head. Due to
limited space, the heater consists of three parallel connected strain gauges
with a resistance of 120 Ohm. The Hall sensor can be mounted according to the
required measurement configuration. The position of the Hall sensor shown in
Fig. 2 is an example for measurements in the presence of a magnetic field
applied parallel to the sample surface.

Figure 3: The reflection (S11) and transmission (S21) coefficients of the
dipper probe with the straight-line shape GCPW mounted. The measurement was
performed at room temperature.

###  II.4 Probe-head using end-launch connector

Although the probe head with straight-line CPW works well in our experiment,
the necessary replacement of CPW due to unavoidable performance fatigue over
time, or for testing new CPW designs can be time consuming. In response, end-
launch connectors (ELC) utilizing a clamping mechanism allow for a smooth
transition from RF cables to CPW. Soldering the launch pin to the center
conductor of CPW is optional and reduces the effort for modifications. In Fig.
4, we show our design of a FMR probe-head using ELC form Southwest Microwave,
Inc. and a homebuilt U-shape GCPW. Similar to the design of Fig. 2, a Au
plated Cu housing is used to mount the GCPW, ELC and the temperature and Hall
sensors. There are two locations for sample mounting. In position A, the
vertical field is used for measurements with the magnetic field applied
parallel to the surface of the thin film sample whereas the horizontal field
is used for measurements with field perpendicular to the sample surface. On
the other hand, measurements for both configurations can be accomplished only
by using the horizontal field if the sample is placed in position B. As shown
in Fig. 4 (b) and (c), to change between configurations simply requires
rotating the dipper probe by 90 degrees. Nevertheless, we prefer to place the
sample in position A for the parallel configuration since the solenoid field
is more uniform. However, we note that the same design with the sample placed
in position B is suitable also for an electromagnet. Furthermore, adding a
rotary stage to the probe enables angular dependent FMR measurements.

Figure 4: FMR probe-head with u-shape GCPW and end-launch connector. (a)
Photograph of the probe-head using end-launch connector and U-shape GCPW.
Sensors are mounted at the backside and at the bottom of the Cu housing.
Simplified sketch of the configuration for measuring with an external field
generated by the split coils (b) parallel and (c) perpendicular to the sample
plane. Rotating the dipper probe in the horizontal plane changes from one
configuration to the other.

##  III Experimental test

In this section, we present data to assess the performance of the FMR probe
head and discuss two sets of magnetic damping measurements, demonstrating the
capabilities and performance of the appratus.

###  III.1 Spring-loaded sample holder

We tested our setup using a Keysight PNA N5222A vector network analyzer with
maximum frequency 26.5 GHz. The output power of the VNA is always 0 dBm in our
test. Note that with 2.4 mm connectors and customized GCPW, our design can in
principle operate up to 50 GHz. The performance of the spring-loaded sample
holder is first studied at room temperature with a 2 nm thick
$\textrm{Co}{}_{40}\textrm{Fe}{}_{40}\textrm{B}_{20}$ film. For direct
comparison, the FMR spectra are recorded with two sample loading methods: One
with a spring-loaded sample holder (Fig. 2) and the other using the common
method(Harward  _et al._ , 2011) which only requires Kapton tape. The magnetic
field is applied parallel to the plane of the thin film sample. Six sets of
data were obtained by reloading the sample for each measurement. In Fig. 5, we
show the amplitude of the power transmission coefficient from Port 1 to Port 2
(S21) at a frequency of 10 GHz and a temperature of 300 K. The open circles
represent a spectrum for a spring-loaded sample mounting whereas the open
squares is the spectrum showing largest signal for the six flip-sample
loadings. The averaged spectra for all six spectra are shown by solid line and
dotted line, for spring and flip-sample loading, respectively. Two
observations are evident: First, the best signal we obtained using the flip
sample method is approximately 20 percent lower compared to the spring-loaded
method. Thus the spring-loaded method gives a better signal to noise ratio and
sensitivity. Second, for the spring-loaded method, the difference between the
averaged spectrum and single spectra is negligibly small. On the other hand
the variation between measurements for the flip-sample method can be as large
as 20 percent. Hence the spring-loaded method has better stability and is
reproducible.

Figure 5: Comparison between S21 signals obtained using spring-loaded sample
holder mounting and flip-sample mounting at 300 K. The sample has a stack of
MgO(3 nm)—CoFeB(2 nm)—MgO(3 nm) deposited on silicon substrate. (Numbers in
parenthesis of the sample composition represent the thickness of the
respective layer.) The frequency is 10 GHz and the FMR center field is at 1520
Oe.

###  III.2 Temperature response

As detailed in the previous section, the probe head is made of Au plated Cu
blocks with high internal thermal conduction and good thermal contact.
Consequently, the response time of the temperature control will be small as
the characteristic thermal relaxation time of a system is $C/k$, where $C$ is
the heat capacity and $k$ is the overall thermal conduction. Also, the
temperature difference between sample and sensor is minimized even with the
heater turned on. Shown in Fig. 6 are the FMR spectra and temperature
variation for a CoFeB film of 3 nm thickness measured at 4.4 K. The external
field was swept at a rate of about -10 Oe/s. For fields close to which FMR
peaks are observed, we detected a temperature rise of a few mK. In fact, the
field values corresponding to maximum temperatures are about 20 Oe lower than
the fields satisfying FMR condition, showing that the characteristic
relaxation time between the sample and cold head is approximately 2 seconds.
The temperature rise of the probe head due to FMR indicates that the magnetic
system absorbs energy from the microwave and dissipates into the thermal bath.
Specifically, at the field satisfying the FMR condition, the damping torque is
balanced by the torque generated by the RF field. However, the dissipation
power of such process is propotional to the thickness of the magnetic film
hence is very small. The successful detection of a temperature rise adds
credence to the high thermal conduction within the probe head and relative low
thermal conduction between different stages. This demonstration shows that the
probe head is capable of measuring samples with phase transitions in a narrow
temperature range, such as a superconducting/ferromagnetic bilayer
system.(Bell  _et al._ , 2008)

Figure 6: Sample temperature variation due to FMR at selected frequencies.
(upper panel) Amplitude of S21 and (lower panel) temperature variation of
MgO(3 nm)—CoFeB(3 nm)—MgO(3 nm) at 4.4 K measured with external field parallel
to the film plane.

###  III.3 Magnetic damping measurements

Although the FMR probe can be used to determine the energy anisotropy of
magnetic materials, our primary purpose is to study magnetic damping
parameter. In the following, two examples of such measurements are briefly
described. Shown in Fig. 7 is FMR response of a Py film of 5 nm thickness
deposited on a silicon substrate, measured at 4.4 K. The sweeping external
magnetic field is parallel to the sample surface. Real and Imaginary parts of
the spectra obtained at selected frequencies are plotted with open circles in
Fig. 7 (a) and (b), respectively. In FMR measurements, the change in the
transmittance, $S_{21}$, is a direct measure of the field-dependent
susceptibility of the magnetic layer. According to the Landau–Lifshitz–Gilbert
formalism, the dynamic susceptibility of the magnetic material in the
configuration where the field is applied parallel to the plane of the thin
film be described as:(Bilzer  _et al._ , 2006)

$\chi_{{\rm{\scriptscriptstyle IP}}}=\frac{4\pi M_{{\rm s}}\left(H_{0}+H_{{\rm
uni}}+4\pi M_{{\rm eff}}+i\,\Delta H/2\right)}{\left(H_{0}+H_{{\rm
uni}}\right)\left(H_{0}+H_{{\rm uni}}+4\pi M_{{\rm
eff}}\right)-H_{f}^{2}+i\,(\Delta H/2)\cdot\left[2\left(H_{0}+H_{{\rm
uni}}\right)+4\pi M_{{\rm eff}}\right]}$ (1)

Here, $4\pi M_{{\rm s}}$ is the saturation magnetization, $H_{{\rm uni}}$ is
the in-plane uniaxial anisotropy, $4\pi M_{{\rm eff}}$ is the effective
magnetization, $H_{f}=2\pi f/\gamma$, and $\Delta H$ is the linewidth of the
spectrum – the last term is of key importance to determine the damping
parameter. As shown by solid lines in Fig. 7 (a) and (b), the spectra can be
fitted very well by adding a background, a drift proportional to time, and a
phase factor.(Nembach  _et al._ , 2011; Kalarickal  _et al._ , 2006) The field
linewidth as a function of frequency – $\Delta H\,(f)$ is plotted in Fig. 7
(c). The data points fall on a straight line. The damping parameter
$\alpha_{{\rm{\scriptscriptstyle GL}}}=0.012\pm 0.001$ is therefore determined
by the slope through(Celinski and Heinrich, 1991; Lenz  _et al._ , 2006):

$\Delta H=\frac{4\pi}{\gamma}\alpha_{{\rm{\scriptscriptstyle GL}}}\,f+\Delta
H_{0}$ (2)

The error bar here is calculated from the confidence interval of the fit.

Figure 7: FMR data of a Py thin film of thickness 5 nm measured at 4.4 K with
magnetic field applied parallel to the sample plane. (a) Real and (b)
Imaginary parts of transmitted signal S21 at selected frequencies. The data
are normalized and the relative strength between the spectra at different
frequencies are kept. (c) FMR linewidth as a function of frequency. The
damping was calculated to be 0.012 $\pm 0.001$, using a linear fit.

We have also tested the setup with a magnetic field applied perpendicular to
the sample plane. The results for a MgO (3 nm)—Co40Fe40B20(1.5 nm)—MgO (3nm)
stack deposited on silicon substrate are shown in Fig. 8. Comparing the
spectra obtained at different temperatures and fixed frequency, two
observations are evident. First, the FMR peak position shifts to higher field
as the temperature is lowered due to changes in the effective magnetization.
Second, the FMR linewidth increases with decreasing temperature. Although the
interfacial anisotropy can be determined by fitting the FMR peak positions to
the Kittel formula,(Sparks, Loudon, and Kittel, 1961) here, we are more
interested in the damping parameter as a function of temperature. The dynamic
susceptibility in this configuration is(Devolder  _et al._ , 2013):

$\chi_{{\scriptscriptstyle{\rm OP}}}=\frac{4\pi M_{{\rm s}}\left(H-4\pi
M_{{\rm eff}}-i\,\Delta H/2\right)}{\left(H-4\pi M_{{\rm
eff}}\right)^{2}-H_{f}^{2}+i\,\Delta H\cdot\left(H-4\pi M_{{\rm eff}}\right)}$
(3)

Following the same procedure as for Py, the real and imaginary part of the
spectra are fitted simultaneously to obtain the linewidth. In Fig. 8 (b), we
plot the linewidth as a function of frequency at the two boundaries of our
measured temperatures. Although the measured linewidth at lower temperature is
larger, the slope of the two curves is in good agreement. The additional
linewidth at 6 K is primarily due to zero frequency broadening, which
quantifies the magnitude of dispersion of the effective magnetization. The
results are summarized in Fig. 8 (c). Gilbert damping is essentially
independent of temperature although there is a minimum at 40 K. The room
temperature value obtained is in agreement with the value for a thicker
CoFeB.(Bilzer  _et al._ , 2006; Ikeda  _et al._ , 2010) On the other hand, the
inhomogeneous broadening increases with lowering temperatue. The value at 6 K
is more than double compared to room temperatue. Notably, neglecting the zero-
frequency offset $\Delta$H(0), arising due to inhomogeneity, would give rise
to an enhanced effective damping compared to the intrinsic contribution.
Cavity based, angular dependent FMR may also distinguish the Gilbert damping
from inhomogeneity effects. A shortcoming however, is the need to take into
account the possible contribution of two magnon scattering, which causes
increased complications in the analysis of the data.(Mizukami, Ando, and
Miyazaki, 2002; Lindner  _et al._ , 2009) On the other hand, broadband FMR
using a dipper probe with the applied magnetic field in the perpendicular
configuration, rules out two magnon scattering making this technique
relatively straightforward to implement.(Arias and Mills, 1999)

The dipper probe discussed here is not limited to measurements of the damping
coefficient. The broadband design is also useful for time-domain
measurements.(Cui  _et al._ , 2010) Furthermore, a spin transfer torque
ferromagnetic resonance(Liu  _et al._ , 2011) measurement on a single device
can be performed with variable temperature using bias Tee and a separate
sample holder.

Figure 8: Temperature dependent FMR measurement for a CoFeB thin film of
thickness 1.5 nm with the magnetic field applied perpendicular to the film
plane. (a) Transmitted FMR signal at 20 GHz obtained at different
temperatures. (b) FMR linewidth as a function of frequency at 6 K and 280 K.
(c) Damping constant and inhomogeneous broadening as a function of
temperature. The solid lines are the guides for the eye.

##  IV conclusion

We have developed a variable temperature FMR to measure the magnetic damping
parameter in ultra thin films. The 3-stage dipper and FMR head with a spring-
loaded sample holder design have a temperature stability of milli Kelvin
during the FMR measurements. This apparatus demonstrates improved signal
stability compared to traditional flip-sample mounting. The results for Py and
CoFeB thin films show that the FMR dipper can measure the damping parameter of
ultra thin films with: Field parallel and perpendicular to the film plane in
the temperature range 4.2-300 K and frequency up to at least 26 GHz.

###### Acknowledgements.

The authors are grateful to Sze Ter Lim at Data Storage Institute for
preparing the CoFeB samples. We acknowledge Singapore Ministry of Education
(MOE), Academic Research Fund Tier 2 (Reference No: MOE2014-T2-1-050) and
National Research Foundation (NRF) of Singapore, NRF-Investigatorship
(Reference No: NRF-NRFI2015-04) for the funding of this research.

## References

  * Slonczewski (1996) J. C. Slonczewski, Journal of Magnetism and Magnetic Materials 159, L1 (1996).
  * Katine  _et al._ (2000) J. A. Katine, F. J. Albert, R. A. Buhrman, E. B. Myers, and D. C. Ralph, Physical Review Letters 84, 3149 (2000).
  * Mangin  _et al._ (2006) S. Mangin, D. Ravelosona, J. A. Katine, M. J. Carey, B. D. Terris, and E. E. Fullerton, Nat Mater 5, 210 (2006).
  * Locatelli, Cros, and Grollier (2014) N. Locatelli, V. Cros, and J. Grollier, Nat Mater 13, 11 (2014).
  * Brataas, Kent, and Ohno (2012) A. Brataas, A. D. Kent, and H. Ohno, Nat Mater 11, 372 (2012).
  * Farle (1998) M. Farle, Reports on Progress in Physics 61, 755 (1998).
  * Heinrich (2006) B. Heinrich,  _Measurement Techniques and Novel Magnetic Properties_ (Springer Berlin Heidelberg, 2006).
  * Heinrich  _et al._ (2003) B. Heinrich, G. Woltersdorf, R. Urban, and E. Simanek, Journal of Applied Physics 93, 7545 (2003).
  * Lenz  _et al._ (2006) K. Lenz, H. Wende, W. Kuch, K. Baberschke, K. Nagy, and A. Jánossy, Physical Review B 73, 144424 (2006).
  * Heinrich  _et al._ (2011) B. Heinrich, C. Burrowes, E. Montoya, B. Kardasz, E. Girt, Y.-Y. Song, Y. Sun, and M. Wu, Physical Review Letters 107, 066604 (2011).
  * Kalarickal  _et al._ (2006) S. S. Kalarickal, P. Krivosik, M. Wu, C. E. Patton, M. L. Schneider, P. Kabos, T. J. Silva, and J. P. Nibarger, Journal of Applied Physics 99, 093909 (2006).
  * Harward  _et al._ (2011) I. Harward, T. O’Keevan, A. Hutchison, V. Zagorodnii, and Z. Celinski, Review of Scientific Instruments 82, 095115 (2011).
  * Wei  _et al._ (2014) J. Wei, J. Wang, Q. Liu, X. Li, D. Cao, and X. Sun, Review of Scientific Instruments 85, 054705 (2014).
  * Montoya  _et al._ (2014) E. Montoya, T. McKinnon, A. Zamani, E. Girt, and B. Heinrich, Journal of Magnetism and Magnetic Materials 356, 12 (2014).
  * Głowiński  _et al._ (2014) H. Głowiński, M. Schmidt, I. Gościańska, J.-P. Ansermet, and J. Dubowik, Journal of Applied Physics 116, 053901 (2014).
  * Bailleul (2013) M. Bailleul, Applied Physics Letters 103, 192405 (2013).
  * Bilzer  _et al._ (2007) C. Bilzer, T. Devolder, P. Crozat, C. Chappert, S. Cardoso, and P. P. Freitas, Journal of Applied Physics 101, 074505 (2007).
  * Neudecker  _et al._ (2006) I. Neudecker, K. Perzlmaier, F. Hoffmann, G. Woltersdorf, M. Buess, D. Weiss, and C. H. Back, Physical Review B 73, 134426 (2006).
  * Denysenkov and Grishin (2003) V. P. Denysenkov and A. M. Grishin, Review of Scientific Instruments 74, 3400 (2003).
  * Bell  _et al._ (2008) C. Bell, S. Milikisyants, M. Huber, and J. Aarts, Physical Review Letters 100, 047002 (2008).
  * Bilzer  _et al._ (2006) C. Bilzer, T. Devolder, J.-V. Kim, G. Counil, C. Chappert, S. Cardoso, and P. P. Freitas, Journal of Applied Physics 100, 053903 (2006).
  * Nembach  _et al._ (2011) H. T. Nembach, T. J. Silva, J. M. Shaw, M. L. Schneider, M. J. Carey, S. Maat, and J. R. Childress, Physical Review B 84, 054424 (2011).
  * Celinski and Heinrich (1991) Z. Celinski and B. Heinrich, Journal of Applied Physics 70, 5935 (1991).
  * Sparks, Loudon, and Kittel (1961) M. Sparks, R. Loudon, and C. Kittel, Physical Review 122, 791 (1961).
  * Devolder  _et al._ (2013) T. Devolder, P.-H. Ducrot, J.-P. Adam, I. Barisic, N. Vernier, J.-V. Kim, B. Ockert, and D. Ravelosona, Applied Physics Letters 102, 022407 (2013).
  * Ikeda  _et al._ (2010) S. Ikeda, K. Miura, H. Yamamoto, K. Mizunuma, H. D. Gan, M. Endo, S. Kanai, J. Hayakawa, F. Matsukura, and H. Ohno, Nat Mater 9, 721 (2010).
  * Mizukami, Ando, and Miyazaki (2002) S. Mizukami, Y. Ando, and T. Miyazaki, Physical Review B 66, 104413 (2002).
  * Lindner  _et al._ (2009) J. Lindner, I. Barsukov, C. Raeder, C. Hassel, O. Posth, R. Meckenstock, P. Landeros, and D. L. Mills, Physical Review B 80, 224421 (2009).
  * Arias and Mills (1999) R. Arias and D. L. Mills, Physical Review B 60, 7395 (1999).
  * Cui  _et al._ (2010) Y.-T. Cui, G. Finocchio, C. Wang, J. A. Katine, R. A. Buhrman, and D. C. Ralph, Phys. Rev. Lett. 104, 097201 (2010).
  * Liu  _et al._ (2011) L. Liu, T. Moriyama, D. C. Ralph, and R. A. Buhrman, Phys. Rev. Lett. 106, 036601 (2011).


# Rapid cargo verification with cosmic ray muon scattering and absorption
tomography

Anzori Sh. Georgadze,11footnotetext: Corresponding author.

###### Abstract

Cosmic ray muon tomography is considered a promising method for the non-
invasive inspection of shipping containers and trucks. It utilizes highly
penetrating cosmic-ray muons and their interactions with various materials to
generate three-dimensional images of large and dense, like inter-modal
shipping containers, typically not transparent with conventional X-ray
radiography technique. The commonly used methods for imaging with muons are
based on muon scattering or absorption-transmission data analysis. Due to
large thickness of cargo material in shipping container substantial scattering
and absorption occur when muons passing through cargo.

One of the key tasks of customs and border security is to verify shipping
container declarations to prevent illegal trafficking, and muon tomography
could be a viable choice for this task. In this paper, we demonstrate through
Monte Carlo simulations using the GEANT4 toolkit that a combined analysis of
muon scattering and absorption data can improve the identification of cargo
materials compared to using scattering or absorption data alone. The
statistical differences in scattering and absorption data for several cargo
materials are quantified. For a particular smuggling scenario where tobacco
declared as paper towel rolls, it is demonstrated that the combined analysis
can accurately distinguish between tobacco and paper towel rolls with
5.5$\sigma$ accuracy for detector spatial resolution (FWHM) of 0.235 mm,
4.5$\sigma$ for 1.175 mm resolution (FWHM), and 3.5$\sigma$ accuracy for2.35
mm spatial resolution (FWHM), in a short scanning time of 10 seconds. This
rapid detection capability has significant implications for anti-smuggling
efforts and cargo inspection.

##  1 Introduction

Muon tomography is an emerging technique with promising applications in
various fields such as non-destructive testing, border security, and
archaeology [1, 2, 3, 4, 5, 6]. The principle behind muon scattering
tomography is similar to medical X-ray imaging, but instead of using X-rays,
it utilizes muons. This technology has several advantages, including its
ability to penetrate dense materials without causing harm and its capability
to provide detailed images of large and massive objects. Traditional methods
of cargo inspection, such as X-ray scanning, require high-intensity sources of
hazardous radiation, which have limited penetration power in dense materials
and may not provide a comprehensive view of container or cargo contents.

Initially, the goal of muon tomography was to detect nuclear material
smuggling to prevent the threat of nuclear terrorism. Later, muon tomography
was proposed for application in border security to scan cargo containers and
vehicles for hidden contraband or illicit materials without the need for
physical inspection, thereby enhancing security while minimizing disruptions
to trade and travel [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17].

Muon tomography utilizes cosmic ray muons, which are the secondary particles
that come from the extensive atmosphere shower of high-energy cosmic rays from
space, mainly protons with a flux about $\approx$ 10000/m2/minute. The
distribution of zenith angle $\theta$ is known to follow a cosine-squared law,
such that $\textit{I}(\theta)$ = $\textit{I}_{0}\cos^{2}\theta$. Muons are a
natural source of radiation that rain down upon the Earth with an average
energy of 3–4 GeV. Their mass is approximately 207 times heavier than the mass
of an electron [18]. As these muons pass through matter they undergo multiple
Coulomb scattering. The degree of scatter observed is dependent on the Z of
material.

When muons interact with matter, they scatter, and by measuring the angles and
intensities of these scattered muons it is possible to create 3D images of the
interior of objects or structures. In addition to scattering, cosmic muons
also absorbed by the materials they pass through. Muon absorption depends on
the density and composition of the material. When a muon interacts with
matter, it can lose energy through various processes such as ionization,
bremsstrahlung, and nuclear interactions. There are two main types of cosmic
ray muon imaging techniques: absorption-based and scattering-based.
Absorption-based imaging exploits the principle that the flux of muons
reaching a particle detector decreases with the density of the material they
traverse. Dense materials, such as metals or high-density objects, attenuate
or absorb more muons than less dense materials. On the other hand, scattering-
based tomography, relies on the scattering of muons by the nuclei of the
materials they pass through. By analyzing the intensity and angular
distribution of muons before and after passing through the sea container or
truck internal structure and density of cargo can be studied.

In our study we have performed detailed GEANT4 simulations of a muon
tomography system that employs detectors with high spatial resolution based on
scintillator technology. To improve the accuracy of verifying the contents of
shipping containers, we developed a novel approach that combines muon
scattering and absorption data. This analysis enabled significantly better
material discrimination compared to analyzing only muon scattering data,
thereby decreasing measurement time, which is critical for maintaining high
throughput in container verification.

##  2 Muon tomography

The scattering muon tomography technique calculates the deflection of muons
from their straight trajectory due to multiple coulomb interactions, which in
turn depend on cargo density and chemical composition. The absorption muon
tomography technique calculates the fraction of muons that were stopped in
cargo. The absorption also depend on cargo density and chemical composition.
Due to the different sizes and loading configurations of the cargo, object
detection techniques can be applied to localize and identify the dimensions of
the cargo in the reconstructed image of the container.

Figure 1: (a) MTS layout built in GEANT4 simulates a tobacco smuggling
scenario. Tobacco fills the entire volume of the container (yellow boxes), and
a row of paper towel rolls (blue boxes) is placed between the tobacco and the
doors to veil the tobacco; (b) Visualization of muon scattering tracks (blue)
and muon stopped (red) tracks in cargo.

Figure 1(b) shows examples of scattering muon tracks and stopped muon tracks.
On the left two scattered muons that went through the cargo are shown in blue
color. The right particles (in red color) are stopped in cargo before reaching
the lower tracker detectors.

###  2.1 Muon scattering tomography

Cosmic ray muon undergoes multiple coulomb scattering while passing through
material. The angular distribution of scattered muon of momentum p is
approximately Gaussian, with zero mean and standard deviation given by:

$\sigma_{\theta}=\frac{13.6MeV}{\beta
cp}\sqrt{\frac{L}{{}_{X_{0}}}}(1+0.038)ln\frac{L}{{}_{X_{0}}})$ (2.1)

where $\beta$ is the ratio between velocity of muon _V_ to velocity of light
_c_ , _X 0_ is the radiation length of the material, _L_ is the length of the
material traversed. _X 0_ is a material property and depends on the density of
the material $\rho$, the atomic mass _A_ and the atomic number _Z_ and can be
expressed as[19]:

$X_{0}=\frac{716.4g/cm^{2}}{\rho}\frac{A}{Z(Z+1)ln(\frac{287}{\sqrt{Z}})}$
(2.2)

For material discrimination purposes in muon tomography we define the
scattering density ($\lambda$), derived from Rossi formula [20]:

$\lambda(X_{0})=\left(\frac{13.6\;\mbox{MeV}}{p_{0}}\right)^{2}\frac{1}{X_{0}}$
(2.3)

where $p_{0}$ is a nominal muon momentum and $\lambda$ is with units in
milliradians2/cm. The approximation $\beta$c $\cong$ 1 for muons is used in
this equation. The scattering density of a material represents the mean square
scattering angle of muons passing through a unit depth of that material. For
generating tomographic image scattering angles are accumulated at different
locations generating a 3D map of scattering densities.

From equations 2 and 3 one can seen that the width of scattering angle
distribution inversely depends on the material radiation length which is
dependent on the atomic number of the material and material density. For
materials with higher atomic numbers Z (high-Z) muons will scatter with larger
scattering angles. It is therefore possible to characterize the density
profile of a large object by measuring muon directions before and after
passing through the target object. If cargo is a composite material muon
scattering depend on its chemical composition. The combined radiation length
of the composite material is given by formula [21]:

$\frac{W_{0}}{X_{0}}=\sum\frac{W_{i}}{X_{i}}$ (2.4)

where W0 is the total mass of the sample in g. X0 is the combined radiation
length of the sample in g$\cdot$cm-2. Wi is the mass of the individual
component in g. Xi is the radiation length of the individual component
g$\cdot$cm-2.

###  2.2 Absorption muon tomography

While passing cargo muons loosing their energy due to ionization and have a
probability to be stopped and decay. The attenuation length depends on the
density of the matter traversed. A muon before being stopped, will be able to
cross on average a certain opacity X defined as X = $\int\cdot\rho\cdot$dL,
where $\rho$ is the density of the material and L is the particle path length.
In the absorption tomography approach, tracks of muon, that were stopped in
cargo can be reconstructed using algorithm described in publication [22]. This
algorithm analyzes muon tracks, that were passed upper detectors but were not
detected in lower detectors but lies on Line-of-Response (LOR), which
connecting upper and lower detectors. Applying a muon absorption algorithm
allows for the identification and counting of muons that have been stopped
within the cargo. Since muon absorption depends on the cargo’s effective
atomic number and density, this algorithm allow identify cargo material.

Following approach, developed in publication [23] we introduce a parameter,
describing muon absorption as the “stopping power” ($SP$) of a cargo material,
which is given by the expression

$SP=\frac{(N_{abs}/V/t)\times\left\langle p\right\rangle}{(N_{stat}/V/t)\times
d_{th}}$ (2.5)

where $\left\langle p\right\rangle$ is the average momentum of the incident
cosmic rays, Nscat the total number of scattering tracks, Nabs total number of
tracks, which had been absorbed (or deflected) in the shipping container
tracks, $d_{th}$ is a thickness of cargo material, $V$ – cargo volume, $t$ –
scanning time. For each material the total number of scattered tracks Nscat
and the number of absorbed tracks Nabs in the shipping container and cargo was
calculated. The ratio of stopping power to scattering, while the latter
expressed in following way

$\lambda=\frac{(\left\langle\theta\right\rangle\left\langle
p\right\rangle)^{2}}{d_{th}}$ (2.6)

where $\left\langle\theta\right\rangle$ is an average scattering angle of
material.

##  3 Monte Carlo simulations

The geometry of a muon tomography station (MST) is shown in figure 1a. MST is
consisting of two muon tracking modules above and below of the shipping
container. Each tracking module includes two position-sensitive detectors,
modeled as a plane made of plastic scintillator with a detection efficiency of
100% and dimension of 8 m $\times$ 4 m $\times$ 1 mm, covering a shipping
container. The distance between two position-sensitive detector planes is 10
cm and distance between upper and lower modules is 3 m. The simulation of CRT
and cargo in sea container is based on GEANT4 toolkit [4]. For generation of
cosmic ray muons at sea level we use the Cosmic-Ray Shower Library (CRY) [24].
The origin points of generated muons were sampled from a horizontal plane
surface of 10 m $\times$ 10 m. Generated muons are interfaced with GEANT4 to
simulate the interaction of muons with the detector, shipping container and
cargo. The simulated data samples were produced by generating 5000 datasets
for each material. Datasets were simulated by tracking 100000 muons, that
corresponds to scanning time of $\approx$10 seconds. We use the standard
physics list for high-energy particle transport named "FTFP_BERT". Synthetic
datasets were analyzed using ROOT data analysis package [25]. The prototype of
tracking detector based on multilayer array of plastic scintillating fibers
readout with Silicon Photomultipliers (SiPMs) was constructed and tested [26].
The obtained spatial resolution in experimental runs was 120 microns (WFHM).
To mimic the material discrimination accuracy of detectors with different
spatial resolutions, the hit positions of muons produced in the detector
plates were smeared using a Gaussian function with FWHM values of 0.235 mm,
1.17 mm, and 2.35 mm.

Contraband materials can be disguised as legal materials in customs
declarations to evade detection. Muon Scattering Tomography (MST) can verify
the real density and atomic number of the material in the container,
discriminate between different products, and allow customs officials to
compare the actual contents of a shipping container with what is declared in
the customs documentation. According to research [27], the mean cargo density
is just under 0.2 g/cm3. As a source of data for the selection of simulation
scenarios one can use a PIERS (Port Import/Export Reporting Service) United
States import data set. The sheer volume and complexity of the flow of
commercial goods make it impractical to model every possible smuggling
scenario. In this work, we have simulated scenarios of uniform cargo loading,
which are most commonly used in international trade traffic. For the
simulations we considered cargo materials in bulk form. The bulk density of a
material is a combination of the density of solid particles, which includes
occluded air, the density of the solids themselves, and the interstitial air.
Data on the density of bulk materials we used in simulations were taken from
Web Tech AutoWeigh Material Bulk Density Chart [28, 29].

Table 1: Bulk material properties [28] of the materials used in the simulations. Bulk material | Bulk density (g/cm3) | Bulk material | Bulk density (g/cm3)
---|---|---|---
Paper (towel rolls) | 0.120 | Sugar | 0.800
Tobacco (cigarette) | 0.190 | Dry pasta | 0.950
Coffee (Roasted (Beans)) | 0.368 | Graphite (granules) | 0.720
Tea | 0.433 | Marble (granular) | 1.089
Polyethylene (Pellets) | 0.561 | Cement | 1.362
Rice | 0.71 | Copper (Fines) | 1.618

##  4 Image reconstruction

Scattering muon tomography is based on the measurement of muon deflections
when passing through the cargo. We calculate deflection using reconstruction
algorithm based on the Point of Closest Approach (PoCA) method [30] which
calculates the closest point of approach between two 3D lines. The PoCA
algorithm makes the simplified assumption that the muon scattering occurs in a
single-point and searches for the point of closest approach between the
incident $\textit{v}_{in}$ and outgoing $\textit{v}_{out}$ reconstructed muon
track directions. The shortest line segment between the incident and outgoing
tracks is estimated by finding the pair of points of closest approach between
the tracks. Midpoint of this line segment is considered as a scattering center
of the muon. The scattering angle between the two tracks is calculated with
the formula [19]:

$\theta_{scatt}=\arccos(\frac{\vec{v}_{in}\times\vec{v}_{out}}{\left|\upsilon_{1}\right|\left|{\upsilon_{2}}\right|})$
(4.1)

Degree of muon deflection while passing the cargo in shipping container depend
on material density and its chemical composition, thus an amount of
reconstructed scattering centers using the PoCA approximation can be used for
the material identified.

Figure 2: (a) Reconstructed tomographic image of cargo with dimensions 3 m
$\times$1.15 m $\times$1.15 m, (b) reconstructed X–Y image; (c) projection
onto the X-axis of the tomographic image.

##  5 Results

The simulation of muon passage through MTS and cargo enables the production of
data sets for statistical analysis. Hits on detector planes are recorded for
each event, enabling the reconstruction of muon tracks using a ROOT analysis
code. Four hits detected in upper and lower trackers, as shown on figure 1(a)
allow calculate scattering angle and scattering center position using PoCA
algorithm. If two hits are detected in upper trackers and LoR is passing the
cargo than trajectory of stopped muon is calculated.

For the calculation of scattering density and stopping power, cargo dimensions
are determined from the reconstructed tomographic image. In the 1D projections
of the reconstructed tomographic image obtained from a 10-second scan, the
edges of the cargo object can be identified using a ROOT script that
calculates the first derivative, as described in publication [17]. In figure
2, tobacco packed into cardboard boxes forms an approximately cubic shape with
dimensions 3 m $\times$1.15 m $\times$1.15 m. The detected cargo edges are
found at positions -147.5 cm and 142.5 cm, while the actual positions are -150
cm and 150 cm. Simulated cargo materials were represented by a
5.85$\times$2.30$\times$2.30 m3 rectangular volume, almost fully filling the
shipping container. In the simulation run, the calculated scattering density
and scattering-to-stopping ratio are recorded. For each cargo material, 5,000
datasets were simulated. As shown in figure 3(c), the two-dimensional
histogram of scattering and scattering-to-stopping ratio density distributions
for different materials are well separated in 2D space and localized in
distinct regions, demonstrating the material discrimination ability. The data
show different behaviors in three regions. For organic materials up to tea,
there is minimal scattering-to-stopping ratio variation but a stronger
scattering density variation. Other organic materials—such as polyethylene
(pellets), rice, sugar, and dry pasta—show a larger variation in scattering-
to-stopping ratio than in scattering rate. Inorganic materials demonstrate a
greater variation in scattering-to-stopping ratio, with scattering density
changing minimally. The impact of chemical composition is visibly distinct.

Figure 3: (a) Histogram of the scattering-to-stopping ratio distributions for
simulated materials for 10 seconds scanning time; (b) histogram of the
scattering density distributions for simulated materials; (c) scatter plot of
the scattering-to-stopping ratio distributions versus scattering density
distributions for simulated materials.

To numerically quantify the accuracy of material discrimination, we considered
a scenario involving the smuggling of tobacco (cigarettes) declared as paper
towel rolls. We simulated datasets for container fully loaded with 5 million
cigarettes and container fully loaded with paper towel rolls. In figure 4 the
scatter plot of scattering density versus scattering-to-stopping ratio are
shown for paper towel rolls and tobacco for detector spatial resolutions of
0.235 mm, 1.17 mm, and 2.35 mm (FWHM). Colors indicate group membership of
distributions. Colors indicate group membership of distributions. On the top
of each 2D histogram in figure 4 we show the 1D histograms of scattering data,
while on the right side we show the histogram of scattering-to-stopping ratio
data. Compared to the scatter versus scattering-to-stopping ratio 2D map, the
1D distributions are limited due to the partial overlap of distributions for
paper towel rolls and tobacco, making materials discrimination less
distinctive. This demonstrates the superiority of 2D data analysis over 1D
data analysis. These plots demonstrate the impact of detector position
resolution on material discrimination analysis. As can be seen, there is
expected degradation of discrimination of accuracy at worse spatial
resolution.

For quantitative analysis of discrimination accuracy, we applied the Gaussian
Mixture Model (GMM) algorithm, which helps categorize data into groups based
on their probability distributions.

Figure 4:  Scatter plot showing the distributions of scattering and
scattering-to-stopping ratio for paper towel rolls and tobacco. Data points
for paper towel rolls are marked in blue, and data points for tobacco are
marked in red. On the top of each 2D histogram we show the 1D histograms of
the scattering distributions for paper towel rolls and tobacco, while on the
right side we show the histogram of the corresponding scattering-to-stopping
ratio distributions. Panels (a), (b), and (c) of this figure show scatter
distributions and confidence ellipses for detector resolutions of 0.235 mm,
1.17 mm, and 2.35 mm (FWHM), respectively.

Figure 5: (a) The scatter distributions of the data samples were fitted with
2D Gaussian functions. The confidence ellipses for each distribution are set
to show 1, 2, 3 and 4 $\sigma$ CL’s. Panels (a), (b), and (c) of this figure
show scatter distributions and confidence ellipses for detector resolutions of
0.235 mm, 1.17 mm, and 2.35 mm (FWHM), respectively.

We initialize a 2D Gaussian GMM with two components and fit it to the data on
the scatter plot (see figure 5). The confidence ellipses for each distribution
are set to show 1, 2, 3 $\sigma$ CL and confidence ellipse at which
distributions are discriminated. As can be seen, the scattering versus
scattering-to-stopping ratio distributions for towel paper and tobacco are
accurately discriminated with more than 4 $\sigma$ confidence level (CL) for
the detector spatial resolution of 0.235 mm, 1.17 mm. For the spatial
resolution of 2.35 mm (FWHM) discrimination accuracy is 3.5 $\sigma$ CL.

##  6 Conclusions

In this work, we used Monte Carlo simulations based on the GEANT4 package to
investigate the ability of cosmic ray muon tomography to perform prompt
characterization of shipping container contents and verification of customs
declarations. We developed a method of combined scattering and scattering-to-
stopping ratio data analysis to improve the material discrimination ability of
cosmic ray muon tomography. This method was applied to analyze simulated
datasets for several cargo materials. We demonstrate how integrating
scattering and scattering-to-stopping ratio data into a two-dimensional
representation enables clearer differentiation between material
classifications. The developed approach leverages the combined information
from both types of data to enhance the accuracy of material discrimination in
muon tomography applications, providing a powerful tool for detecting smuggled
goods disguised as other items. For the specific case of tobacco smuggling,
the combined analysis of scatter-absorption data allows for the accurate
discrimination between tobacco and paper towel rolls with 5.5 $\sigma$
accuracy for detector spatial resolution (FWHM) of 0.235 mm, 4.5 $\sigma$ for
1.175 mm resolution (FWHM), and 3.5 $\sigma$ which corresponds to 99.95%
discrimination accuracy for 2.35 mm spatial resolution (FWHM) within a
10-second scanning time.

The ability to rapidly characterize cargo materials comes from the analysis of
muon scattering and absorption data for large amounts of cargo material
weighing several tons, resulting in high statistical precision. Thanks to the
muon penetration power and the developed method of scatter-absorption combined
analysis, cosmic-ray muon tomography can be an effective solution for non-
invasive cargo inspections, aiding in the fight against contraband
trafficking.

##  7 Acknowledgments

This work was partially funded by the EU Horizon 2020 Research and Innovation
Programme under grant agreement no. 101021812 (“SilentBorder”).

## References

  * [1] L. Bonechi, R. D’Alessandro and A. Giammanco, _Atmospheric muons as an imaging tool_ , _Reviews in Physics_ 5 (2020) 100038.
  * [2] S. Barnes et al., _Cosmic-ray tomography for border security_ , _Instruments_ 7 (2023) 13.
  * [3] K.N. Borozdin et al., _Radiographic imaging with cosmic-ray muons_ , _Nature_ 422 (2003) 277.
  * [4] L. Cuéllar et al., _Soft cosmic ray tomography for detection of explosives_ , in _2009 IEEE Nuclear Science Symposium Conference Record (NSS/MIC)_ , pp. 968–970, 2009, DOI.
  * [5] Z. Yifan, Z. Zhi, Z. Ming, W. Xuewu and Z. Ziran, _Discrimination of drugs and explosives in cargo inspections by applying machine learning in muon tomography_ , _High Power Laser and Particle Beams_ 30 (2018) 086002.
  * [6] J. Chen et al., _Towards a muon scattering tomography system for both low-Z and high-Z materials_ , _Journal of Instrumentation_ 18 (2023) P08008.
  * [7] L.J. Schultz, _Cosmic ray muon radiography_ , Portland State University (2003).
  * [8] E. Åström, G. Bonomi, I. Calliari, P. Calvini, P. Checchia, A. Donzella et al., _Precision measurements of linear scattering density using muon tomography_ , _Journal of Instrumentation_ 11 (2016) P07010.
  * [9] P. Checchia, _Review of possible applications of cosmic muon tomography_ , _Journal of Instrumentation_ 11 (2016) C12072.
  * [10] V. Antonuccio, M. Bandieramonte, U. Becciani, D. Bonanno, G. Bonanno, D. Bongiovanni et al., _The muon portal project: Design and construction of a scanning portal based on muon tomography_ , _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 845 (2017) 322.
  * [11] C. Pugliatti, V. Antonuccio, M. Bandieramonte, U. Becciani, F. Belluomo, M. Belluso et al., _Design of a muonic tomographic detector to scan travelling containers_ , _Journal of Instrumentation_ 9 (2014) C05029.
  * [12] C. Morris, J. Bacon, K. Borozdin, H. Miyadera, J. Perry, E. Rose et al., _A new method for imaging nuclear threats using cosmic ray muons_ , _AIP Advances_ 3 (2013) .
  * [13] E. Preziosi, F. Arcieri, A. Caltabiano, P. Camarri, S. Casagrande, D. Cavicchioni et al., _Tecnomuse: a novel, rpc-based, muon tomography scanner for the control of container terminals_ , _Journal of Physics: Conference Series_ 1548 (2020) 012021.
  * [14] A.S. Georgadze and V.A. Kudryavtsev, _Geant4 simulation study of low-Z material detection using muon tomography_ , _Journal of Instrumentation_ 18 (2023) C12014.
  * [15] A. Georgadze, A. Giammanco, V. Kudryavtsev, M. Lagrange and C. Turkoglu, _A simulation of a cosmic ray tomography scanner for trucks and shipping containers_ , _Journal of Advanced Instrumentation in Science_ (2024) .
  * [16] A.S. Georgadze, _Simulation study into the detection of low- and high-Z materials in cargo containers using cosmic ray muons_ , _Acta Physica Polonica B Proceedings Supplement_ 17 (2024) 1.
  * [17] A.S. Georgadze, _Automated object detection for muon tomography data analysis_ , _Journal of Instrumentation_ 19 (2024) C07004.
  * [18] Particle Data Group collaboration, _Review of particle physics_ , _Phys. Rev. D_ 98 (2018) 030001.
  * [19] G.R. Lynch and O.I. Dahl, _Approximations to multiple coulomb scattering_ , _Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms_ 58 (1991) 6.
  * [20] V.L. Highland, _Some practical remarks on multiple scattering_ , _Nuclear Instruments and Methods_ 129 (1975) 497.
  * [21] ATLAS Collaboration, _Calculation of radiation length in materials_ , _PH-EP-Tech-Note-2010-013_ (2010) .
  * [22] S. Vanini, P. Calvini, P. Checchia, A. Rigoni Garola, J. Klinger, G. Zumerle et al., _Muography of different structures using muon scattering and absorption algorithms_ , _Philosophical Transactions of the Royal Society A_ 377 (2019) 20180051.
  * [23] G. Blanpied, S. Kumar, D. Dorroh, C. Morgan, I. Blanpied, M. Sossong et al., _Material discrimination using scattering and stopping of cosmic ray muons and electrons: Differentiating heavier from lighter metals as well as low-atomic weight materials_ , _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 784 (2015) 352.
  * [24] C. Hagmann, D. Lange and D. Wright, _Cosmic-ray shower generator (CRY) for monte carlo transport codes_ , in _2007 IEEE Nuclear Science Symposium Conference Record_ , vol. 2, pp. 1143–1146, 2007, DOI.
  * [25] R. Brun and F. Rademakers, _Root — an object oriented data analysis framework_ , _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 389 (1997) 81.
  * [26] G. Anbarjafari et al., _Atmospheric ray tomography for low-Z materials: implementing new methods on a proof-of-concept tomograph_ , _arXiv preprint arXiv:2102.12542_ (2021) .
  * [27] M. Descalle, D. Manatt and D. Slaughter, _Analysis of recent manifests for goods imported through US ports_ , Tech. Rep. Lawrence Livermore National Lab.(LLNL), Livermore, CA (United States) (2006), DOI.
  * [28] W.T. AutoWeigh, _Material Bulk Density Chart material bulk density chart_ , 2024.
  * [29] Bulk-density-guide, _Bulk-density-guide_ , 2020.
  * [30] R. Hoch, D. Mitra, K. Gnanvo and M. Hohlmann, _Muon tomography algorithms for nuclear threat detection_ , _Opportunities and challenges for next-generation applied intelligence_ (2009) 225.


# X-Ray bone abnormalities detection using MURA dataset

Anna Solovoya, Igor Solovyov, TraumAI

###### Abstract

We introduce the deep network trained on the MURA dataset from the Stanford
University released in 2017. Our system is able to detect bone abnormalities
on the radiographs and visualise such zones. We found that our solution has
the accuracy comperable to the best results that have been achieved by other
development teams that used MURA dataset, in particular the overall Kappa
score that was achieved by our team is about 0.942 on the wrist, 0.862 on the
hand and o.735 on the shoulder (compared to the best available results to this
moment on the official web-site 0.931, 0.851 and 0.729 accordingly). However,
despite the good results there are a lot of directions for the future
enhancement of the proposed technology. We see a big potential in the further
development computer aided systems (CAD) for the radiographs as the one that
will help practical specialists diagnose bone fractures as well as bone
oncology cases faster and with the higher accuracy.

##  1 Introduction and Motivation

Muskuloskeletal traumas affect more than 1.7 billion people a year. Missing or
misdiagnosing pathology after X-Ray or other radiology examination is the
widespread problem that directly connected with human-factor as a result of
huge workload due to increased number of patients and shorten the time that
practical specialist has to evaluate the radiographs. These mistakes have
serious consequences and often result in delayed treatment, worse long-term
prognosis, increase of treatment costs and time, and can lead to disability of
patients. The success of detection of the bone abnormalities combines both
specialists experience and quality of the image. And while the quality of the
equipment that are used for the X-ray imagine steadily improves, human factor
in the bone fracture detection remains today an unsolved problem. There are a
lot of factors that lead to the incorrect diagnosis, such as timing, level of
experience of the radiologist and personal issues that might affect the
quality of the work. For example, according to the data from the emergency
departments there is a common peak in the biggest level of the radiographs
interpretation error between 8 pm and 2 am (47 percent error in compare to 20
percent in the control group) [1]. Another problem that lead to the
misdiagnosis is a big proportion of the young specialists, who lack the
experience to detect the complicated injuries and usually need assistance from
the more senior specialists to make the final diagnosis. Moreover, there are
some areas of the skeleton, where detection of the abnormalities is very hard
due to the limitation of the X-ray equipment and the lack of time to detect
the right bone position of the image classification. The opposite situation
when every x-ray image or other imaging examinations are assessed by 2-3 and
more specialists has a negative effect on the speed and quality of the
diagnostics. All experienced medical doctors in all specialities rely on
radiologists and described findings, but always need to analyse all images
themselves for the purpose of better understanding and better planning of
treatment. Without systems that help and, for instance, highlight pathologie,
each specialist does the same work and image-assessment from the beginning .

Understanding all the challenges that face the clinical radiology branch
today, our team decided to build the simple computer-aided diagnose system
(CAD) to detect bone fractures on different parts of the human body and, thus,
reduce the workload for the practical specialists and reduce the percentage of
the possible falls negative cases.

The field of deep technology has dramatically shifted during past 10 years and
the field of medicine is the one that may more than any other area of our life
benefit from machine learning (ML)-powered tools. Today a lot of the consumers
healthcare apps incorporated into the portable devices such as smartphones or
smartwatches already use some kinds of the algorithms that track the physical
condition of the user, store the data and create prediction about hearth
attacks, strokes and diabetics complications. But in the clinical practice the
level of using such kind of the technology solutions remains at the minimum
level. This facts can be explained by the specific of the medical field which
has a high level of safety requirements to the tools that might be used in the
practical environment, as well as the strict requirements to the safety of the
patients personal data. From the other side, the problem is that the number of
the effective solutions available to the clinics right now are too small and
usually this solutions are tailored specific for the need of the particular
hospital and cannot be applied in the other clinics. This gap between the
potential to use technological solution in the diagnoses and the current
politics of the hospitals affects negatively both the patients and practical
specialists, as the work pressure for the doctors remains too high and the
speed of the patients diagnose remains to low.

###  1.1 Contribution

The huge pressure and work-overload that faced medical specialist in the last
year showed the urgent need in development the computer aided diagnose systems
that will significantly reduce the human biased errors and help medical
specialist create the prioritization lists for the emergency patients. Despite
the big progress in the field of computer vision in the medical area, some
branches of the medical diagnostics still lack the competent CAD solutions for
the purpose of the diagnostics and further treatment. Our team applied the
current state-of-art technology to the one of the most difficult and human-
biased field of radiology - the problem of bone abnormalities detection using
standard X-ray images. Using the proposed network architecture we have
gathered good experiment results; this indicates that the further development
of the CAD system for bone abnormalities detection based on the described in
the article technology approach will lead do the shifting the quality of the
real-time diagnostics.

##  2 Previously works

Because of the shared interest of the academia and practical specialists, the
field of the medical ML solutions is rapidly growing for the last 10 years.
The quality of the algorithm predictions has significantly raised and some of
the areas of expertise are performed today better by algorithms rather by
human specialists. One of the particular area, when today algorithms are
outperforming practical specialists is image-based diagnoses, for example, ML
algorithms detects spotting malignant tumours better than experienced
radiologists and can guide doctors how to construct cohorts for costly
clinical trials [2]. Such fields of medicine as oncology, lung radiology,
cardiology and brain surgery are already benefited from advanced Ml solutions,
that used in the practical diagnostics[3] [4]. However, the better diagnoses
is not the only issue that can be solved using algorithms, foe instance,
optimizing the work load to prioritizing of the most difficult and urgent
cases, pre-analyses of the cases in the high-loaded circumstances, where human
biases such as a fatigue can lead to the negative effects and improving the
quality of the medical images - are the areas that can get a significant
support from the ML-tools incorporation[4].

It must be mentioned, that the area of the image-based diagnoses has largely
benefited from rapid development of the computer vision technology, for
instance, such initiatives as ImageNet and Kaggle helped to shift the image
recognition to the new level. But a lot of image recognшtion medical tasks
remained unsolved today. The field of trauma diagnostic is one of the example
of the medical branch that relies heavily on the image-based diagnostics and
it is one of the most promising directions of the ML-powered healthcare
solutions development. But today the CAD technology in trauma radiology still
haven’t reached the acceptable level of accuracy, which would allow to use
algorithms in the practical diagnosis.

There is a substantial amount of work on the analysis in the field of
computerized decision support systems for trauma patients. Thus, the main
problem that remains is still a lack of the quality data. However some good
results even with the small amounts of the training data have been achieved so
far. For example, in his work Gene Kitamura reached 76 per cent accuracy of
the ankle fracture predication using a small dataset (a total of 596 normal
and abnormal ankle cases) [5]. But using the bigger datasets can lead to even
better results, that’s why using the additional sources of the trained data as
well as usage of the augmentation approach to the dataset enlargement become
the crucial steps in developing the algorithm for CAD. For instance, Takaaki
Urakawa reported the accuracy of 95.5 per cent in detection of the
intertrochanteric hip fractures using the data from 1773 patients, who were
enrolled in this study[6]. Good results were also showed in detection of the
distal radius fractures - using the dataset of 4,476 augmented anteroposterior
hands X-ray images Erez Yahalomi and colleagues have achieved 96 accuracy in
the radius bone fracture prediction[7]. Chi-Tung Cheng used 25,505 limb
radiographs gathered between January 2012 and December 2017 and 3605 PXRs
gathered between August 2008 and December 2016 for the algorithm pretraining
and training to achieve 91 per cent accuracy in the limb fracture
prediction[8].

However, the size of the available training data is not the only problem that
connected with the radiographs. The quality of the images as well as picture
size requirements can differ significantly even within one country. The
attempt to reduce the X-ray radiation doses affects the quality of the images
and so called quantum noise appeared on the images, making them blurry and
difficult to understand both for human-specialists and CAD [9]. Moreover, the
absence of the universal requarements to the image quality also makes the
diagnose process harder - images can be grayscale or color, high-resolution or
low-resolution, hard copy or soft copy, uncompressed or compressed. All these
difference in the radiographs caused up to 30 per cent error rate in the
diagnosis every year, that lead to the complications in the patients
treatment[10]. So, another direction of the research and development with a
great potential for the healthcare field is the area of image processing,
which can significantly shift the quality of the X-ray images both for the
needs of human doctors and intelligent systems.

The other serious drawback of the current datasets is the absence of the
universal approaches to the radiographs among different hospitals and private
practical centers. The difference in the image specifications that used in
different clinics even withing the territory of one city makes the process of
CAD implementation very difficult. To achieve better results additional steps
such as analysis and processing of the bone fracture images should be used,
however; there is different image issues that cant be solved even after the
image processing steps.

Besides the size of the training data, the technology and chosen architecture
used in the study can affect the achieved level of accuracy. The other crucial
factor of the study results is the type of CNN that were used in the CAD. For
example, using RestNet50 for the relatively small dataset of the femure bones
fractures showed a good result of 90 per cent accuracy for the fracture
prediction [11]. Using the DensNet architecture for the hip fractures
diagnosis showed also a promising results both on the small and large
datasets- more than 90 per cent accuracy in the fracture detection using X-ray
images[12], [13]. InceptionV3 architecture was used to detetct fractures on
the wirst and femur bones with relatively small data sets (1389 images) and
showed a good result during the testing - 0.95 AUC [14].

The great diversity of the different metrics used to track the progress of the
neural networks architecture described in the studies in the field of
musculoskeletal radiographs abnormality detection should also be mentioned.
Among the most widespread metrics of success that use different authors is the
level of accuracy of prediction [15], [16], [17]. The other popular
measurements of the fracture detection success are sensitivity (the percentage
of true positive cases to the general amount of the positive predictions) and
specificity (the percentage of true negative cases to the general amount of
the negative predictions) [18], [12],[15]. As for the competitions with MURA
dataset the main index of the success of the algorithm is the Cohen’s kappa
statistic, which expresses the agreement of the model with the gold standard.
The use of Cohen’s kappa statistic is determined by the ability to provide
more valuable information in the case of musculoskeletal studies [19]. Today
the best showed Kappa score is 0.843, that was achieved in November 2018 year.

##  3 Dataset

As the source of the train and test data for our experiment MURA dataset of
the musculosceletal abnormalities was used. MURA is one of the largest public
radiographic image datasets published in Jan 2018 by Stanford ML Group. MURA
dataset consists of musculoskeletal radiographs of 14,982 studies collected
from 12,251 patients, with a total of 40,895 multi-view radiographic images.
Each radiographs belongs to one of seven standard upper extremity radiographic
study types: elbow, finger, forearm, hand, humerus, shoulder, and wrist. All
images have standard format of 8-bit png. Each study was manually labeled by
board-certified radiologists from the Stanford Hospital at the time of
clinical radiographic interpretation in the diagnostic radiology environment
between 2001 and 2012. All the images were divided into the positive (with
abnormality) and negative (without abnormality) cases and present the binary
classification problem; no additional labeling from the practical specialists
from our side on this stage was used due to the aim of the experiment.

To achieve better result through the diversity of the training and tested data
we used data augmentation, for the purpose of this stage the next steps were
made: random horizontal flip, random image rotation up to 30 degrees,image
scaling randomly in the range of 95-130 per cent, brightness changing in the
range 80-120 per cent.

##  4 Model

###  4.1 Image pre-processing

MURA dataset of the bone abnormalities consist of the large variety of the
different pictures with the different formats and sizes. Thus, one of the
first problems that our team has faced working with this dataset was the
necessity to make the data more homogeneous. That’s why for the purpose of the
research work we used image pre-processing tools to achieve better quality of
the images. The first step that we have made was to find the RoI in the image
by binary threshold as well as to find contours of the image. This procedure
allowed us to detect the region of the interest for the classification and to
crop only RoI.

Data augmentation as an important part of the algorithm learning was also used
during our research stage. To enlarge dataset we used different data
augmentation approaches such as:

  * •

Random flipping the images horizontally;

  * •

Random rotation of the images up to 30 degrees;

  * •

Scaling randomly in the range 95-130 per cent;

  * •

Randomly changing brightness in the range 80-120per cent;

Before using radiographs to the first step of the learning session they were
normalized to have the same mean and standard deviation of images in the
ImageNet training set.

After the data augmentation normalization was used to DenseNet indexes of mean
and std. This step was made to allow us change values of mean and std on the
next stages of our project.

###  4.2 Network architecture and training

The first step of our work was to choose between different architecture of
Deep learning models to find the most appropriate one for the purpose of our
project. The most widespread approaches for the image detection in the medical
file is the usage of Depp Neural Networks (DNN) and Convolutional Neural
Networks (CNN). Thus, both of them showed pretty good results in the previous
works, we decided to use CNN architecture based on two main chose criteria:
good performance of the architecture with fast model training and the
possibility to eliminate error propagation during the process of the learning
sessions. As a result, 169-layer convolutional neural network was used to
predict the probability of abnormality for each radiograph in MURA dataset.
The specifics of the chosen network is that it uses a Dense Convolutional
Network architecture to connect each layer to every other layer in feed-
forward fashion to make the optimize the deep networks results. For the
purpose of the experiment the weights of the network were initialized with
weights of the ImageNet pretrained model. Taking into account the specifics of
the MURA dataset and importance of the precise results for the medical
diagnostics, during the learning sessions the batch size was chosen to be 8
images per session. Learning rate was 1e-4 -¿ and it was multiplied every 5-7
epochs by 0.1. In average the number of epochs was 35-50 from ImageNet
pretrained DenseNet-169. In the first pretraining session denselayer with
freezed encoder was used (for first 5 epochs), then the whole network was
unfreezed and it was trained 20-30 epochs more.

Figure 1: Example of the fracture detection on the leg bone

##  5 Results

To understand the level of success of the proposed architecture normally the a
new proposed system is compared to the already existing solution with the
similar functionality. However, there are not so many working in the practical
conditions bone abnormalities detection systems, that why we decided to
compare achieved score with the results that were achieved in the similar
studies.

As the measure of success of our solution we decided to chose three basic
parameters - Kappa score, AUC ROC score and accuracy of the prediction. Kappa
score as the main measurement index of the MURA dataset challenge was included
to compare the achieved results with the results that were achieved directly
on this data. The best results based on the Kappa score index that were
achieved by own solution belongs to the wrist (0.942), hand (0.862) and
shoulder (0.735) radiographs. Thus, the result of our solution is slightly
better than the maximum indexes showed on the official challenge web-site
(0.931, 0.851 and 0.729 accordingly). AUC ROC has the index of 0.870 - this
result showed the positive progress in the algorithm prediction, however, some
issues connected with the image quality and diversity of the clinical
conditions presented into the dataset remains that need to be solved to
achieve a higher results. The accuracy of the prediction was estimated of
0.863, which is also a pretty good result for the CAD that were trained only
using MURA dataset, but we see a lot of rooms for further improvement that we
plan to work with on the next stages of our RnD work.

Graphical representation of the achieved results you can see on the Pic.1 and
Pic. 2

Figure 2: Prediction over study using 4 densenet-169 models ensemble (Cohen’s
kappa statistic)

##  6 Further Work

Despite the good results achieved on this stage of our research project there
are however some possibilities to increase the accuracy of the algorithm. The
first of them is the further improvement of the data used for training the
model. The quality and diversification of the data presented in MURA dataset
to some degree affects negatively the results of the model training. The deep
learning technique as the one that rely heavily on data, is very sensitive to
the quality of the training data. All the radiographs presented in MURA
dataset were made in clinics from 2001 to 2012 with the resolution of 600x600
pixels. The low resolution of the images almost eliminates the possibility to
detect tiny fractures without precise labeling from the radiologists with a
big practical experience. For instance, the current state of the radiograph
resolutions used in the clinics is an average of 2000x2500, which is far more
better that the ones used in MURA dataset [20]. To achieve the better results
on the next stages of our work we plan to use radiographs with the current
standard matrix-size collected from the last three years by the local clinics.

The volume of the training data affects also the level of the accuracy that
the model can achieved and solve one of the main factors of the project
success. Thus, the MURA dataset present quite a big amount of the data in
general, it gives not so many data to train the algorithm to train on the
separate body parts. To gain the better result the train dataset should be no
less than 10 000 samples, this will allow to get precise forecast about the
fracture localisation. Pic 2 reflects the proportion of the cases in the MURA
dataset according to the different parts of the body. As we can see on the
diagram the number of the cases for each segment is far too small in
comparison to the data needed to archive the credible results.

Figure 3: Number of the case studies according to each type of the fracture

The other problem that we faced during this stage is the big variety of
different clinical cases. Thus a great diversity of different types of the
abnormalities in the combination with the list of different body parts used in
the dataset creation can be a cause of the вecrease in success rates of the
algorithm. On the next stages of our research work we plan to focus on the one
type of the abnormality (fractures) and narrow down the dataset to only one
type of the bones (arm).

The absence of the labeling from the practical specialist is also an issue we
plan to eliminate on the next stages. Extension of the initial information
from a simple binary classification to the determination of a site with a
specific bone abnormality will increase the accuracy of the abnormality
prediction. To the purpose of this issue we made an agreement with the
clinical experts with the average practical experience of 10 years that will
work on the labeling of the radiographs for the next stage of our work.

##  7 Conclusion

To conclude, this article presents bone abnormalities detection model trained
on the Stafford’s MURA dataset. Muskoskelital abnormality detection is an
important phase of the clinical diagnoses, that can be affected by the human
biases such as tiredeness, big workload, lack of expertise, difficult-to-see
cases of the bone fractures and muscoskeletal abnormalities, overworking, lack
of concentration and multitasking of the radiologist.To eliminate this factors
and to make the abnormality detection process faster and more reliable
computer-aided diagnose systems should be incorporated into the working tools
of the practical specialist. Ability to identify a zone of the abnormality and
further building of the prioritization list of the patients according to the
degree of the danger of the possible diagnoses will significantly increase
such crucial factor as the speed of the speed of medical care and will
decrease the unnecessary work load for the doctors. However, such computer
aided systems cannot replace radiologists in any way, but will significantly
shifted the quality of their work.

Our solution has showed a good results on the bone abnormalities detection,
however it was focused only on the binary problem classification (existence or
absence of the abnormality). To make it useful for the practical experience
the model should be further trained to be able to distinguish different types
of abnormalities with the stable high level of the accuracy.

It should also be mentioned, that besides the technological difficulty of the
CAD system development, there is also an ethical problem, that directly
affects the further development of this field. Despite the fact that ML-
powered tools significantly reduce the work load for the practical
specialists, as well as they increase the speed of the diagnoses, the ethic
aspect of using algorithms to predict some deviations in the health condition
of the patient remains a big issue that directly affects the incorporation
these solutions into the daily practice of the radiologists. There are two
main approaches in developing ML-based medical applications - supervised and
unsupervised learning. An while, the first one is more simple an
understandable for human-operators, the second one is more promising, as it is
the one that can significantly reduce the human biases into the algorithms
work and, thus, shift the accuracy of the prediction. But the mechanism of
unsupervised algorithm can be very hard to understand and in the healthcare
branch, as the one of the most regulated ones, its almost impossible to get
the permission for its usage in the clinical conditions. This means that
further steps to the cooperation between medicine and technology should be
made to make the research projects a part of the clinical practice.

## References

  * [1] Peter Hallas and Trond Ellingsen.  Errors in fracture diagnoses in the emergency deparment - characteristics of patients and diurnal variation.  BMC emergency medicine, 6:4, 02 2006.
  * [2] Thomas Davenport and Ravi Kalakota.  The potential for artificial intelligence in healthcare.  Future Healthcare Journal, 6, 2019.
  * [3] Ian Boon, Tracy Au Yong, and Cheng Boon.  Assessing the role of artificial intelligence (ai) in clinical oncology: Utility of machine learning in radiotherapy target volume delineation.  Medicines, 5(4):131, Dec 2018.
  * [4] James H. Thrall, Xiang Li, Quanzheng Li, Cinthia Cruz, Synho Do, Keith Dreyer, and James Brink.  Artificial intelligence and machine learning in radiology: Opportunities, challenges, pitfalls, and criteria for success.  Journal of the American College of Radiology, 15(3, Part B):504 – 508, 2018.  Data Science: Big Data Machine Learning and Artificial Intelligence.
  * [5] Gene Kitamura, Chul Chung, and Barry Moore.  Ankle fracture detection utilizing a convolutional neural network ensemble implemented with a small sample, de novo training, and multiview incorporation.  Journal of Digital Imaging, 32, 04 2019.
  * [6] Takaaki Urakawa, Yuki Tanaka, Shinichi Goto, Hitoshi Matsuzawa, Kei Watanabe, and Naoto Endo.  Detecting intertrochanteric hip fractures with orthopedist-level accuracy using a deep convolutional neural network.  Skeletal Radiology, 48, 06 2018.
  * [7] Erez Yahalomi, Michael Chernofsky, and Michael Werman.  Detection of Distal Radius Fractures Trained by a Small Set of X-Ray Images and Faster R-CNN, pages 971–981.  06 2019.
  * [8] Chi-Tung Cheng, Tsung-Ying Ho, Taoyi Lee, chih chen Chang, Ching-Cheng Chou, Chih-Chi Chen, I-Fang Chung, and Chien-Hung Liao.  Application of a deep learning algorithm for detection and visualization of hip fractures on plain pelvic radiographs.  European Radiology, 29, 04 2019.
  * [9] Til Aach, Ulrich Schiebel, and Gerhard Spekowius.  Digital image acquisition and processing in medical x-ray imaging.  Journal of Electronic Imaging, 8(Special Section on Biomedical Image Representation):7–22, 1999.
  * [10] Elizabeth Krupinski.  Current perspectives in medical image perception.  Attention, perception & psychophysics, 72:1205–17, 07 2010.
  * [11] Amelia Jimenez Sanchez, Anees Kazi, Shadi Albarqouni, Chlodwig Kirchhoff, Peter Biberthaler, Nassir Navab, Diana Mateus, and Sonja Kirchhoff.  Towards an interactive and interpretable cad system to support proximal femur fracture classification.  Int J Comput Assist Radiol Surg, 15, 02 2019.
  * [12] Chi-Tung Cheng, Tsung-Ying Ho, Taoyi Lee, chih chen Chang, Ching-Cheng Chou, Chih-Chi Chen, I-Fang Chung, and Chien-Hung Liao.  Application of a deep learning algorithm for detection and visualization of hip fractures on plain pelvic radiographs.  European Radiology, 29, 04 2019.
  * [13] Justin Krogue, Kaiyang Cheng, Kevin Hwang, Paul Toogood, Eric Meinberg, Erik Geiger, Musa Zaid, Kevin Mcgill, Rina Md, Jae Sohn, Alexandra Md, Bryan Darger, Kevin Padrez, Eugene Phd, Sharmila Majumdar, and Valentina Pedoia.  Automatic hip fracture identification and functional subclassification with deep learning.  Radiology: Artificial Intelligence, 2, 09 2019.
  * [14] Daniel Kim and T. MacKinnon.  Artificial intelligence in fracture detection: Transfer learning from deep convolutional neural networks.  Clinical Radiology, 73, 12 2017.
  * [15] Seok Chung, Seung Han, Ji Lee, Kyung-Soo Oh, Na Kim, Jong Pil Yoon, Joon Yub Kim, Sung Moon, Jieun Kwon, Hyo-Jin Lee, Young-Min Noh, and Youngjun Kim.  Automated detection and classification of the proximal humerus fracture by using deep learning algorithm.  Acta Orthopaedica, 89:1–6, 03 2018.
  * [16] Jakub Olczak, Niklas Fahlberg, Atsuto Maki, Ali Razavian, Anthony Jilert, Andr√© Stark, Olof Sk√ґldenberg, and Max Gordon.  Artificial intelligence for analyzing orthopedic trauma radiographs: Deep learning algorithmsвАФare they on par with humans for diagnosing fractures?  Acta Orthopaedica, 88:1–6, 07 2017.
  * [17] Sher Ee Lim, Yage Xing, Ying Chen, Wee Kheng Leow, Tet Sen Howe, and Meng Ai Png.  Detection of femur and radius fractures in x-ray images.  In IN: PROC. 2ND INT. CONF. ON ADVANCES IN MEDICAL SIGNAL AND INFORMATION PROCESSING, 2004.
  * [18] Robert Lindsey, Aaron Daluiski, Sumit Chopra, Alexander Lachapelle, Michael Mozer, Serge Sicular, Douglas Hanel, Michael Gardner, Anurag Gupta, Robert Hotchkiss, and Hollis Potter.  Deep neural network improves fracture detection by clinicians.  Proceedings of the National Academy of Sciences, 115:201806905, 10 2018.
  * [19] A.F.M. Saif, Celia Shahnaz, Wei-Ping Zhu, and M.O. Ahmed.  Abnormality detection in musculoskeletal radiographs using capsule network.  IEEE Access, PP:1–1, 06 2019.
  * [20] Walter Huda and R Abrahams.  X-ray-based medical imaging and resolution.  AJR. American journal of roentgenology, 204:W393–W397, 04 2015\.