Equalized MR Grayscale Mapping

Department of Science and Technology

Institutionen för teknik och naturvetenskap

Linköping University Linköpings Universitet

SE-601 74 Norrköping, Sweden

601 74 Norrköping

LiU-ITN-TEK-A--08/044--SE

Equalized MR Grayscale Mapping

Anders Hagvall

LiU-ITN-TEK-A--08/044--SE

Equalized MR Grayscale Mapping

Examensarbete utfört i medieteknik

vid Tekniska Högskolan vid

Linköpings universitet

Anders Hagvall

Examinator Anders Ynnerman

Detta dokument hålls tillgängligt på Internet – eller dess framtida ersättare –

under en längre tid från publiceringsdatum under förutsättning att inga

extra-ordinära omständigheter uppstår.

Tillgång till dokumentet innebär tillstånd för var och en att läsa, ladda ner,

skriva ut enstaka kopior för enskilt bruk och att använda det oförändrat för

ickekommersiell forskning och för undervisning. Överföring av upphovsrätten

vid en senare tidpunkt kan inte upphäva detta tillstånd. All annan användning av

dokumentet kräver upphovsmannens medgivande. För att garantera äktheten,

säkerheten och tillgängligheten finns det lösningar av teknisk och administrativ

art.

Upphovsmannens ideella rätt innefattar rätt att bli nämnd som upphovsman i

den omfattning som god sed kräver vid användning av dokumentet på ovan

beskrivna sätt samt skydd mot att dokumentet ändras eller presenteras i sådan

form eller i sådant sammanhang som är kränkande för upphovsmannens litterära

eller konstnärliga anseende eller egenart.

För ytterligare information om Linköping University Electronic Press se

förlagets hemsida

http://www.ep.liu.se/

Copyright

The publishers will keep this document online on the Internet - or its possible

replacement - for a considerable time from the date of publication barring

exceptional circumstances.

The online availability of the document implies a permanent permission for

anyone to read, to download, to print out single copies for your own use and to

use it unchanged for any non-commercial research and educational purpose.

Subsequent transfers of copyright cannot revoke this permission. All other uses

of the document are conditional on the consent of the copyright owner. The

publisher has taken technical and administrative measures to assure authenticity,

security and accessibility.

According to intellectual property law the author has the right to be

mentioned when his/her work is accessed as described above and to be protected

against infringement.

For additional information about the Linköping University Electronic Press

and its procedures for publication and for assurance of document integrity,

please refer to its WWW home page:

http://www.ep.liu.se/

Abstract

MR images are crucial for health care today and the application areas are con-tinuously increasing. A major problem regarding visualization of these data sets is that there is no absolute scale for the data values. Even for the same type of examination the scale varies from patient to patient, from time to time, and from scanner model to scanner model.

This thesis addresses the challenge of automatically generating visualization parameters for these data sets, eqalizing visual appearance and interactions to diversities in data set distributions. Main objectives are to achieve reasonably good starting visualizations and to ensure consistent interaction behaviour.

Acknowledgments

I would like to give special thanks to my supervisor Magnus Ranlöf for constructive discussions about ways to extend the thesis work, as well as a good introduction to Sectra’s current visualization pipeline and development routines. I would also like to thank Claes Lundström, whose Ph.D thesis has been the ground of the tissue evaluation algorithms, and Erik Sjöblom, for helpful hints about implementation issues at Sectra.

Thanks also to my mother Gunilla Hagvall, opponent Ajden Towfeek and ex-aminer Professor Anders Ynnerman for valuable thesis feedback.

Contents

1 Introduction 1 1.1 Problem Description . . . 1 1.2 Project Description . . . 2 1.3 Recommended Prerequisites . . . 2 2 Background 3 2.1 Medical Imaging . . . 3 2.1.1 Clinical Use . . . 3 2.1.2 Modalities . . . 4

2.2 Magnetic Resonance Imaging . . . 5

2.2.1 Spin Physics . . . 5 2.2.2 Principle . . . 6 2.2.3 Imaging . . . 9 2.2.4 Weighting . . . 10 2.2.5 Image Contrast . . . 12 2.2.6 MRI vs. CT Radiation . . . 12 2.3 Scientific Visualization . . . 13 2.3.1 Intensity Mapping . . . 14 2.3.2 2D Slice Viewing . . . 14

2.3.3 Direct Volume Rendering . . . 16

2.3.4 Window Level Interaction . . . 17

2.4 Statistical Tools . . . 18

2.4.1 Histograms . . . 18

2.4.2 Statistical Properties . . . 19

2.5 Tree Structures . . . 21

2.5.1 Quad Trees . . . 21

3 MR Grayscale Mapping Challenges 25 3.1 MR Uncertainties . . . 25

3.1.1 Lack of Absolute Scale . . . 26

3.1.2 Weighting . . . 27

3.2 Tissue Detection . . . 27

3.2.1 Range Weights . . . 28

3.2.2 Neural Networks . . . 28 ix

3.2.3 Wavelet Features . . . 28

3.3 Spatial Coherence . . . 29

3.3.1 Alpha Histograms . . . 31

3.3.2 Partial Range Histograms . . . 32

4 Equalized MR Grayscale Mapping 35 4.1 Subdivision . . . 36

4.1.1 Fixed Splitting . . . 36

4.1.2 Adaptive Splitting . . . 37

4.2 Tissue Evaluation . . . 39

4.2.1 Statistical Property Weighting . . . 39

4.2.2 Alpha Histograms . . . 41

4.2.3 Partial Range Histograms . . . 42

4.3 Window Level Generation . . . 43

4.3.1 Percentiles . . . 43

4.3.2 Weighting . . . 44

4.3.3 Mapping . . . 45

4.4 Window Level Interaction . . . 47

4.4.1 Absolute . . . 48

4.4.2 Percentile . . . 48

4.4.3 Region of Interest Tool . . . 49

5 Results 51 5.1 Tissue Evaluation . . . 51 5.1.1 Robustness . . . 51 5.1.2 Speed . . . 54 5.2 Subdivision . . . 55 5.2.1 Accuracy . . . 55 5.2.2 Speed . . . 57 5.3 Mapping . . . 59 6 Conclusion 63 6.1 Results . . . 63 6.1.1 Radiologist Feedback . . . 63 6.1.2 Final Product . . . 64 6.1.3 Implementation . . . 68 6.2 Clinical Limitations . . . 68 6.3 Future Work . . . 69 Bibliography 71

Chapter 1

Introduction

This thesis’ overall research topic lies within Scientific Visualization, or more fo-cused Grayscale Mapping for MRI (Magnetic Resonance Imaging). Lots of re-search effort has been put into this topic, but science and usage are two sides of the same coin, and have very different preferences when it comes to the result. For a scientist, advanced mathematical models can be used to an unlimited extent, since the goal is the result, whereas for a radiologist, the goal is to work in a fast environment with intuitive interaction possibilities.

This thesis work takes place on both sides, demanding result through the tech-nical point of view of a Master’s Thesis, and speed through the usability point of view of implementing the results in Sectra’s Image Display System application IDS7. This introductory chapter gives an overview to the project, describing the problematic situation and initial aims. Reader prerequisites are also stated.

The outline is to give background information about the thesis subject in chap-ter 2, followed by a more in depth focused chapchap-ter 3, explaining visualization chal-lenges, recent research topics considered in the implementation of this thesis work and also recent work in the area of MR Grayscale Mapping. Chapter 4 presents the implementation of this thesis, combining tools from chapter 2 and 3 into the solution implemented in IDS7 at Sectra, proceeding into chapter 5 and 6, address-ing implementation results, radiologist feedback and conclusions from and around the thesis.

1.1

Problem Description

MR images are crucial for health care today and the application areas are con-tinuously increasing. A major problem regarding visualization of these data sets is that there is no absolute scale for the data values. Even for the same type of examination the scale varies from patient to patient, from time to time, and from scanner model to scanner model. When performing the diagnostic review of these images the radiologists are forced to manually adjust the grayscale mapping for each data set, which is enormously time-consuming, especially considering that the interaction might be difficult.

1.2

Project Description

The project aims at finding a solution to the problem described above, and should be carried out in close cooperation with radiologists specializing in MR having problems with Sectra’s current solution. There are basically two problems origi-nating from the fact that MR images lacks an absolute scale, thus there are two objectives with this project:

• To achieve automatic generation of window level presets that provides a reasonably good starting point for the grayscale window.

• To ensure that manual adjustments have a consistent behavior. For instance, a certain mouse movement should correspond to a fixed relation to the dy-namic MR range instead of a fixed change of the grayscale window.

1.3

Recommended Prerequisites

In order to understand the contents of this thesis, the reader should preferably have an engineering background with knowledge in Image Processing and Analysis. Understanding of Fourier Transforms, familiarity with the Frequency Domain and basic knowledge in Electromagnetism is preferable to appreciate the Background about MRI, but in less need for the latter chapters.

Regarding the problematic situation this project aims to solve, challenges and impact of the solutions will also be addressed from a radiologist’s point of view, cre-ating a non-theoretical path were understandability is limited only by the reader’s ability to imagine working as a radiologist.

Chapter 2

Background

Although this thesis focuses on visualization and clinical use of Magnetic Reso-nance Imaging, in order for the reader to fully appreciate its contents, an introduc-tion to the more broad area of Medical Imaging is needed. This beginning chapter provides an overview of the subject, and explains the basic techniques needed to understand the meaning of the challenges, and impact of the solutions presented.

2.1

Medical Imaging

Referring to the techniques used to create images of the human interior body, med-ical imaging has been a continuously increasing area of research ever since Wilhelm Conrad Röntgen discovered X-rays in 1895 [12]. Analogically, these examinations have been printed on plastic paper to be viewed in front of an area light source, but over the last 15 years however, the possibility of digitizing medical imaging has been developed, and nowadays there are lots of film free hospitals. There are also several other modalities, or other available apparatuses than X-ray scanners to create medical imaging data sets. A Magnetic Resonance (MR) Scanner is one of them, producing data sets for MR Imaging (MRI), in which category this thesis work falls. The following sections describe the most frequently used modalities along with clinical use and workflow of medical imaging at hospitals.

2.1.1

Clinical Use

The purpose of medical imaging in clinical use is to reveal, examine or diagnose diseases, and the practitioner responsible for interpreting the data is a radiologist. At a hospital, examination requests are sent to the radiology department, including diagnostic questions to be answered. According to these questions, a suitable modality is chosen, through which the data set is then acquired. A radiologist interprets the data and writes a report back to the requester.

At film free hospitals there is a data management system, also known as PACS (Picture Archiving and Communication System). The PACS is the hub in the

Figure 2.1. Slice of a head scan, CT (left), MR (right).

digitization of medical imaging, providing visualization tools for these data sets, handling storage and distribution, replacing film archives and light cabinets.

The benefits of using a PACS are revolutionary to medical imaging [6]. Faster examination pipeline, unlimited access to examinations, less risk of losing or dam-aging images and also abilities to physically interact with the data set by zooming, panning, measuring or with the visualization by changing brightness, contrast, transfer functions etc. Sectra is a provider and developer of a PACS and a side of this thesis work is for the algorithms to run in Sectra’s current development application IDS7.

2.1.2

Modalities

In medical imaging, a modality is an apparatus from which it is possible to obtain a medical image data set. There are many modalities available at hospitals today, acquiring the data in various different ways. Several of these make use of X-rays, sending ionizing radiation into the body to measure the shadow or absorption on the opposite side. Computed Radiography and Direct Radiography are 2D imaging modalities making use of this technique. Another is Computed Tomography (CT), where the X-ray source and detector rotates around and moves along the patient, measuring the X-rays passing through it. In this case tissues don’t shadow each other, making 3D volumetric data sets possible, where each value describes the attenuation at a single point in space [6]. Although CT and other X-ray modalities have developed, a major drawback is the dangerous radiation dose [3].

MRI is based on nuclear magnetic resonance (NMR), and is to this day known not to have any harmful effects on the patient, citehornak:2003. An in depth description of MR is given in 2.2. Figure 2.1 shows a head scan captured with CT and MR respectively. The amount of noise is significantly higher in the MR scan, but also the amount of detail for soft tissue, brain substance for instance.

Ultrasound is another technique without negative side effects, [6], which mea-sures the tissue reflectance of high frequency sound waves in the interval 2-10 megahertz to produce 2D images. It is not as accurate and detailed as CT or

2.2 Magnetic Resonance Imaging 5

MR, but it is very fast and it is possible to study moving structures in real time, whereas CT and MRI are slow in the data acquisition.

Nuclear Imaging adds another dimension to medical imaging, displaying phys-iological activity instead of reflectance. Radioactive substances which are more readily absorbed by biologically active regions, such as fracture points or tumors, are administered to the patient, and then gamma monitors are used to produce images. There are 3D modalities employing Nuclear Imaging, such as Positron Emission Tomography (PET). Data acquired through PET are often viewed along with CT scans of the same area, [6]. This way, diseases detected by the PET scan and tissue detected by the CT scan can be viewed simultaneously.

2.2

Magnetic Resonance Imaging

MRI is an imaging technique originally invented by Felix Bloch and Edward Purcell in 1946 [3], but it wasn’t until more than 20 years later, when Raymond Damadian made revolutionary findings about differences between healthy tissue and tumors [2], that MRI became interesting enough to motivate the expenses it brought, and started to develop. The technique is based on principles of NMR, but due to the negative associations to the word nuclear, N was dropped leaving only MR combined with I for imaging when naming. To help understanding this section, some physics and properties of spin are brought to the subject.

2.2.1

Spin Physics

The essence of NMR when used to create medical image data sets is the property of spin. Spin is a fundamental part of nature, just like electrical charge or mass, and it comes from individual unpaired electrons, protons or neutrons [3]. Two particles with opposite spin can pair up to eliminate its presence, making only individual unpaired particles measurable.

The property of interest in NMR is that spin can be thought of as a tiny magnet moment vector, causing the particle to behave like a magnet. Therefore, when the particle is placed in a magnetic field, the spin vector will align itself with the field. There are two configurations in which the poles of the magnetic field and the particle with spin can align - one low energy state with north-south, north-south alignment, and one high energy state with north-north, south-south alignment, see figure 2.2.

In order for the particle to translate between the low and high energy state, it must either absorb the energy difference through a photon to go to high energy state, or release energy and drop to the low energy state. The photon frequency needed is related to the magnetic field strength B0 and gyromagnetic ratio γ of

the particle as Eq. 2.1. The energy of a photon is related to its frequency f and Planck’s constant h as eq. 2.2. Combined, these define the particular difference between high and low energy state as eq. 2.3 [3].

Figure 2.2. Spin alignment configurations in a magnetic field. Low energy state N-S,

N-S (left) and high energy state N-N, S-S (right).

E = h · f (2.2)

∆E = h · γ · B0 (2.3)

In NMR and MRI, the frequency is called the resonance frequency or Larmor frequency, and when a photon of that frequency hits the particle, carrying the energy difference between the particle states, energy will be absorbed and the particle will jump to the high energy state [1]. The gyromagnetic ratio differs from particle to particle, inducing different Larmor frequencies for different particles in the same magnetic field.

Basically all elements in the periodic table have isotopes with a non-zero nu-clear spin, but NMR can only be performed on isotopes with high prevalence. In the human body, the best suitable element is Hydrogen H, since we have a lot of it. All Hydrogen isotopes have non-zero nuclear spin, but the natural abundance, or fraction of Hydrogen being1H is 99.985 %, making a study of the other isotopes impossible. The magnetic field strength varies from MR Scanner to MR Scanner, but overall Larmor frequencies typically lie in the 15-80 MHz interval, which cor-respond to radio frequencies [3]. Different tissues have different1H density, and simplified, MRI is a method of measuring these changes.

2.2.2

Principle

An MR Scanner consists of a great bore magnet creating a static magnetic field B0

in the Z-direction, see figure 2.3. This is the field in which the spin of the Hydrogen nuclei aligns. The net magnetization in the system has two components, one along the z-axis, corresponding to particles in the low energy spin state, and one perpendicular to the Z-axis, i.e. in the XY -plane, corresponding to particles in the

2.2 Magnetic Resonance Imaging 7

Figure 2.3. Bore maget placement and static magnetic field direction Z.

Figure 2.4. Magnetization components at equilibrium. Left, low energy state particles

aligning spins along the Z-axis. Center, high energy state particles contribute orthogonal to the Z-axis, however incoherent phase in the XY -plane causes the contributions to cancel out, resulting in pure Z-directional net magnetization (right).

Figure 2.5. Net magnetization behaviour during its return to equilibrium. Left, the

radio pulse causes particles to translate to high energy state, rotating along with the radio pulse at angular frequency ω = 2πf , f being the Larmor frequency. This induces coherent XY -phase and a net magnetization in the XY -plane. When the radio signal is turned off (center), particles return to the low energy state and a dephasing in the XY -plane starts. The net magnetization returns to equilibrium (left), but faster in the XY -plane due to the dephasing, compare T1and T2.

high energy state. The property of spin will make the low energy state alignment fixed, but rotating around the z-axis at the Larmor frequency for particles in the high energy state. However, there is no coherent pattern in the rotation at equilibrium, causing the net magnetization in the XY -plane to cancel out, see figure 2.4.

Radio waves at the Larmor frequency corresponding to B0 and the

gyromag-netic ratio γ of1_{H is sent into the system, creating the effect of particles translating}

to the high energy state, thus changing the net magnetization. In the Z-direction, it changes since less particles remain in the low energy state, and in the XY -plane, we introduce a net component. This is because the radio frequency signal is also rotated at the Larmor frequency, and particles in the high energy spin state align with the radio waves, inducing coherent spin behavior in the XY -plane. When this radio frequency signal is turned off, the net magnetization will return to its equilibrium. During this time, there are two measurable things happening, which are used to create the MRI data sets, see figure 2.5. [3]

First, the time it takes for the net magnetization along the Z-axis to recover is interesting. This is called the spin-lattice relaxation time, and it is denoted T1.

Eq. 2.4 governs the behavior of spin-lattice relaxation as a function of time t.

Mz= Mz0(1 − e

−t

T1_{) (2.4)}

Second, the loss of phase coherence in the XY -plane is measurable. This causes the net magnetization to decrease, since the spins loose their coherence, and start to cancel the net XY -component out. The time it takes for this to return to equilibrium is called the spin-spin relaxation, and is denoted T2. Eq. 2.5 is the

corresponding equation for spin-spin relaxation as a function of time t.

Mxy= Mxy0e

−t

2.2 Magnetic Resonance Imaging 9

Figure 2.6. Gradient coils in X, Y and Z directions, used to disturb the static magnetic

field.

There is a third parameter also commonly used, which is T₂∗. This parameter is a combination of T2 and the static field B0 not being completely homogenous

and it is defined in Eq. 2.6: [3]

1 T∗ 2 = 1 T2 + 1 T2inhomo (2.6) The relationship between the three parameters is always T1> T2> T2∗, which

shows that the spin-spin relaxation is delayed because of B0not being homogenous,

since T2∗ is shorter than T2, and the loss of phase coherence makes the process of

recovering the net magnetization faster in the XY -plane than along the Z-axis, since T1> T2.

In order to convert these parameters into a medical image data set, MR Scan-ners are equipped with receiver coils. While the radio signal is turned off and the magnetic field is returning to its equilibrium, the magnetization produces an oscillating magnetic field, which emits energy at rates received by these coils. This signal is called a Free Induction Decay (FID), and it contains information about the scanned material [3].

2.2.3

Imaging

The signal we get from the process explained in 2.2.2 is however still no good to us, since we have no idea from exactly where it is coming. The static magnetic field affects the entire body, and there are lots of possible locations from where a received signal could have come. This induces a need to control the magnetic field. By making the magnetic field only affect certain parts of the body, for example a very thin slice, it would be possible to put the received signal in the correct location, thus being able to interpret it and form an image out of it.

The solution to this is to equip MR Scanners with magnetic gradient coils [3]. There are three coils, one for each dimension in space, and they are used to create a magnetic field gradient, i.e. to disturb the static magnetic field, see figure 2.6. If different regions of the body experience a different static magnetic field, different Larmor frequencies also apply to the different regions, creating the possibility to focus the signal on a particular region, for example a very thin slice of the body. Since there are coils for all dimensions, the orientation of a slice is completely free. The next step in the current MRI process was individually found by two scien-tists in 1983, Ljunggren and Tweig [16], and made a number of seemingly complex

tasks easy, such as interpreting the received signal. They both showed that the modulated MR signal received at the receiver coil equals the Fourier transform of the spin density, i.e. the spatial frequencies at their orientations in the Fourier do-main, also called k-space in MRI. The magnetic gradient coils control the position k = (kx, ky, kz) in k-space, and then the signal, or spin density, can be expressed

as Eq. 2.7:

S(t) = ˆρ(k(t)) = Z

ρ(r)ei·k(t)·rdr (2.7)

With r being the position in image space, ρ being the spin density, ˆρ its Fourier Transform and S being the received signal. The Fourier Transform is a well known transformation between the spatial domain and the frequency domain. It is pre-sented here for comparison, and a formal proof can be found in [11] as an example. The continuous 1D Fourier Transform and its inverse are defined in Eq. 2.8 and 2.9 respectively. ˆ f (ω) = ∞ Z −∞ f (t)e−iωtdt (2.8) f (t) = 1 2π ∞ Z −∞ ˆ f (ω)e−iωtdω (2.9)

To solve Eq. 2.7 for ρ, thus extracting intensities for slice coordinates complet-ing the task of obtaincomplet-ing a data set of the imaged anatomy for display, Eq. 2.9 is applied in 2D, transforming the frequency distribution into spatial information about spin density, i.e. intensities at pixel position.

2.2.4

Weighting

The term weighting originates from the fact that we want to emphasize contribu-tions from T1 or T2 in the received signal. Certain tissues may have similarities

in either of T1 or T2, thus making them hard to distinguish from each other [3].

This is done by using different pulse sequences when sending radio waves into the system. Not only does it allow emphasizing T1 or T2 contribution, it is also good

to repeatedly receive what should be the same signal in order to reduce signal to noise ratio (SNR), [3].

There are several sequences used in MRI, which all serve particular purposes. The differences between them lie in how and when the radio frequency waves are sent into the system. Extractions of data from the spins always work as described in the previous sections. There is no need to go deep into this subject, but a more in depth explanation can be found in [3].

The regular way described before corresponds to a sequence called the 90-FID sequence, see figure 2.7. Simplified, the signal received doesn’t have any T2

2.2 Magnetic Resonance Imaging 11

Figure 2.7. Top, a 90-FID sequence. The radio signal rotates the net magnetization

90 degrees into the XY -plane. During recovery, a FID is received depending on T1 and

the repetition time T R, see eq. 2.10. Bottom, a spin-echo sequence, where an additional radio pulse is used to rotate the net magnetization 180 degrees, down the negative Z-axis. During return to equilibrium, an echo of the FID is created, introducing T2 and

dependencies, since what we receive is just how the magnetization recovers, and not any information about the dephasing. Eq. 2.10 governs the 90-FID signal related to spin density ρ, T1and the repetition time T R.

S = ρ(1 − e−T RT1 _{) (2.10)}

A more versatile sequence to use is the spin-echo sequence. The FID is not considered here, but instead we apply a second radio signal during the relaxation time, a radio signal designed to invert the magnetization along the Z-axis. This doesn’t affect the Z-recovery, but instead partially causes the XY -magnetization to rephase, and produce an echo, see figure 2.7. This echo reaches is maximum intensity T E seconds after the sequence started, and it has both T1 and T2

depen-dencies as: S = ρ(1 − e −T R T1 _)e −T E T2 _(2.11)

There are three typical ways of alternating T R and T E, thus there are three kinds of weightings associated with spin-echo sequences, T1, T2, and ρ weighted

images. As can be seen by studying Eq. 2.11, using a T R much greater than T1 will cause the parenthesis to approach 1, eliminating the effect of T1. Using

a very quick T E will cause the T2 effect to vanish, as the term depending on T2

approaches 1 for such a configuration. From the same equation it’s easy to see the opposite effect, i.e. when T1 and T2 will contribute more to the solution. Table

2.1 show T R and T E values to use in order to emphasize ρ, T1or T2 for a certain

tissue.

Weighting T R T E

ρ  T1  T2

T1 ≤ T1  T2

T2  T1 ≥ T2

Table 2.1. T R and T E values to use for different weightings.

2.2.5

Image Contrast

Since we know the parameters controlling the intensity of the signal, we can adjust these to fit the needs of each particular examination. T1, T2and ρ values are known

to a certain extent for each tissue, and therefore Eq. 2.11 can be used to get data sets weighted towards the tissues of interest. Figure 2.8 demonstrates the results of varying T R and T E to achieve T1 and T2 weighted images respectively.

2.2.6

MRI vs. CT Radiation

A relevant question rising from this chapter is where the fact that MRI has no known negative side effects on humans comes from, since radiation doses are clearly

2.3 Scientific Visualization 13

Figure 2.8. Two MR images scanned using a spin-echo sequence. Left, T R = 500ms,

T E = 20ms, T1 weighted image. Right, T R = 2000ms, T E = 80ms, T2weighted image.

as big as for CT. This orginates from the fact that the radiation lie in radio frequencies, which are not classified as ionizing radiation. As an example, let’s look at the energy ratio of an X-ray photon and a typical MRI photon. The static magnetic field strength B0in an MR Scanner may be 1.5T , and the gyromagnetic

ratio γ of Hydrogen is 42.58M hz_T , correspondsonding through Eq. 2.1 and 2.2 to a Larmor frequency f and energy E of:

f = γ · B0= 42.58 · 106· 1.5 = 63.87M Hz

E = h · f = 6.662 · 10−34· 63.87 · 106_{= 4.32 · 10}−26_J

For an X-ray photon of frequency 20EHz the energy then becomes E = 1.33 · 10−14J . The X-ray photon carries 1.33·10_4.32·10−14−26 = 3.14 · 1011 times the energy of the

MRI photon, telling us that X-rays are of order 1011 _{times more ionizing than}

radio photons, which is why MRI is not considered hazardous to patients.

2.3

Scientific Visualization

Scientific Visualization is the process of transforming data set intensity values into pixel intensities on screen. As explained earlier, see 2.1.2, there are several modalities capable of producing 3D data sets, however the highly dominant method is to view the 2D slices in the format they were delivered from the modality, i.e. by scrolling through a stack of images, interacting by zooming, panning and changing the intensity mapping to highlight areas of interest in the images. As explained in section 2.2.3, MRI has no orientation limitations, meaning that the slices are probably already delivered in a thought through orientation, however it is still a drawback to not be able to manually adjust this during visualization.

There are methods employing this interaction possibility, one being Multi-planar Reconstruction (MPR). A slicing plane of arbitrary orientation and curve is

selected, and then intensities are interpolated from the closest slices to the positions in space corresponding to the selected slicing plane. There are also methods to visualize entire 3D volumes, providing even more user freedom, with Direct Volume Rendering (DVR) being the most versatile and relevant for diagnostic work. The problematic situation this thesis aims to solve lies within grayscale mapping for 2D slice viewing, but adding a third dimension would only improve the algorithms, to a cost of more computations. To provide a wider view of the subject, DVR also has a paragraph in this section.

2.3.1

Intensity Mapping

Simplified, this is the central step in the visualization pipeline. When a viewing plane and for 3D additional settings such as camera position and orientation have been specified, the next step is to find a transition between the intensity values in the data set and the color output for display, intensity being the weighted spin density for MRI, as explained in section 2.2.3 and 2.2.4.

Computer displays today are capable of outputting 28_{different intensity levels}

in three different channels - red, green and blue, for a total of (28₎3 _{= 2}24 ₌

16777216 different colors. However, since the information in the data set consists only of one component, it’s impossible to connect the different color channels to different visual dependencies, limiting the use to one intensity channel, or 28= 256 different intensity levels. Using all three channels simultaneously will result in an intensity scale from black to white, also referred to as the grayscale.

MR Scanners usually produce 12-bit data sets, meaning we have 212 _{= 4096}

different possible intensity levels, but both higher and lower bit representations exist, as well as signed and unsigned data, meaning for a 12-bit data set that the intensities either range from [0, 212_{] or from [−2}11_{, 2}11_{]. Using a uniform mapping}

of intensity values is of course not possible, since that would mean displaying the intensity level 127 equal for a signed 8-bit and unsigned 16-bit data set. 127 is the highest possible intensity, 100.0% for the signed 8-bit case whereas for the unsigned 16-bit set, 127 represents only 127₂16 =

127

65536 = 0.19% of the possible range.

Because of properties of MRI, the amount of the possible scale actually occu-pied by intensities varies a lot, so there will be cases where an 8-bit and 16-bit data set could share grayscale mapping. However, this only adds to the fact that a uniform mapping is impossible, since it also tells us that two different 16-bit data sets might need very different mappings. Uncertainties in MRI will be explained further in 3.1.

Mapping of intensity values become very different in the 2D and 3D case, since in 2D the image plane is specified, and there can be no tissues on top of each other, clouding one or the other’s appearance.

2.3.2

2D Slice Viewing

With image plane specified, thereby having the data set to be visualized narrowed down from 3 to 2 dimensions, the transition between intensity and display output only needs to map the intensity range to a color range. The grayscale preferably

2.3 Scientific Visualization 15

Figure 2.9. Example of a window level function.

Figure 2.10. Three different renderings and corresponding window levels. Left, a min/max window mapping the minimum intensity and maximum intensity as 0 and 255, i.e. c = 500, w = 1000. Center, a narrow window in the tissue range, c = 750, w = 250. Right, c = 750, w = 1 providing white for intensities above center and black for intensities below.

Figure 2.11. Example of an optical model.

since we don’t have any multi-channel coloring information about the anatomy. This is done with a function called the window level.

The simplest and almost unanimously used form of a window level is a linear function, mapping intensity values to the display output range [0, 255] in a linear fashion. These window levels are described using the terms center and width, corresponding to the intensity mapped to half the display output range and the entire range of intensities mapped with contrast respectively. Intensities out of the contrast range are clamped to either the max display output for intensities to the right or to the min display output for intensities to the left of the contrast range, see figure 2.9.

The goal for a window level is to provide contrast in the intensity range corre-sponding to interesting tissue in the anatomy. This is found very tough though, also due to uncertainties, dealt with in 3.1. Figure 2.10 show some examples of different window levels presented along with resulting renderings of a fixed MR data set.

2.3.3

Direct Volume Rendering

DVR is a process of generating images directly from a volumetric data set. In the 2D case we have a slice of intensities to display, and positioning just corresponds to fitting the data set values into a position corresponding to screen resolution and image window on the display. In 3D however, in order to know which volume elements, also known as voxels, corresponding to each display pixel, an optical model has to be evaluated, see 2.11. This can be thought of as simulating a camera at a certain position and orientation with objective properties in space

2.3 Scientific Visualization 17

Figure 2.12. Examples of DVR renderings for a consistent data set using different

camera positions and applied transfer functions.

taking a picture of the 3D data set. The resulting picture is then viewed on the display, to compare to a slice in 2D slice viewing.

A major difference from 2D slice viewing though is that tissues are clouding each other in the view of the camera, making the intensity mapping used for 2D slice viewing impractical. Instead DVR make use of multi-channel coloring and transparencies, specified by transfer functions. A transfer function can be thought of as similar to the window level, except each transfer function is associated to a color, and maps intensities to transparencies [5]. When calculating the display output, a version of the volume rendering integral is evaluated, for each pixel position, summing up contributions from the corresponding voxels in the data set through ray casting, [10], blending different colors together into a display output using a technique called alpha blending [6].

Through these transfer functions, the possibility to completely shut out inten-sity ranges out of interest is made possible, see figure 2.12.

DVR has nothing to do with the scope of this thesis; however it is good for a total understanding of the subject, and several techniques and algorithms used and explained originates from and could easily be used for 3D purposes as well. For a more detailed description of DVR and techniques such as ray casting and alpha blending, see [10].

2.3.4

Window Level Interaction

Another side of this thesis is to improve interaction possibilities with these visual-izations. Interaction for the 2D case corresponds to manually changing or tuning the window level to better be suited for the interesting tissue in the data set. The traditional way, currently used in Sectra’s IDS7 application is to make an either absolute or relative change in the window level. An absolute change would be letting an interaction correspond to a fixed, absolute change in either center or width for the window level. A relative change would mean to interact with a factor used to scale the window level, introducing a dependency to its previous value.

Interactions in the 3D case are not considered here, but would correspond to manual adjustments in or among the transfer functions to induce better image contributions from interesting tissues, or to separate tissues using different transfer functions instead of utilizing one to cover the range of two different tissues, making them tough to distinguish.

Figure 2.13. Grayscale image and corresponding histogram, displaying numbers of pixels for each intensity.

2.4

Statistical Tools

In order to understand the latter chapters in this thesis, addressing challenges around MR Grayscale Mapping and solutions to the problems described in 1.1, some basic statistical concepts are brought to the subject. If the reader feels comfortable in properties of data set statistics and use of histograms, this section is only going to bring concepts already of his/her knowledge to the subject.

When analyzing a medical image data set computationally, such as weighted spin densities from an MR scan, very little can be figured out just by looking at intensities corresponding to pixel positions. For humans, having a very intelligent vision system, it’s easy to put an image in a context and distinguish different tissues from each other, for example to look at a head scan from above and use common sense and experience to figure out that substances a little bit inside the outline of the head must correspond to brain tissue. If trained and experienced in anatomy, such as a radiologist, the task becomes trivial, along with several other, more advanced identifying tasks through use of recognition.

In order for recognition and data set classification to run digitally, there is a need to find properties representing a data sets in absolute numbers, thus having something to compare different data sets by. These absolutes and rules are created through histograms and statistical properties.

2.4.1

Histograms

In a statistical sense, a histogram is a tool to graphically display tabulated data, divided into uniform intervals, also called bins. For a set of collected data, there are usually several ways to make use of histograms. For a group of 100 people, age could be divided in several intervals, and then display people count for each of these age categories. Another possibility would be to split height or weight into intervals, and display people count for these. Both would be correct histgorams.

In the mathematical sense, histograms map the number of observations falling into k disjoint, uniform bins xi, with x1...xk∈ [xminxmax]. When used in imaging,

2.4 Statistical Tools 19

and use one bin for each particular intensity. A histogram function H(N, x) then maps the instances of value x in the set of intensities N , see eq. 2.12. |S| is the cadinality of set S.

H(N, x) = |N ∩ x| (2.12)

Histograms are of great interest in visualization and design of window level, since they provide information about the distribution of intensities; absolutes to use for recognition and classification [6] [3]. An image with corresponding his-togram of displayed intensities is shown in figure 2.13.

2.4.2

Statistical Properties

Data sets, or mathematical distributions in general have properties which are going to be used a lot further on in this thesis. A couple of them are presented here, along with correspondence to histograms.

Mean Value

The mean value of a distribution is defined as the average observated value, in imaging the average intensity, see eq. 2.13.

xmean= 1 n n X i=0 xi (2.13)

The mean value is found by summing every element’s intensity xi together

and divide by the total number of elements n, however if we have a histogram of the distribution, the mean value corresponds to the weighted sum of each bin, or represented intensity xi, multiplied by its count H(xi), see eq. 2.14.

xmean= 1 n k X i=1 xi· H(xi) (2.14)

Unless the relationship xmax− xmin ≥ n is true, x being observated intensities

and n the number of elements in the set, the mean value is less computationally costly to extract through a histogram, since it would correspond to less operations. The mean value gets influenced by every single element in the data set, and extreme outliers sometimes make the mean value a bad choise to represent the data set as a whole. Consider a multiset {1, 1, 2, 2, 2, 3, 14}. The mean value is 3.57, which is well over the majority of the numbers in the set. For a medical example, think of two MR images of a knee, one being natural and one with gathered fluid due to some abnormality, occupying 5% of the amount of elements in the set. Fluids contain huge amounts of Hydrogen, which will introduce extreme intensities in the data set. These two will have very different mean values, even though 95% of the data set will be almost identical.

Median Value

The median is defined as the number separating the higher and lower half of an arranged population or distribution, leading the median value in a distribution to be the observation at the median number position. For a data set, this would correspond to the intensity value separating the higher and lower half of the data set. The math befind the median orginates from probability theory, but with histograms, it can be expressed as the intensity, or bin xmedian of the histogram

satisfying inequality 2.15, H(i) and n being the histogram value at bin i and total number of elements respectively.

xmedian−1 X i=0 H(i) < n 2 ≤ xmedian X j=0 H(j) (2.15)

Recall the examples from 2.4.2, brought to the subject to express a drawback about the mean value. The median value of {1, 1, 2, 2, 2, 3, 14} is 2, and would have stayed the same no matter how low the beginning or high the ending values would have been. Same thing in the MR knee case, where the high intensities from the fluid would not have affected the median value, showing its robustness to outliers.

Percentiles

For an image, a percentile is defined to be the value below which a certain percent of the intensities fall. It is related to the median in the sense that the median actually corresponds to the 50th percentile. This similarity is highly noticable in the corresponding inequality 2.16, evaluated the same way as the 2.15 in order to extract a certain percentile intensity xpercentile.

xpercentile−1 X i=0 H(i) < n · percentile 100 ≤ xpercentile X j=0 H(j) (2.16)

When using percentiles however, it’s more common to evaluate the percentile from a known intensity bin l, i.e. the opposite way:

percentile = 100 n l X i=0 H(xi) Standard Deviation

The term deviation is a measure of the difference between an observed value and the mean value in a data set. In order to measure statistical dispersion however, the term variance is more videly used, measuring squared distances to the mean instead, averaging to obtain the variance of a data set. While the mean value of a data set is a way of expressing the location of a distribution, the variance is used

2.5 Tree Structures 21

to express its spread. One reason to deal with squared deviations is to always use positive numbers. This way the different components can’t cancel each other out to eliminate display of a present spread. The standard deviation is defined as the square root of the variance, and is a more commonly used term, since its results are in the same units as the data, therefore easier to interpret.

Eq. 2.17 shows the mathematical expression of standard deviation, with n being the number of elements and xi the observed intensity.

SD = v u u t 1 n n X i=0 (xi− xmean)2 (2.17)

As an example, look at the multi-set {3, 9}. The mean value is 6, resulting in the variance 1₂ (−3)2+ 32
= 9, and then standard deviation SD =√9 = 3. The corresponding equation from a histogram H of the data set is shown in 2.18.

SD = v u u t 1 n k X i=0 H(xi) · (xi− xmean)2 (2.18)

2.5

Tree Structures

Another technique used to technical depth in this thesis is tree structures. A tree structure is a data structure emulating a tree using linked nodes, connected to each other in a similar way to a real tree, with a stem separating into twigs and twigs separating into leaves. The first node, corresponding to the stem is called the root node, and it doesn’t have any superior nodes, also called parents in the tree. Nodes having a parent node are child nodes to the parent, and nodes not having any child nodes are called leaf nodes, since they are at the end of the structure, just like a leaf, see figure 2.14.

There are numerous different tree structures and the applications where tree structures come in handy are broad. In the scope of this thesis though, they are to be used for splitting 2D data sets for classification. Then the structure of a quad tree is suitable.

2.5.1

Quad Trees

A quad tree is a tree structure in which parent nodes each have four children, mak-ing it perfect for the task of recursively subdividmak-ing a 2D data set into quadrants, see figure 2.15. A region, such as the pixels in an image or positions in space is set to be the root node. According to an evaluation of the regions, such as estimating tissue likelyness in an image region, the algorithm either decides to classify the region and stop the splitting, or to split it and evaluate each child region in order to get an easier decision.

Figure 2.14. Basic binary tree structure.

2.5 Tree Structures 23

The applications are wide even within the scope of using quad trees. Examples of areas where they are used are image representation, collision detection in two dimensions and view frustum culling in computer graphics. In three dimensions the corresponding structure to a quad tree is the octree, with each node having 8 child nodes, and the additional 4 corresponding to a split in the additional dimension.

Chapter 3

MR Grayscale Mapping

Challenges

Through the process of translating spin physics into frequencies and further into Hydrogen densities at positions in the body, see 2.2, we obtain data sets for visu-alization through mapping procedures explained in 2.3. This however, is just the pipeline describing the workflow. In both theory and practice, there are lots of uncertainties and parameters affecting the outcome, making grayscale mapping of MR data sets a very complex task.

This chapter brings diversities and problematic issues causing the difficulties of MR grayscale mapping to the subject. It also presents rescent research and specialized techniques used in this thesis, orginating from [6] and theory presented in chapter 2.

3.1

MR Uncertainties

By taking a closer look at the parameters affecting display output of a 2D MR visualization, we immediately see that many of them have similar effects. Different pulse sequences cause the received spin density signal to vary. The same signal is amplified by the strength of the magnetic field and the actual anatomical spin density, which is related to the current water level in the anatomical region. On top of this, weightings are applied affecting the appearance of different tissues and intensity mapping finally decides display output for the resulting data set, causing all previous dependencies to be summed up into one variable - the received signal. To take the challenge even further, the same display output will appear differently on different monitors and in different lighting environments [14], inducing a further need of analysis in the visualization pipeline.

This thesis work takes place late in the process, i.e. at the window level and interaction states, and therefore the challenges orginates from the signal output, i.e. factors causing the received signal to vary. These can be summed up into uncertainties of weighting and the lack of an absolute scale.

Figure 3.1. Example of poor result with a min/max window level.

3.1.1

Lack of Absolute Scale

The fundamental issue with MR is that there is no calibrated value range, and in order to achieve accurate diagnostic review, the received intensities have to be interpreted in the context of the protocol used to capture them. But due to anatomical dependencies mentioned above, this won’t solve the problem, since two different patients or scans using the same protocol and presets may create different data sets.

This is not the case with CT for example, where the values always describe X-ray attenuation calibrated in Hounsfield units, (HU) [6]. Air corresponds to -1000 HU and water to 0 HU. Given tissue types correspond to a fairly constant HU value, making visualization of certain tissue trivial.

A common solution to the lack of absolute scale is to visualize MR data sets with a min/max window, meaning to adjust the window to display the highest and lowest represented intensities as 255 and 0 respectively. This is inadequate in many situations though, first of all because the maximum and minimum intensity values are very instable quantities [14], since they are storngly influenced by spikes and background noise. Second, the value range of interest in the intensity domain may occupy only a portion of the dynamic range [14], and a min/max window will often represent the tissue of interest flat and without contrast, see figure 3.1.

3.2 Tissue Detection 27

Figure 3.2. Range weight combines data values and spatial neighborhoods. Here, the

range weight for range Φ = [3, 4] are shown for each 2x2 neighborhood.

3.1.2

Weighting

Another challenge in MR grayscale mapping is the differences in data sets due to different weightings. As mentioned before, see 2.2.4, different pulse sequences, T R and T E values will emphasize contribution from tissue T1, T2or pure spin density

ρ. These may need to be viewed, and hence visualized differently, because posi-tioning of interesting tissue range and appearance of abnormalties in the data set may differ. These facts indicate that a uniform window level algorithm is imprac-tical, and induce a need for a window level algorithm to adapt to the weighting [14].

3.2

Tissue Detection

In order to address the challenges presented in 3.1, MR data sets have to be further analyzed. Lots of research efforts have been put into this subject, both for 2D and 3D cases, and the overall starting point is to look at the data set histogram, see 2.4.1. Histograms provide good information about the data set to a low computational cost of extraction. The idea behind using histograms is that relevant tissues will stand out as peaks [6]. Due to the nature of MR scans though, this is seldom the case. Usually the major parts of the data set correponds to air and noise, clouding the prescence of interesting tissues. The amount of these non-interesting parts differ a lot from examination to examination, inducing a need to detect which parts of the data set that actually corresponds to tissue.

The challenge of automatic adaption of visualization parameters for MR have been widely addressed before, [17] [9] as examples. These are both based on the shape of the global histogram, with the limit of only being applicable for examinations where the shape of the global histogram is consistent from patient to patient [6].

This thesis work employs two methods for extraction of data corresponding to intersting tissue, 4.2.2 and 4.2.3, both orginating from techniques proposed by Lundström et al. [6], and explained in detail in 3.3.1 and 3.3.2.

3.2.1

Range Weights

A range weight is a statistical property introduced by Lundström et al [6], and it is a measure of a neighborhood footprint in a given intensity range. In MR and medical imaging in general, features to be studied and visualized are often spatially concentrated. To this end, the range weight is used in order to incorporate spatial intensity relations in the data analysis. For a given partial intensity range, the range weight measures how well a local neighborhood is described by that range, see figure 3.2. For a neighborhood N and set of elements NΦ⊆ N in the intensity

range Φ, the definition of a range weight wr is given in eq. 3.1, |A| being the

cardinality of set A.

wr(Φ, N ) =

|N ∩ NΦ|

|N | (3.1)

3.2.2

Neural Networks

As mentioned earlier, founding automatic adaption of visualization parameters purely on the global histogram is often impractical, see 2.2.4. Data sets with simi-lar histograms may have very different spatial distributions, hence it is desirable to have some adaptation capability biult-in. A tool suitable to accomplish this, also used in rescent research is an aritificial neural network, (ANN), providing “learn by example” capability.

An ANN is a simulated artificial group of connected neurons, which uses a mathematical model for information processing [13]. There are three ways of learning in ANN - supervised, unsupervised and reinforced learning. If used in MR, an ANN could learn by studying how a radiologist manually adjusts the window level of certain data sets, and then find patterns to later use for automatic adjustments.

There are several methods proposed in the area of automatic window level ad-justment utilizing neural networks. In [15], features are extracted from the global histogram, and a neural network is used to learn the relationship between these features and desired window level output of a radiologist. Inventions by [14] gen-erats data set features from a combination of the global histogram and computed spatial information, generated from statistical properties of data set elements clas-sified as tissue. However, while these methods both solve the versatility of an MR scan, they also share the drawback of being incompatible with MR images they havn’t been trained for [14]. An example of this would be if a new scanning pulse sequence, see 2.2.4, was utilized. In such a case, the software would have to be re-programmed and re-trained, which would be impractical for a hospital.

3.2.3

Wavelet Features

Wavelets are mathematical functions used to split time signals into different fre-quency components. Similar to the fourier transform, the wavelet transform is the representation of a function by these wavelets, exploting frequency information.

3.3 Spatial Coherence 29

Figure 3.3. Grid showing wavelet coefficient positions for corresponding filtering.

For the purpose of imaging, a 2D discrete wavelet transform (DWT) is used, having two major advantages over the fourier tranform. First of all by running in linear time n instead of nlog(n) as the Fourier Transform, and second since it also preserves spatial information, thus representing both spatial and frequency content at the same time [4]. During transformation, the wavelets work as high and low pass filters, splitting and downsampling the frequency content into wavelet coefficients, which are stored according to the combination of filtering applied, see figure 3.3. Further, the transform is usually applied recursively on the LL quadrant, decomposing the low frequency content even more. A decomposition and resulting storage of wavelet coefficients is demonstrated in figure 3.4.

Wavelet features are properties of these wavelet coefficients [4], similar to how mean value and standard deviations are properties of a spatial data set. There are many ways of defining wavelet features, but in order to use them for classification, a neural network is utilized and trained on a wide variety of cases. Wavelet fea-tures have been successfully utilized in [14] for example, however in the scope of generating a window level for visualization, their use is limited. The classification of features looks for non-usual patterns, hence it finds abnormalties rather than interesting tissue in general.

3.3

Spatial Coherence

Histograms of data sets are great complements when analyzing the contents of an MR examination, providing an intuitive summarization of the distribution. The main drawback though, is that a histogram is only dependent on intensities, totally ignoring shape and texture, i.e. all spatial relations are lost [6]. Two histograms can potentially be identical for two MR data sets with different content. Contrary,

Figure 3.4. A level 2 wavelet transform of an image. Low (L) or high (H) pass filtered

content are stored as demonstrated in figure 3.3

Figure 3.5. Two images of the same content, but very different intensity distributions,

3.3 Spatial Coherence 31

Figure 3.6. The effect of the alpha histogram.

two MR scans of identical content can be indistinguishable using histograms, due to different weightings, see figure 3.5.

As mentioned before, the fundamental assumption motivating histogram use in the context of automatic window level adaption is that interesting tissue will show up as peaks. This is seldom the case, but such can be obtained through tissue classification. This thesis proposes algorithms utilizing 3 methods in order to accomplish classification of tissue, 2 of them proposed by Lundström et al. [8] [7]. Introducing spatial information, they emphasize coherent data, usually corresponding to tissue of interest.

3.3.1

Alpha Histograms

Through the loss of spatial relations, the global histogram doesn’t have any in-formation about where in the data set the different intensities come from. By studying the data set in smaller fractions however, this problem can partially be overcome. The alpha histogram is a method exploiting this, by splitting the data set in fractions and retrieve local histograms of these regions to operate on. Within each smaller fraction, spatial relations are of course still lost in the local histogram, but since a location restriction in the data set is introduced, each local histogram exploits the global spatial relations to a certain extent, depending on fraction size. After operation on local histograms, these are summed up and normalized into a new global histogram, in which spatially coherent regions are amplified, see figure 3.6.

The operation performed on local histogram is a raise of each bin count to the power of α. Coherent tissue will cause peaks in these local histograms, and through a raise of α > 1, the amplification of different bin counts will be disproportional, and coherent tissue will emerge [8]. These local histograms are then summed up and the new global histogram is raised to the power of _α1. This operation is disproportional to the amplification of the local raises of α, hence an additional normalization is done, to make the new histogram integral equal the initial, i.e. the summed bin counts should equal the sum of elements in the data set. This last operation is a proportional amplification though, and the shape of the histogram is sustained [8].

functions Hn(N, x), counting the instances of intensity x in a local neighborhood,

see eq. 3.2 [8], which is a modified version of eq. 2.12. N is the set of data values in the neighbourhood, Dx is the set of data values equalling x in the data set D,

and |S| is the cardinality of set S.

Hn(N, x) = |N ∩ Dx| (3.2)

The neighborhood subdivision of the data set creates spatial regions N1, ..., Nk.

The alpha histogram Hα(x) is the sum of all the enhanced local histogram

func-tions, raised to the power of _α1, see eq. 3.3. An important case for understanding purposes is what would happen when α → ∞. The resulting bin counts of the al-pha histogram would equal each bin’s maximum count among the local histograms, since a raise to the power of alpha would cause all lower represented values to van-ish in comparison to the maximum, i.e. the alpha histogram for α → ∞ equals a maximum operator among the local histograms.

Hα(x) = k X i=1 Hn(Ni, x)α !α1 (3.3)

Applications of the alpha histogram and modified versions of it are many, though the original purpose was to ease transfer function design for DVR [8]. In such case, the new histogram would simplify the task of separating overlapping tissue, as well as to distinguish peaks from coherent tissue in the histogram. In the scope of this thesis, the property of amplifying coherent tissue can be used to determine a suitable range for the window level, see chapter 4 and 5. Results from amplifying coherent material are demonstrated in figure 3.6.

3.3.2

Partial Range Histograms

The Partial Range Histograms (PRH), is also a method utilizing a subdivision of the data set. The definition of a PRH is the histogram for a set of neighbor-hoods that are typical for a given intensity range [7]. For each local histogram, a range weight, see 3.1, is calculated, and more specifically, a PRH groups all blocks exceeding a threshold in terms of these range weights.

Mathematically, the local histogram functions are calculated according to eq. 3.2, and the PRH Hprh(Φ, x), Φ being the intensity range, is the histogram of

local regions Ni having a range weight wr(Φ, Ni) exceeding a threshold τ , see eq.

3.4. Hprh(Φ, x) = k X i=1 |(Ni∩ (wr(Φ, Ni) ≥ τ )) ∩ x| (3.4)

A problematic task is to find a relevant value range Φ. PRH was originally used in DVR for tissue detection [7], in which a manual selection is possible during

3.3 Spatial Coherence 33

Figure 3.7. Partial Range Histogram (PRH) example: Left, subdivision into block

neighborhoods. Center, pixels within a given intensity range are highlighted. Right, the PRH is the histograms for the group of neighborhoods exceeding a range weight threshold. The 80% threshold is shown in yellow and additional blocks for the 50% theshold in orange.

interaction [6]. That would not be possible for use in 2D slice viewing, due to its time consumption. An automatic tissue detection scheme based on the PRH has also been developed [7], but it doesn’t exactly translate to the context of this thesis, since all interesting tissue are of interest instead of particular tissues. This induces a need for a statistical analysis in order to determine Φ, explained in depth in 4.2.3. A demonstration of how a PRH works and can be used is shown in figure 3.7.

Chapter 4

Equalized MR Grayscale

Mapping

This chapter presents the proposed algorithms used to solve the problem descrip-tion, see 1.1. These consist of theories presented in chapter 2 and 3 put together into techniques, which were implemented in IDS7 during thesis work.

As explained earlier, when visualizing an MR data set, it is a huge and often completely false assumption to think that all the data are of interest. There are many applications of MRI, but it is significantly used to detect abnormalities in soft tissue. Fat, and other highly hydrogen concentrated anatomy are seldom of interest, and noisy air never corresponds to intensities in need of contrast from the window level, and so it would be of great interest if we could identify those areas in the data set, thus not taking them into account when designing parameters for visualization.

Due to the lack of an absolute scale, different purposes of MR scans and also the different types of representations, see 3.1, it is not possible to find a general case how to handle MR data sets. Also, the amount of air, or noise, as it most often appears in data sets, differs a lot from scan to scan making clever conclusions about region importance a very tough task without a data set subdivision into smaller regions, i.e. to exploit spatial relations in the analysis, as explained in 3.3. Basic calculations and operations such as equalities to histogram equalization, i.e. min/max windowing can be done without splitting, but these are very simple methods without much thought about correlated, smooth or interesting tissue regions. Depending on the data set they may also lack contrast due to large spans of intensities, and may show a great deal of noise if the dynamic range of the data set is short. On the other hand they are very fast and therefore useful for comparison, both for resulting window level and algorithm speed, see 5.

The objectives of this thesis, deinfed in 1.2, aim at finding good window level presets and intuitive, data set related manual window level interactions through the PACS. To structure the algorithmic solutions proposed to solve these problems, there are four ground areas in need of investigation - subdivision, evaluation, generation and interaction, each brought to the subject in this chapter.