WO2017183958A1 - Method for quantifying lung kinematics - Google Patents

Abstract

Computer-implemented method for quantifying lung kinematics. A sequence of three dimensional lung surface meshes (Lxyz) representing a recorded lung motion is calculated. A first lung surface mesh (L1) is selected at a first recording time (t1). Part of the selected first lung surface mesh (L1) is identified as a first subsurface mesh (S1). A second lung surface mesh (Lend) is selected at a second recording time (tend). Part of the second lung surface mesh (Lend) is selected as a second subsurface mesh (Send). The identification of the second subsurface mesh (Send) corresponds to the identification of the first subsurface mesh (S1) propagated through time (t1…tend) during the recorded lung motion. A kinematic parameter (A) is calculated based on a comparison of the respective corresponding subsurface meshes (S1,Send) at different recording times (t1,tend).

Description

Title: METHOD FOR QUANTIFYING LUNG KINEMATICS

TECHNICAL FIELD AND BACKGROUND

The present disclosure relates to a method, system, and computer readable medium for quantifying lung kinematics.

Respiratory failure, characterized by impaired oxygenation (oxygenation failure) and/or carbon dioxide elimination (ventilatory failure), can have multiple reasons. One specific cause is dysfunction of the respiratory musculature. Different methods and modalities to measure and assess weakness or dysfunction of thoracic and diaphragmatic musculature and respiratory dysfunction in general are described herein. For example, Pompe disease (glycogen storage disease, type 2) is an inherited

neuromuscular disorder characterized by progressive limb-girdle weakness and pulmonary insufficiency. Enzyme replacement therapy (ERT) may positively alter the course of Pompe disease and stabilize or improve motor function. However, the therapy may be mostly effective when the disease is detected in its early stages. Accordingly, quantitative methods for

measuring diaphragm impairment are desired to determine the optimal start of ERT and to monitor and predict treatment response. Other examples for which quantitative methods are desired may include patients suffering from chronic obstructive pulmonary disease (COPD) and children born with a defect of the diaphragm.

Pulmonary function tests (PFT) can be used to evaluate general respiration performance by breathing into a spirometer. However, in lung patients, the contribution of the diaphragm in breathing cannot be evaluated using PFT due to compensatory efforts of the intercostal, accessory and abdominal musculature. Other useful modalities to assess the structure and function of respiratory muscles include high resolution muscle magnetic resonance imaging (MRI) and computed tomography (CT). In particular, dynamic imaging may enable direct functional analysis of the respiratory system. In dynamic imaging, multiple images (referred to as frames) acquired at different time points are combined to an image sequence. This can be achieved in 2D (e.g. at mid-sagittal plane in right and left lung or mid-coronal plane), or in 3D. The characterization of the respiratory system based on dynamic images may involve e.g. the extraction of geometry or motion features. However, manual characterization of relevant features is time-consuming especially in a sequence with many three-dimensional frames.

There is a need for automatic characterization of lung features in dynamic MRI sequences of three dimensional images to determine relevant kinematic parameters.

SUMMARY

Geometry features of the lung can in principle be obtained from a single frame, independent of the other frames in the sequence. For example, automated lung segmentation in static 3D MRI data may involve, 3D region growing, confidence-connectedness segmentation, or statistical shape models in combination with deformable mesh segmentation. Extracted geometry features may include e.g. lung volume, lung size (in anterior- posterior, cranial-caudal and left-right direction), diaphragm length (2D images) or diaphragm surface area (3D images), both divided into zone of apposition and diaphragm dome. However, the structural analysis may not capture the actual functional effects on the respiratory system.

Some aspects of the present disclosure therefore provide a method for quantifying lung kinematics, i.e. movement of the lungs. In one embodiment, the method comprises receiving a sequence of three- dimensional images of a lung motion. The sequence may be recorded as a function of time during a breathing maneuver. For example, the three- dimensional images may be obtained using magnetic resonance imaging (MRI). Other or further aspects of the present method may comprise calculating, based on a sequence of three-dimensional images, a sequence of three dimensional lung surface meshes representing the recorded lung motion. Motion features can be extracted in different ways, e.g.

segmentation- and registration-based. For example, deformable mesh segmentation can be applied successively to the 3D frames of an image sequence while incorporating temporal coherence. Yet, the robustness of this approach with respect to image artifacts (e.g. ghosting at lung tissue borders), spatial/temporal accuracy, and smoothness of the segmentation remains challenging.

Preferably, the sequence of meshes is provided by a single three dimensional mesh and a transformation which models deformation of the mesh at different times. For example, the deforming mesh may comprise the same number of points and the same layout of connections between the points, even though coordinates of the points may change as a function of time. Alternatively, or in addition, the sequence of meshes may comprise different meshes, e.g. different numbers of points and/or different

connections for different times. Correspondence between subsurface meshes may be established e.g. on criteria of the individual meshes or based on a more general modelling of transformation between the different meshes.

Aspects of the present method may involve selecting a first lung surface mesh in the sequence of three dimensional lung surface meshes corresponding to a first recording time. For example, the first recording time may be at the start of the breathing maneuver, e.g. at end-inspiration when lung volume is maximal. Part of the selected first lung surface mesh may then be identified (categorized) as a first subsurface mesh. For example, the first subsurface mesh may correspond to lower part of the lung associated with diaphragm induced breathing motion.

Further aspects may involve selecting a second lung surface mesh in the sequence of three dimensional lung surface meshes corresponding to a second recording time. For example, the first recording time may be at the end of the breathing maneuver, e.g. at end-expiration when lung volume is minimal. Part of the second lung surface mesh may then be identified as a second subsurface mesh. The identification of the second subsurface mesh corresponds to the identification of the first subsurface mesh propagated through time during the recorded lung motion. In other words, a

representation of the same part of the lung surface is tracked through time. For example, if the first subsurface mesh corresponds to a diaphragm surface at the first recording time, so the second subsurface mesh will correspond also to the same diaphragm surface but at the second recording time. Typically, the second subsurface mesh will be displaced and/or deformed with respect to the first subsurface mesh as a result of the breathing movement.

By propagating the identification of a specific subsurface through time, the identification need only be done once and automatically yields a corresponding subsurface at a different time. By comparing the respective corresponding subsurfaces at different times, additional information can be retrieved compared e.g. to general evaluation of respiration performance by breathing into a spirometer. For example, the functioning of the diaphragm can be specifically monitored and compared to the rest of the lung function.

Comparison of the respective corresponding subsurface meshes at different recording times may result in the calculation of a kinematic parameter. The kinematic parameter may be used to quantify the lung movement. For example, performance of a diaphragm can be quantified for analysis. In some embodiments, the quantification may be visualized, e.g. as a number or graph. For example, diaphragm displacement distance and/or displaced volume can be shown. The kinematic parameter can also be based on a combination of lung kinematics, e.g. comparing the displaced lung volume of the diaphragm to the volume displaced by other lung surfaces, e.g. costal surface and/or the total displaced volume. Comparing lung surfaces at the points of end-inspiration and end expiration, may yield the most prominent result as the difference is typically maximal. These points may also be more reproducible and comparable, for the same patient and between patients. For example, the kinematic parameter may be calculated based on a diaphragm contribution to displaced volume between the diaphragm surface at end-inspiration and the diaphragm surface at end-expiration. Calculation of displaced volume of a subsurface may e.g. comprise closing the volume by connecting edges of the same subsurface at different times during the breathing maneuver.

It will be appreciated that also other or further subsurfaces can be used. For example, a third subsurface may be identified in the sequence and the calculation can be based on three subsurfaces. Also four, five, or more subsurfaces can be used in the calculation. In particular, the whole sequence of subsurfaces can be used as basis for the calculation. For example, some parameters may be based on an average or extreme (minimum or

maximum) value in the subsurface at any point on time during the sequence. For example, a maximum angle of the diaphragm surface during the breathing motion may be calculated.

Also other relevant parameters may be defined from the diaphragm surface, e.g. a change in diaphragm orientation. Relevant parameters can also be derived from other subsurfaces or combinations of subsurfaces. For comparison between patient groups, it may be useful to normalize the results by total lung volume, which can be either measured in the image or estimated from patient height and gender (predicted total lung capacity). Advantageously, the total volume may be determined by using the whole lung surface. The calculation of the lung volume at different times can also be used to determine the points of end-inspiration or end- expiration, e.g. when the total volume is maximal or minimal, respectively.

Image segmentation methods can be used to partition the image into meaningful segments. The segment boundaries (e.g. lung surface or diaphragm) may be obtained by manually tracing object silhouettes, semi- automatically including little user interaction, or fully automatically.

Manual segmentation of the diaphragm and lungs is a time-consuming and tedious task, especially in 3D images or images sequences. On the other hand, automatic extraction of the diaphragm segment has hitherto remained challenging.

It is presently recognized that, from an anatomical point of view, the lung surface can be divided into three regions and two boundaries. The inferior region of the lung adjacent to the diaphragm is referred herein as the base or diaphragmatic surface. It is enclosed by the inferior border and may generally be characterized as a completely concave surface. The large convex region close to the ribs is referred to as the costal surface and consists of the anterior, posterior and lateral sides of the lung surface. The apex is defined as the superior part of the lung which extends over the first rib and can also be considered as part of the costal surface. The residual region of the lung surface lies next to the mediastinum and, hence, is called medial or mediastinal surface. It is bounded by the anterior border and the posterior border.

A lung surface mesh is typically built from nodes and edges which e.g. form interconnecting triangles in three dimensional space. To improve robustness of the subsurface classification algorithm, parts of the mesh surface may be first grouped into patches. As used herein, the terms subsurface and patch are used which in some cases may both describe a connected subset of nodes and triangles of the input mesh. The subsurface may generally refer to larger meaningful regions and patches may be small collections of triangles and nodes which may e.g. occur in the context of watershed transformations.

To make identification of the diaphragmatic surface easier, it is found preferable, in some embodiments, to first make a general distinction based on the concavity/convexity of the surface. For example, a maximum set of interconnected convex patches (or larger subsurfaces) may be identified as the costal surface. The non-identified residual surface may then be further subdivided into the diaphragm surface and the mediastinal surface. The further subdivision of the residual surface may be based on other anatomical observations or criteria. For example, it is found that the most superior patch of the residual surface typically belongs to the mediastinal surface and the most posterior patch on the residual surface typically belongs to the diaphragm surface. Also other criteria may be employed, e.g. that patches in the diaphragm surface are entirely visible when viewing the residual surface from the inferior side, e.g. from a center point of the generally dome shaped diaphragm surface. Terms such as superior, posterior, and inferior are medical terms which may refer to the frame of a patient body from which the images were acquired, e.g. top-side, back-side and bottom-side, respectively. This orientation may be recorded at the same time as acquiring the images, inferred from a standard position of the patient with respect to the recording device (e.g. MRI), or inferred from the images themselves (e.g. shape, orientation, and relative position of recorded organs)

The initially identified patches can be used as a starting point to find further patches. Adjacent patches may be merged based on different criteria, e.g. the similarity in (average) orientation, a distance between vertices at an edge between the neighboring patches, and a curvature of the neighboring patches being concave or convex. By selecting initial seed points for the random walker at a triangle in the center of each patch, the chance of wrong classification may be lower. For example, a random walker algorithm can be used to fill in the unidentified pieces. The random walker algorithm can be used e.g. to generate a probability map indicating a probability that a part of the residual surface belongs to the diaphragm surface and mediastinal surface. The random walker algorithm can be iterated to improve confidence of the classification. For example seed points can be added based on a previous iteration. The iterations may be

terminated when there are no more significant changes.

Once a specific subsurface is identified for one recording time, it can be advantageous to propagate the identification to other recording times. To achieve this, it is useful to establish a correspondence between subsurfaces at different times. By modelling the lung motion, the movement of the surface can be mapped as a function of time. For example, a deformation field or other transformation can be calculated based on the three dimensional (voxel based) images. The deformation field may then be applied to one of the surface meshes to calculate surface meshes at other times. In this way also the identification of a subsurface can be propagated to corresponding subsurfaces.

The deformation field can be calculated e.g. based on a registration algorithm which models the lung motion according to a non- rigid transformation model to take into account deformability of different organs. The registration algorithm may be based on an optimization problem where a cost function is minimized with respect to parameters that describe a transformation between three-dimensional images at different times in the sequence. The transformation may also be calculated between images of the sequence and a reference image, e.g. average image. To calculate transformation between images, the transformation to and back from the reference image may then be concatenated.

The cost function can be based, e.g. on a variance of intensity values of voxels in the three-dimensional images over time. For example, the registration algorithm may utilize a four-dimensional free-form, e.g. B- spline, deformation model to fit the data. To guide the registration algorithm, a mask can be used to select the most relevant parts of the images. To further improve performance of the registration algorithm it is sometimes best to preprocess the images, e.g. correct any intensity inhomogeneities, where necessary. This is also referred to as bias field correction. Inhomogeneities may e.g. be caused by non-ideal behavior of the measuring equipment (e.g. MRI) and/or transmission through tissues.

Some aspects of the present disclosure may be embodied as software instruction on a computer-readable medium. Some aspects of the present disclosure may also be embodied as a system for quantifying lung kinematics. For example, the system comprises a computer-readable medium storing instructions that causes the system to execute methods as described herein. The system may also provide other or further

functionality. For example, the system may comprise (or be coupled to) a display for displaying a representation of the kinematic parameter or other information resulting from execution of the software instructions. For example, the system may comprise (or be coupled to) means for recording the sequence of three-dimensional images such as an MRI apparatus. BRIEF DESCRIPTION OF DRAWINGS

These and other features, aspects, and advantages of the apparatus, systems and methods of the present disclosure will become better understood from the following description, appended claims, and accompanying drawing wherein:

FIG 1 shows an overview of an embodiment for the image analysis pipeline;

FIG 2 shows an example embodiment for mesh representation of 3D lung segmentation;

FIG 3 shows an overview of an embodiment for a lung surface partitioning procedure;

FIG 4 shows example results of deviation between manual and automatic lung size computation;

FIG 5 shows example results of chest wall and diaphragm contribution to overall lung volume change; FIG 6 shows example results of diaphragm orientation and volume displacement during exhalation.

FIG 7 shows example results of cranial excursion of the right hemidiaphragm: anterior vs. posterior;

FIG 8 show example results of maximum excursion of points on costal and diaphragm surface.

FIGs 9A-9C illustrate steps of an embodiment for quantifying lung kinematics;

FIG 10 illustrates a flow diagram of an example procedure to distinguish the costal surface (convex) from the rest of the lung surface (concave), namely the base and mediastinal surface;

FIG 11 illustrates a flow diagram of an example procedure for automatically segmenting the lung base and mediastinal surface in the concave parts of the lung surface mesh;

FIG 12 shows a table with detailed results for the proposed image-based features and the spirometry features.

DETAILED DESCRIPTION

In some instances, detailed descriptions of well-known devices methods including algorithms may be omitted so as not to obscure the description of the present systems and methods. Terminology used for describing particular embodiments is not intended to be limiting of the invention. As used herein, the singular forms "a", "an" and "the" are intended to include the plural forms as well, unless the context clearly indicates otherwise. The term "and/or" includes any and all combinations of one or more of the associated listed items. It will be understood that the terms "comprises" and/or "comprising" specify the presence of stated features but do not preclude the presence or addition of one or more other features. It will be further understood that when a particular step of a method is referred to as subsequent to another step, it can directly follow said other step or one or more intermediate steps may be carried out before carrying out the particular step, unless specified otherwise. All publications, patent applications, patents, and other references mentioned herein are incorporated by reference in their entirety. In case of conflict, the present specification, including definitions, will control.

Some aspects of the present disclosure are supported by a patient study as presented herein. While, the study makes reference to particular patients, instruments, and algorithms, the skilled artisan will appreciate that the disclosure is not limited to specific details. In the study we investigate e.g. the suitabihty of dynamic MRI imaging in combination with automatic image analysis methods to assess respiratory muscle weakness. For example, diaphragm weakness is found to be a main reason for respiratory dysfunction in patients with Pompe disease, a progressive metabolic myopathy affecting respiratory and limb-girdle muscles. Since respiratory failure is the major cause of death among adult patients, early identification of respiratory muscle involvement is necessary to initiate treatment in time and possibly prevent irreversible damage.

The diaphragm can be considered as a flexible membrane that is attached to the ribs and separates the thoracic from the abdominal cavity. As a consequence of contraction and relaxation of the diaphragm during breathing, a part of this big muscle rests against the rib cage (referred to as zone of apposition) and the other part rests against the lungs and heart (referred to as diaphragm dome). Normally, the diaphragm is the primary breathing muscle and diaphragmatic movement normally accounts for approximately 75% of lung volume change during breathing. The remaining 25% is then performed by movement of the rib cage, via the intercostal muscles. In patients with Pompe disease the diaphragm movement may be impaired which can manifest e.g. as reduced diaphragm volume

displacement and/or different orientation of the diaphragm which can be quantified by the present methods. Some embodiments, described herein, may include image registration and lung surface extraction to automatically quantify lung kinematics during breathing. This allows for the extraction of geometry and motion features of the lung e.g. characterizing the independent contribution of the diaphragm and the thoracic muscles to the respiratory cycle.

Advantageously, dynamic 3D MRI may provide data for analyzing the contribution of both diaphragm and thoracic muscles independently. Some aspects of the present image analysis thus have the potential to detect less severe diaphragm weakness and can be used to determine the optimal start of treatment in adult patients with Pompe disease in prospect of increased treatment response.

The present study includes 16 adults, of which 10 were diagnosed with late-onset Pompe disease (referred to as P01-P 10) and 6 age- and gender-matched controls (referred to as C01-C06). At data acquisition, the median age of the patients was 46 years (range: 32-66) and for controls 43 years (range: 27-55). Both groups consisted of 50% females. Three patients required nocturnal ventilation and one was wheelchair dependent. The median time since diagnosis was 16 (range: 9-30) years and all patients received enzyme replacement therapy (range: 0-7 years, median duration: 5.5 years). While it is appreciated that the present methods can provide particular advantage to the identification and analysis of Pompe disease, the disclosure may also be applied in other circumstances for quantifying lung kinematics.

In the study, subjects were scanned on a 3T GE Signa 750 MRI scanner (General Electric Healthcare, Milwaukee, USA) using a 32 -channel torso coil. Static, full-inspiration and full-expiration scans may be acquired for example during 12-seconds breath-hold using a 3D RF-spoiled gradient echo sequence (3 mm slice thickness, 1.5 mm slice spacing, 1.4063 x 1.4063 mm in -plane resolution). While it is appreciated that the present methods can provide particular advantage to the analysis of 3D MRI images, the disclosure may also be applied to data acquired from other sources. More information on a possible implementation of the acquisition procedure, as used herein can be found for example in Wens et al. ("Lung MRI and impairment of diaphragmatic function in Pompe disease". BMC Pulmonary Medicine. 2015; 15(1): 1-7).

In the study, results in 16 3D+t MRI scans (10 Pompe patients and 6 controls) of a slow expiratory maneuver show that kinematic analysis from dynamic 3D images reveals important additional information about diaphragm mechanics and respiratory muscle involvement when compared to conventional pulmonary function tests. Pompe patients with severely reduced pulmonary function showed severe diaphragm weakness presented by minimal motion of the diaphragm. In patients with moderately reduced pulmonary function, cranial displacement of posterior diaphragm parts was reduced and the diaphragm dome was oriented more horizontally at full inspiration compared to healthy controls.

For comparison, the present study also features a dynamic sequence which was acquired while the subjects performed a slow

exhalation maneuver from total lung capacity (TLC) to residual volume (RV). This sequence contains 48 3D volumes which were continuously acquired at a rate of approximately 2.5 volumes per second. Fast acquisition was achieved using the time resolved imaging of contrast kinetics (TRICKS) technique and a reduced spatial resolution (12 mm slice thickness, 6 mm slice spacing, 1.875 x 1.875 mm in-plane resolution). TRICKS is described for example by Korosec et al ("Time-resolved contrast-enhanced 3D MR angiography." - Magnetic Resonance in Medicine. 1996;36:345-351.) To ensure correct execution of the maneuver during static and dynamic imaging, both the preliminary training and the scanning were controlled with an MR-compatible spirometer and instructions were given by a certified lung function technician who also supervised the quality of the maneuver. In some embodiment, the dynamic sequences can be manually cropped, e.g. discarding the tidal breathing after full exhalation. In the present study, the sequence length T (expressed as mean±std) for all subjects thus decreased to 36±8.9 frames (patients: 33.6±9.8, volunteers: 40±5.3). In the study, a trained medical observer manually segmented right and left lung in both static scans (inspiration and expiration) of each subject by tracing the lung surface in every second axial slice using a 3D slicer. Slices without annotations were filled automatically via interpolation.

Pulmonary function tests (PFT) were also performed prior to the MRI scan to measure vital capacity (VC), forced vital capacity (FVC), forced expiratory volume in one second (FEV1), peak expiratory flow (PEF) and maximum inspiratory and expiratory pressure (MIP, MEP). The

measurements VC, FVC, FEV1 and PEF were acquired both in sitting and supine position and expressed as absolute values as well as percentage with respect to predicted values. The postural drop AFVC was measured as the relative difference between FVC in sitting and supine position. A drop of >25% is commonly considered an indicator of diaphragm weakness.

The invention is described more fully hereinafter with reference to the accompanying drawings, in which embodiments of the invention are shown. This invention may, however, be embodied in many different forms and should not be construed as limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the invention to those skilled in the art. The description of the exemplary embodiments is intended to be read in connection with the accompanying drawings, which are to be considered part of the entire written description. In the drawings, the absolute and relative sizes of systems, components, layers, and regions may be exaggerated for clarity. Embodiments may be described with reference to schematic and/or cross-section illustrations of possibly idealized embodiments and intermediate structures of the invention. In the

description and drawings, like numbers refer to like elements throughout. Relative terms as well as derivatives thereof should be construed to refer to the orientation as then described or as shown in the drawing under discussion. These relative terms are for convenience of description and do not require that the system be constructed or operated in a particular orientation unless stated otherwise.

FIG 1 shows an overview of an embodiment for the image analysis pipeline. The gray boxes present general steps to process the input images (dynamic sequence, static image and manual segmentation of the static image) and to derive motion and geometry features of the diaphragm and lung surface. The main framework can be divided into four steps: preprocessing, deformation field estimation, lung surface partitioning, and motion analysis. In the following, the steps are presented in more detail.

In order to automatically extract characteristics of respiratory muscle dynamics, the motion of the lungs and surrounding structures may be captured e.g. via non-rigid image registration of the dynamic image sequence resulting in a dense deformation field. Based on these deformation fields, the contribution of the different muscular groups involved in the respiratory maneuver can be derived. By tracking the bottom of the lung surface the diaphragm motion can be investigated in detail, while the chest wall contribution can be inferred e.g. by evaluating the movement of the anterior-posterior and left-right surfaces of the lung.

Step 1 - Image preprocessing / Bias field correction

In one embodiment, intensity inhomogeneity that are typical for

MRI scans is corrected. This is also referred to as "Bias Field Correction".

For example, the present study uses the algorithm N4ITK described by Tustison et al ("N4ITK: Improved N3 Bias Correction." - IEEE Transactions on Medical Imaging. 2010;29(6): 1310-1320.) Also other algorithms for correcting inhomogeneities can be used.

In one embodiment, an iterative approach is used to generate masks that define the image regions from which the bias field is estimated. The mask includes the body but excludes all air -filled regions, specifically the lungs and the image background.

In some embodiments, to generate such masks, the image is thresholded and the intensity bias field is computed using only the bright regions. For example, the image can be thresholded using Otsu's method described by Otsu ("A Threshold Selection Method from Gray-Level

Histograms." IEEE Transactions on Systems, Man, and Cybernetics.

1979;9(l):62-66.) In another or further embodiment, an image with reduced intensity inhomogeneity is then thresholded again to obtain a new, more accurate, mask. In the next iteration this new mask is used to correct the original image, e.g. with the N4ITK algorithm. In the present study, three iterations were used, as visual inspection showed that more than three iterations produced no further improvement, and the resulting images were of sufficient quality for the next image processing steps. During bias field correction of the static scans, the corresponding manually generated lung mask can be subtracted from the mask obtained via Otsu's method.

Step 2 - Motion estimation / 3D+t registration

The correct estimation of the lung deformation in the dynamic scans is important to determine the motion. For example, image registration is a method which estimates parameters of a transformation model in order to align two or more images. For this general optimization procedure, some typical components can be defined which determine, for example, how the similarity between the images is measured, what kind of transformation model underlies the data and how the optimization is performed. Typically dynamic MRI scans of the lung motion comprise multiple time points and exhibit a non-rigid deformation.

In one embodiment, the motion of the lungs in a dynamic sequence is estimated via 3D+t groupwise non-rigid registration (+t refers to the extra time dimension). The registration approach is explained for example by Metz et al ("Nonrigid registration of dynamic medical imaging data using nD+t B-splines and a groupwise optimization approach." - Medical Image Analysis. 2011 Apr; 15(2):238-249) using Elastix, described by Klein et al ("Elastix: A toolbox for intensity-based medical image registration." IEEE Transactions on Medical Imaging. 2010;29(1): 196-205). In this embodiment, the underlying deformation model uses a four- dimensional free-form B-spline that provides smoothness in the three spatial dimensions and in the time dimension. In one embodiment, 3D volumes of a sequence are aligned in an implicitly defined mean reference frame by minimizing the variance of intensity values over time. For example, registration can be achieved using an adaptive gradient descent optimizer and a multi-resolution scheme, e.g. with four resolution levels.

In order to restrict the registration to use image information in and around the lung, a lung mask can be created. In one embodiment, manual segmentations of the static full -inspiration and full-expiration scans, Min^s and M_ex ^s (with s referring to the static scans), are non-rigidly registered e.g. onto the first and last frames of the dynamic sequence respectively. For example, in the present study, mutual information and a 4- level coarse-to-fine pyramid approach with a final B-spline spacing of 50x50x50 mm was used, resulting in the masks Mi^d and Μτ (with d referring to the dynamic scans, 1 denoting the first and T the last frame). The union of the two binary masks was then dilated by 50 mm.

The above explained 3D+t non-rigid registration of the dynamic sequence estimates the B-spline coefficient vectors of the so-called forward (B-spline) transformation Τ_μ, which transforms each 3D volume to a mean reference frame. Additionally, an inverse transformation Τ_μ ^_1 can be computed to transform volumes of the mean reference frame to any time point in the original coordinate system. When both forward and inverse transformation are concatenated, a deformation field can be established which can be used to transform points or images of one specific time point, for instance the last frame (full-expiration), to any other time point.

Also other methods can be envisaged for motion field generation. For example, 4D segmentation instead of a 4D registration. For example, the first frame of the dynamic sequence is segmented using a shape based segmentation approach. The resulting 3D-surface model (e.g. triangle mesh) is then propagated to the next frame and adapted (e.g. via deformable models) to fit to the image data. The result can then again be used to segment the next frame and so forth. With this method points on the lung surface can directly be tracked throughout the sequence.

Step 3 - Mesh generation and partitioning of the lung surface. In one embodiment, a lung surface mesh, typically consisting of a set of vertices V and a set of triangles F, is obtained from the voxel-based image, e.g. a full-expiration lung segmentation, denoted herein as Μτ . In another or further embodiment, vertices of a 3D triangle mesh are propagated from one time point to other time points by applying a

deformation field which may be separately computed. This procedure may result in a 3D+t lung segmentation in which points on the lung surface are directly traceable throughout the breathing maneuver. For convenience of further processing, outward-pointing normals can be computed for each triangle.

Computer implemented methods can in principle be performed in any suitable program or device. For example, algorithms presented herein were implemented in MATLAB (R2013b, MathWorks, Inc.) using the CGAL 3.5 3D mesher (based on Delaunay refinement) which can be accessed through the mesh processing toolbox iso2mesh. See Fang et al ("Tetrahedral mesh generation from volumetric binary and grayscale images." - In: IEEE International Symposium on Biomedical Imaging; 2009. p. 1142-1145). Of course also other algorithms, operating on similar or other principle, can be used to obtain the lung surface mesh.

To assess lung function it may be desirable to determine the contribution of the different respiratory muscles during breathing. Thereto, chest wall and diaphragm motion may be investigated separately by analyzing the movement of lung surface regions adjacent to the respective anatomical structures.

Step 4 - Analysis and visualization.

When a complete identification of the subsurfaces in the 3D+t images is achieved, various parameters may be calculated for analysis of lung function. The data may also be visualized in various ways, e.g. as described in the following.

FIG 2 shows a mesh representation of 3D lung

segmentation. A 3D mesh of the lungs is presented in the top image. The surface is colored (shaded) to distinguish between the different surface segments. On the bottom, the same mesh is shown after unfolding the individual surface segments into a plane. The orientation labels indicate the viewing direction.

In one embodiment, the lung surface mesh, e.g. at full exhalation, is subdivided into three major subsurfaces associated with the

• chest wall: anterior (A), posterior (P), lateral (L,R), and superior (S),

• mediastinum: medial parts (M) adjacent to the heart), and

• diaphragm: inferior part (D).

In a preferred embodiment, this subdivision is achieved using the steps illustrated in Fig. 3 and detailed as follows. FIG 3 shows an overview of an embodiment for a lung surface partitioning procedure. The diagram shows steps to subdivide the surface mesh of a binary lung segmentation into six parts. In the two main parts the costal surface is extracted and the residual surface is divided into diaphragmatic and medial surface. An example with further details for extraction of the costal surface is shown in FIG 10. An example with further details for extraction of the diaphragmatic and medial surface is shown in FIG 11.

Extraction of the costal surface.

In one embodiment, the lung surface region adjacent to the chest wall, called costal surface, is extracted based on an anatomical observation that it is typically the most prominent large surface area which is convex. Both the diaphragmatic and medial lung surface regions are normally non- convex (concave) due to adjacent organs (heart and liver respectively).

In another or further embodiment, an algorithm is applied to the whole mesh in order to generate a subdivision into patches. For example, a watershed algorithm can be used wherein the patches are also referred to as watershed basins. See for example: Beucher S, Lantuejoul C. Use of

Watersheds in Contour Detection. In: International Workshop on Image Processing: Real-time Edge and Motion Detection/Estimation. Rennes, France; 1979. p. 12-21. Also other patch forming algorithms can be used. Typically, each patch has a relatively homogeneous curvature and nodes with high curvature (i.e. surface ridges) only he on the boundary between patches. For example, the maximal principal curvature of the surface, which can be computed at each vertex of the mesh, may serve as base feature for the watershed algorithm to divide the patches. In one embodiment, the costal surface is extracted by merging the patches which have a relatively small distance to the convex hull of the lung, e.g. at least one quarter of its triangles closer than 10 mm.

For visualization and analytic purposes the costal surface may be further subdivided into anterior, posterior, superior and lateral parts, e.g. based on which of the six faces (anterior, posterior, lateral, medial, superior, inferior) of the lung's bounding box are intersected by the triangles normal vectors. For example, a box projection approach can be applied wherein the lung is mapped to a box. The different sides of the box can be used as labels for the nodes associated with this side. Geometrically, rays are cast from each node in the direction of the node normal. The nodes with rays intersecting either the anterior, posterior, lateral or superior box side are assigned the corresponding label. In a post-processing step a connected component analysis can be applied to retain only the largest connected component for each of the four designated regions. Costal surface nodes which are unlabeled may e.g. receive the label of the nearest labeled node.

Extraction of the diaphragmatic and medial surface.

Due to high intersubject variability of lung and heart shape and size, the differentiation between medial and diaphragmatic lung surface may be inherently difficult. A flexible and robust approach that may deal with such difficulty, is the random walker strategy which was first described by Leo Grady ("Random walks for image segmentation." - IEEE Transactions on Pattern Analysis and Machine Intelligence. 2006;28: 1768- 1783). Also other algorithms and parameters for identifying subsurfaces can be used.

The random walker operates on arbitrary graphs consisting of a set of nodes P and a set of edges E. Given seed nodes for each target region, it outputs a probability vector for each non-seed node which expresses the probability that the node belongs to any of the specified target regions. The following four steps describe in detail, a preferred embodiment how the random walker is applied to separate the medial from the diaphragmatic surface.

(1) Graph definition - A graph is constructed from the lung surface mesh, excluding the previously determined costal surface. Mesh triangle center points are used as graph nodes (pi £ P) and graph edges (e_¾ £ E) are established between adjacent triangles. After assigning labels to the graph nodes based on the presented procedure, node labels are transferred back to the corresponding triangles of the original mesh.

(2) Edge weight computation - Each edge e_¾ of the graph is assigned a weight Wij between 0 and 1 expressing the probability that node pi and pj belong to the same target region. The applied weight function wij = exp(-di_j ²) uses a distance measurement d j which is adopted from Zhang et al. ("Interactive Mesh Cutting Using Constrained Random Walks." - IEEE Transactions on Visualization and Computer Graphics. 2011; 17(3):357- 367.). In one embodiment, it takes the Euclidean distance between the nodes Pi and pj , the deviation of their surface normals m and n, , and the normal curvature Κρφ_} into account:

The coefficients wi,W2 and W3 are for example set according to Zhang et al. ("Interactive Mesh Cutting Using Constrained Random Walks." - IEEE Transactions on Visualization and Computer Graphics.

2011; 17(3):357-367.). Contrary to the adopted method where concave edges are considered relevant partition boundaries, in our case convex edges need to be emphasized which can be achieved for example by setting:

„¼)— with normal ai v&iure κ {2}

(3) Initial seed node selection - The initial seed node selection can be based on anatomical knowledge and the previously defined mesh patches from the watershed transform (see previous section "Extraction of the costal surface"). After patch generation, patches are selected which can be assigned to one of the two target regions with high certainty. In a preferred embodiment, one or more of the following criteria is used. The most posterior patch is assigned to the label diaphragmatic surface. The most superior patch and all patches which are at least partially occluded when viewed from a point centrally located below the lung are assigned to the medial surface. To further decrease the number of patches with an unknown membership, an unlabeled patch acquires the label of an adjacent patch if their normals diverge less than 30°and if no high ridge separates the patches into different regions (based on average weight Wij between nodes of both patches). This procedure can be repeated until no more patches can be labeled. In the present study, in 15 out of 32 lung meshes, all patches were labeled, rendering the succeeding random walker step redundant. At last, a single seed node is obtained from every labeled patch by determining the graph node which is associated with the central triangle of the respective patch.

(4) Iterative random walker surface partitioning - Based on the edge weights defined in step (2) and initial seed nodes from step (3), the random walker algorithm is executed and the resulting probability maps are thresholded with high values (pdiaphragm > 0.8 and p_mediai > 0.6) in order to derive a first partitioning of the diaphragmatic and medial surface. In the next iteration, seeds are randomly placed in the last-obtained partitions and added to the initial set of seed nodes. The iteration process is terminated when the partitions no longer change significantly (number of triangles is partition changed less than 3% to previous partition). Ultimately, the final partitioning results from assigning nodes with Pdiaphragm > 0.8 to the diaphragmatic surface and nodes with Pdiaphragm < 0.8 to the medial surface. Motion analysis and visualization

Using the deformation field and the partitioned surface of the lung, it is possible to extract for example one or more of the following features to assess the motion of different respiratory muscles and highlight the diaphragmatic weakening caused by Pompe disease. In the present study, all features listed below are extracted from both left and right lung separately.

Basic features.

In one embodiment, for each time point t G [1...T ] the volume V

(Mt^d) is computed from the lung mesh. For example, the volume can be computed using the method proposed in Zhang et al. ("Efficient feature extraction for 2D/3D objects in mesh representation." In: Proceedings 2001 International Conference on Image Processing, vol. 3. Thessaloniki, Greece; 2001. p. 935-938.)

In another or further embodiment, other values are computed such as time-dependent volume measurement, image-derived values for inspiratory volume (V (Mi^d )), residual volume (Bl: V (Μτ ), also referred to as expiratory volume) and vital capacity (B2: V (Mi^d) - V (Μτ^ά)).

Furthermore, the ratio between total lung capacity and residual volume (B3: V (Mi^d )=V (M_T ^d )) can be determined.

In another or further embodiment, as a coarse estimate of independent muscle contribution, the change of lung size is determined separately in the craniocaudal (CC) direction, which is mostly driven by the diaphragm, and in the anteroposterior (AP) direction, which is mostly determined by motion of the chest wall. The CC size is defined by the vertical distance between the most superior point of the whole lung (lung apex) and the most superior point of the diaphragmatic surface (diaphragm apex). The AP size is defined as the distance between the most anterior point in the lung surface to the most posterior point in the AP direction. Time-dependent size measurements are derived by computing the CC and AP size of the lungs at each time point t G [1...T ]. In another or further embodiment, the increase in CC and AP size is computed (B4:

CC(Mi^d)=CC(M_T ^d) and B5:

and in a further embodiment, the ratio between these two features (B6: B4/B5) can also be computed.

Advanced features.

Diaphragm and chest wall contribution can be measured by directly computing the volume displaced by the diaphragmatic surface and the costal surface. At first, the set of triangles D c F and C c F,

representing the diaphragmatic and costal surface respectively, are eroded by two triangle strips to eliminate unreliable measurements at the edges. The volume displaced by a single triangle f (f = AviV2V3 e F and Vi £ V ) is obtained by computing the volume of the polyhedron that connects the vertices of the triangle at time point t with the vertices of the same triangle at time point t + 1. The polyhedron volume can be negative or positive depending on whether the triangle moves in the direction of its normal vector or not. The volume displaced by the whole diaphragmatic surface between time points t and t + 1, denoted as AVD*, is computed by

accumulating the polyhedron volumes of all triangles in D:

At last, the volume which is displaced during the entire breathing maneuver is obtained by summing the displaced volumes AVD* over all time points.

The volume displaced by the costal surface, AVc* and AVc (A2), is computed in an analogous manner. Note that for inter-subject comparison, the features Al and A2 are normalized based on the residual volume in order to compensate for differences between subjects (e.g. related to age, height and gender). Furthermore, we define the contribution of the diaphragm as:

It is found that in healthy subjects the cranial excursion of the diaphragm is posteriorly greater than anteriorly. This may especially apply to the right hemidiaphragm (half of the diaphragm). Excursion variations within the diaphragm can be observed with various imaging modalities, such as ultrasonography, fluoroscopy and MRI. In order to quantify anterior-posterior kinematic variations in the present data, the mesh triangles of the diaphragmatic surface D can be divided into two groups based on their position at full-expiration. For example, the mid-coronal plane of the lung bounding box can be utilized as a threshold. Consequently, the maximum CC displacement of each triangle can be as calculated with respect to the full-inspiration state and a weighted average was computed for the anterior and posterior part separately, taking the triangle area into account. For example, the ratio between mean anterior and posterior diaphragm displacement can be calculated for the diaphragmatic surface of the right lung (A5). The left lung may be discarded from this analysis since a large portion of the left hemidiaphragm is typically missing due to the heart position. Additionally, the mean displacement of the whole

diaphragmatic surface can be reported (A4).

In another or further embodiment, the orientation of the right diaphragm dome can be determined as a surrogate measurement of the whole diaphragm shape. For example, at each time point, a weighted average of triangle orientation (angle between normal vector and CC axis) is computed in the sagittal plane, using the triangle area as weights.

Diaphragm orientation at full-inspiration (A6) and the largest signed difference between the orientation at the start of motion and any time point between start and end of motion (A7) can also be obtained. The start and end of motion can be defined e.g. as the time points at which 10% and 90% of the volume change was reached, respectively. FIG 4 shows example results of deviation between manual and automatic lung size computation. The boxplots summarize the ASE and MARE error for CC (left plot) and AP (right plot) lung size

measurements. Measurements are compared between observer 1 (01), observer 2 (02), mean of observer 1 and 2 (mO), and the automatic method (A). Median values are displayed as horizontal lines in the box, the box limits express the 25th and 75th percentile and outliers are plotted as crosses.

In the present study, for the B-spline based 3D+t forward transformation isotropic B-spline grid resolutions of 10, 20, 40 and 80 mm in spatial domain and a spacing of 1, 2, 3 and 5 frames in time domain were evaluated. The computed transformations were used to transform the manual segmentations Mi^d and Μτ of the first and last frame (respectively) to the mean reference frame. For every subject the Dice overlap scoring between both masks was computed and mean values and standard deviations were calculated for the whole population. A perfect overlap of both masks would result in a Dice score of 1. The results show that the best spatial resolution is 20x20x20 mm which gives a mean Dice overlap of 0.87±0.015. However, temporal spacing has very little effect on the deformation fields of the first and last frame (the only time-points with manual segmentations) and could not be optimized in this fashion. Visual inspection of the registration results showed that with a very large temporal spacing (e.g. 5 frames) the start and end of the main respiratory motion are not captured correctly. On the other hand, temporal smoothing is necessary to compensate for ghosting artifacts. A temporal spacing of 2 frames was chosen as a trade-off between accuracy and ghosting compensation. The presented method for lung motion estimation was evaluated quantitatively by comparing AP and CC lung size measurements based on the 3D+t lung meshes and landmarks extracted manually from the dynamic image sequence by two observers (Mogalle and Perez-Rovira). One observer selected the axial shce in the middle between lung apex and highest point of the whole diaphragm at full-expiration. The maximal AP lung size was then assessed by both observers using this predefined subject-specific slice. The CC lung size was measured from the lung apex to the diaphragm apex. For every subject with T time points, the disagreement between observations f(t) and g(t) is split into two aspects: The absolute systematic error, that captures the systematic bias between observations and is measured at full- expiration as:

and the mean absolute residual error, that captures the variations once t

The figure shows the inter-observer difference and the difference between the automatic method and the mean observer. When determining the CC size at full-expiration, the two observers differ only marginally (mean 0.76 mm, median 0 mm) indicating that the diaphragm edge is easily identifiable by a human observer. The highest point of the diaphragm is determined slightly differently for the automatic method compared to the manual annotation, which results in a median ASE of 1.57 mm (less than a voxel). In both comparisons the median MARE lies below one voxel (inter- observer: 0.99 mm, mean observer vs. automatic method: 1.23 mm). A paired, two-sided Wilcoxon signed rank test with significance level of 5% showed that there were statistically significant differences between the interobserver variability and the manual vs. automatic error. The exact and reliable identification of the anterior lung edge in this study was more difficult due to strong ghosting artifacts occurring at the chest region, resulting in a median interobserver ASE of 2.8 mm.

Median ASE between the automatic method and the mean observer was 3.6 mm. Median inter-observer MARE was computed as 1.9 mm and median MARE of mean observer vs. automatic method as 2.3 mm. For the AP size measurement, no significant difference was observed between the interobserver error and the error between the automatic method and the mean observer (PASE = 0.26, PMARE = 0.41).

FIG 12 shows a table with example results for the proposed image-based features and the spirometry features. The extracted features indicate differences between diseased and controls. In patients where the diaphragm is severely affected by Pompe disease, almost no change in CC size (B4) can be observed during exhalation, indicating an impaired diaphragm. In 6 out of 10 patients the CC size at full-inspiration is less than 10% increased versus the full-expiration state. In contrast, 5 out of 6 controls show an increase in CC size of at least 30%. The increase of AP size (B5) is however not significantly different between patients and controls and ranges from 10% to 36%. This indicates that Pompe patients conserve the capabilities to partially inflate the lung using the chest muscles despite having an impaired diaphragm. The CC-AP -ratio (B6) shows that 5 out of 6 controls have a larger increase in CC than in AP size, while only 4 out of 10 patients present the same behavior.

FIG 5 shows example results of chest wall and diaphragm contribution to overall lung volume change. The bar plot shows the amount of volume displaced by the costal surface (upper bars) and

diaphragm surface (lower bars) for all subjects. Patients (P01-P10) and controls (C01-C06) are sorted within their group in descending order with respect to supine FVC (% of predicted).

FIG 6 shows example results of diaphragm orientation and volume displacement during exhalation. The orientation of the right diaphragm dome (black line: mean value, gray area: variance) is plotted for the exhalation maneuver of two patients (P01,P03) and a healthy control (C02). Start and end of motion are indicated by vertical, dotted lines. Frame- to-frame volume displacement (dashed red line: costal surface, dashed blue line: diaphragm surface) and total lung volume change (solid green line) are accumulated over the sequence and normalized with respect to total lung capacity.

The volume that is displaced by the diaphragmatic surface and the costal surface is presented both integrated over the whole maneuver (Fig. 5), and as a time-dependent cumulative increase (Fig. 6). In all but one subject (C05), the total volume displacement was higher for the costal surface (A2) than the (smaller) diaphragmatic surface (Al). Fig. 5 shows that in 8 out of 10 patients less than 217 ml were displaced by the

diaphragm, whereas all controls displaced more than 617 ml. In contrast, the volume displaced by the costal surface is not significantly different between patients and controls, ranging from 581 ml to 1457 ml and 467 ml to 1986 ml respectively. The contribution of the diaphragm to volume change (A3) is significantly smaller in patients than in controls. In patients, lung volume changes are thus strongly based on chest wall movement. The table in FIG 12 presents further statistics for the features A1-A3.

Evaluation of basic image features (B1-B6), advanced image features

(A1-A7) and PFT features (S1-S6). Mean values and standard deviations are displayed for patients and controls separately. The two-sided Wilcoxon rank sum test (with 5% significance level) indicates significant differences between patient and control group. The table is sorted in descending order by p-value. Normalization with respect to residual volume was performed for features labelled with a star (*).

In Fig. 6 the changing orientation of the diaphragm dome is shown for three selected subjects. Both a representative patient (P03) and control (C02) are chosen as well as a patient with relatively good respiratory performance (P01). As seen in the example of subject C02, the diaphragm dome is tilted backwards at full-inspiration in healthy subjects. In the group of controls, the average diaphragm orientation (A6) was 22°±5° (range: 15°-29°). During exhalation the diaphragm rises and the posterior part of the muscle abuts the chest wall. The angle indicating diaphragm dome orientation thus decreases, resulting in a maximum orientation difference (A7) of -14°±6° (range: -23°— 8°) in the healthy control group. In Pompe patients the mean orientation at full -inspiration was 4°±5° (range: -2°-17°). Compared to the healthy controls, the orientation change during exhalation behaved in the opposite manner for 8 out of 10 patients. On average, the maximum orientation change with respect to the start of exhalation was 4°±5° (range: -5°-10°). Patients P04 and P01 showed characteristics of both groups (Fig. 6). The course of diaphragm orientation was similar to healthy controls but diaphragm orientation at full -inspiration was considerably reduced with respect to healthy subjects.

FIG 7 shows example results of cranial excursion of the right hemidiaphragm: anterior vs. posterior. The light and dark boxplots for each subject represent the cranial excursion of all points in the anterior and posterior part of the right diaphragm dome, respectively.

Significant differences between anterior and posterior excursion are marked with triangles (down-pointing triangle: larger anterior excursion, up- pointing triangle: larger posterior excursion) and circles indicate no significant difference. The maximum craniocaudal excursion of diaphragm points was investigated. Average cranial displacement (A4) in healthy controls was 61 mm (range: 51-81 mm) and 14 mm (range: -4-48 mm) in Pompe patients. The hypothesis that anterior and posterior part of the diaphragm move differently in patients and healthy controls (A5), was investigated for the right hemidiaphragm (Fig. 7). In accordance with findings from literature, the median displacement was significantly larger in the posterior than in the anterior diaphragm part (difference in the range of 1.7-12.8 mm). The patient group presented the opposite behavior such that displacement is either in the caudal direction (P10, P09), no significant difference between anterior and posterior diaphragm parts can be observed (P04), or the anterior part moves significantly more cranially than the posterior part (P01-P03, P05-P08). FIG 8 show example results of maximum excursion of points on costal and diaphragm surface. The outer six segments of subjects C02, POl and P03 represent the left, right, left/right-anterior and left/right-posterior parts of the costal surface. The maximum amplitude of surface excursion is measured in right, left, posterior and anterior direction, respectively. The left and right diaphragm surface segments of subjects C02, P03, POl and P02 display cranial excursion. More displacement corresponds to a darker color. The dashed lines in P02 indicate the separation line for the anterior-posterior comparison presented in Fig. 7.

The 2D displacement maps in Fig. 8 enable the visual

investigation of maximal excursion in all relevant lung surface regions. Due to the supine position during the scan, the subjects back is resting on the scanner table which leads to a minimal movement of the posterior costal surface. The three displayed maps represent common motion patterns observed in the study population. Patients with strongly reduced pulmonary function (e.g. P03) show less motion in the diaphragm parts than in the anterior costal surface. Even paradoxical motion can occur, which means that diaphragm points move caudally during exhalation (P03 in Fig. 8). In healthy controls, the diaphragmatic surface displaces the most. Especially in the right hemidiaphragm a gradient can be observed with maximum displacement in the posterior and lateral parts, and less displacement in the anterior part. In patient P01, overall displacement of the diaphragmatic surface is comparable with the healthy controls. However, the posterior parts in both left and right hemidiaphragm show less movement than the anterior parts.

In some embodiments, a visual representation is constructed to show the magnitude of displacement over the whole lung surface. For example, to generate a 2D representation, the diaphragmatic surface and the anterior, posterior, left lateral and right lateral parts of the costal surface can be associated with the bottom, front, back, left and right side of a rectangular box respectively. Subsequently, the box may be unfolded by rotating its sides into the bottom plane resulting in a 2D map which gives an overview of the entire lung surface. The maximum excursion of each point along one of the three major axes (LR, AP, CC) is visualized in color on the 2D surface maps. For example, the eight displayed lung surface parts can show one or more of the following type of movement:

• lateral lung surface parts: maximum excursion along LR axis

(positive towards mid-sagittal plane),

• anterior/posterior lung surface parts: maximum excursion along AP axis (positive towards mid-coronal plane),

· diaphragm parts: maximum excursion along CC axis (positive in

cranial direction).

The relation between pre-scan spirometry features (S1-S6), basic image features (B1-B6) and advanced image features (A1-A7) was investigated by computing the Pearson correlation coefficient between all combinations of these features over the 16 subjects. From this it was found that the spirometry features FVC and FEV1 highly correlate with each other (r=0.95), indicating that similar clinical conclusions can be drawn from both measurements. On the other hand, it was found that he residual volume (B 1) does not correlate well with any other features, likely because it is highly dependent on the subject's gender, height and age. However, the image-derived volume ratio (B3: TLC/RV) correlates well with FVC (r=0.88), indicating that both features measure a similar aspect of pulmonary function.

The image-based features related to chest wall movement (B5,A2) correlated best with MIP (r=0.67, r=0.56). Compared to other spirometry features, only poor correlation was observed between MIP and features describing diaphragm motion (r<0.65). Very high correlation (r=0.94, r=0.95, r=0.97) was found between the image-derived features CC size ratio (B4), volume displaced by the diaphragm (Al) and mean CC diaphragm excursion (A4). All of these features use different computation methods, yet all express the overall diaphragmatic motion, revealing diaphragm impairment in Pompe patients (see table in FIG 12). Of the image-based features, Al and A4 correlate best with spirometry features. However, the table in FIG 12 shows that A5 and A7 have better discriminative capabilities. This observation suggests that spirometry and image-based features may capture different aspects of diaphragm weakness.

The present study shows that dynamic 3D MRI of the thorax is suitable to evaluate the functional interaction of respiratory muscles during breathing. Automatic, quantitative analysis of diaphragm motion during slow exhalation was introduced to assess diaphragm weakness and characterize respiratory muscle involvement in patients with Pompe disease. Of the proposed features, diaphragm orientation and displacement difference between anterior and posterior parts of the diaphragm were the features with highest sensitivity to subtle pathophysiological changes. Our image-based results indicate that at an early disease stage patients can compensate for diaphragm weakness with increased chest wall movement leading to relatively normal pulmonary function. For instance, patients P01 and P02 exhibit the highest supine FVC in the patient group (69% and 62% of predicted values respectively), however patient P02 extensively uses rib cage expansion to maximally increase lung volume during inspiration as seen in Fig. 5.

In some embodiments, the calculation of CC and AP lung expansion can enable the assessment of general diaphragm and chest wall motion, though this method may not always be sensitive enough to detect early pathophysiological changes. For instance, patient P01 and control subject C02 both present a CC lung size increase of 30% and comparable AP expansion of 17% and 21%, respectively. Still, by evaluating the diaphragm orientation (Fig. 6) and diaphragm displacement feature (Fig. 8), moderate diaphragm weakness could be detected in patient P01. Note that all proposed image-based features can be obtained for left and right lung separately, e.g. revealing significant different elevation of right and left hemi diaphragm in patient P02.

Next to dynamic MRI, ultrasound can also be used as a radiation- free technique to quantitatively evaluate diaphragm kinematics. Some of the frequently presented advantages are its relatively low costs, portability, accessibility, applicability to a wider range of subjects and high temporal resolution. However, in order to assess kinematics of the respiratory system it is desirable to evaluate chest wall and diaphragm motion simultaneously and especially analyze diaphragm motion at defined locations in a reproducible manner, which can be achieved with MRI.

The presented study may provide further insight compared to static breath-hold MRI images. Based on manual segmentations, significant differences in CC and AP lung expansion may be observed between patients and controls. The present methods enable measuring properties of the diaphragm and lung directly and more accurately due to the 3D mesh representation of the entire lungs. Furthermore, it allows intermediate respiratory phases to be investigated in and the movement of different parts of the diaphragm can be tracked and evaluated.

FIGs 9A-9C illustrate steps of an embodiment for quantifying lung kinematics.

One embodiment comprises receiving a sequence of three- dimensional images Ixyz of a lung motion recorded as a function of time t during a breathing maneuver. Another or further embodiment comprises calculating, based on the three-dimensional images Ixyz, a sequence of three dimensional lung surface meshes Lxyz representing the recorded lung motion. A first lung surface mesh LI can be selected at a first recording time tl in the sequence of three dimensional lung surface meshes Lxyz. Part of the selected first lung surface mesh LI can be identified as a first

subsurface mesh SI.

In another or further embodiment, a second lung surface mesh Lend is selected at a second recording time tend in the sequence of three dimensional lung surface meshes Lxyz. Part of the second lung surface mesh Lend is identified as a second subsurface mesh Send. The identification of the second subsurface mesh Send corresponds to the identification of the first subsurface mesh Si propagated through time [tl...tend] during the recorded lung motion. In some embodiments, a kinematic parameter A is calculated based on a comparison of the respective corresponding subsurface meshes Si and Send at different recording times tl and tend.

In one embodiment, the method comprises identifying, within the sequence of three dimensional lung surface meshes, an end-inspiration lung surface mesh and an end-expiration lung surface mesh. For example, the first recording time tl coincides with end-inspiration of the lung motion where lung volume is maximal For example, the second recording time tend corresponds coincides with end-expiration of the lung motion where lung volume is minimal. Accordingly, the kinematic parameter can be calculated based on a comparison of the respective subsurface meshes at recording times tl,tend coinciding with end-inspiration and end-expiration of the lung motion.

In one embodiment, the respective subsurface mesh is a diaphragm surface mesh corresponding to a moving part of the lung surface associated with diaphragm movement. In a further embodiment, the kinematic parameter is based on a diaphragm contribution A3 to displaced volume obtained by calculating a volume between the diaphragm surface at end-inspiration and the diaphragm surface at end-expiration. In another or further embodiment, the kinematic parameter is based on a maximum diaphragm orientation change A7 obtained by calculating a maximum of angle differences between the diaphragm surface at end-inspiration and the diaphragm surface end-expiration. Preferably, the kinematic parameter is normalized by a total displaced lung volume, e.g. obtained by calculating a difference between the maximum total volume enclosed by the lung surface mesh at end-inspiration and the minimum total volume enclosed by the lung surface mesh at end-expiration.

One embodiment comprises calculating total lung volumes enclosed by the three dimensional lung surface meshes Lxyz and

identifying, within the sequence of three dimensional lung surface meshes an end-inspiration lung surface mesh wherein a maximum total volume is enclosed, and an end-expiration lung surface mesh wherein a minimum total volume is enclosed.

FIG 10 illustrates a flow diagram of an example procedure to distinguish the costal surface convex from the rest of the lung surface concave, namely the base and mediastinal surface. FIG 11 illustrates a flow diagram of an example procedure for automatically segmenting the lung base and mediastinal surface in the concave parts of the lung surface mesh.

In one embodiment, identification of the first subsurface mesh S 1 comprises division of the first lung surface mesh LI into a costal surface mesh and a residual surface; and subdivision of the residual surface into a diaphragm surface and a mediastinal surface. For example, the costal surface is identified based on a convexity of the respective subsurface. In another or further embodiment, the diaphragm surface is identified based on a concavity of the respective subsurface.

In some embodiments, subdivision of the residual surface into a diaphragm surface and a mediastinal surface is based on a random walker algorithm. For example, the random walker algorithm starts with initial seed points on the residual surface assigned to the mediastinal surface or diaphragm surface. For example, the initial seed points are based on concave patches generated by dividing the first lung surface mesh LI using a watershed transform. In one embodiment, the initial seed points are automatically generated based on a criterion that the most superior patch belongs to the mediastinal surface and the most posterior patch on the residual surface belongs to the diaphragm surface. In another or further embodiment, the initial seed points are automatically generated based on a criterion that patches of the diaphragm surface are entirely visible when viewing the residual surface from the inferior side.

In one embodiment, patches are expanded to adjacent patches based on a criterion that the adjacent patches are likewise oriented with surface normal angles smaller than a threshold angle, e.g. less than thirty degrees. In another or further embodiment, patches are expanded to adjacent patches based on an edge weight cost. For example the edge weight cost between neighboring patches is calculated based on distance between vertices at an edge between the neighboring patches, a difference between surface normal of the neighboring patches, and a curvature of the

neighboring patches being concave or convex. In one embodiment, initial seed points are selected by identifying a triangle in the center of each patch.

In one embodiment, a random walker algorithm is used to generate a probability map indicating a probability that a part of the residual surface belongs to the diaphragm surface and mediastinal surface. In another or further embodiment, subsequent seed points are generated based on the probability map and the random walker algorithm is executed again including the subsequent seed points. For example, iteration of the random walker algorithm is terminated when the result does not change significantly.

In a preferred embodiment, the method for quantifying lung kinematics, comprises one or more steps of

- calculating, based on a sequence of three-dimensional images, a three dimensional lung surface mesh and a transformation which models deformation of the three dimensional lung surface mesh to represent a recorded lung motion as a function of time;

- calculating a first lung surface mesh e.g. based on application of the transformation to the three dimensional lung surface mesh at a first time;

- identifying part of the first surface mesh as a first subsurface mesh, e.g. based on criteria of the mesh itself and/or based on different parts of the sequence;

- identifying part of a second lung surface mesh at a second time as a second subsurface mesh, wherein the identification of the second subsurface mesh corresponds to the identification of the first subsurface mesh propagated through time during the recorded lung motion, e.g. based on application of the transformation to the three dimensional lung surface mesh and/or the subsurface mesh at the second time. In one embodiment, the identifying of the second subsurface mesh Send comprises calculating a deformation field T which models the lung motion through time t; and applying the deformation field to the first subsurface mesh Si to propagate identification of the first subsurface mesh Si from the first recording time tl to the second recording time tend to obtain identification of the second subsurface mesh Send corresponding to the first subsurface mesh Si. For example, the deformation field T is calculated based on a registration algorithm which models the lung motion according to a non-rigid transformation model. For example, the registration algorithm utilizes a four-dimensional free-form B-spline deformation model.

In some embodiments, the registration algorithm is based on an optimization problem where a cost function C is minimized with respect to parameters that describe a transformation Tu between a first three- dimensional image II in the sequence of three-dimensional images Ixyz and a reference three-dimensional image. In other or further embodiments, transformation between the first three-dimensional image II and a second three-dimensional image lend is calculated by concatenation of

transformation between the first three-dimensional image II to the reference three-dimensional image and inverse transformation between the second three-dimensional image and the reference three-dimensional image. For example, the reference three-dimensional image is based on an average of the sequence of three-dimensional images Ixyz. Preferably, the cost function is based on a variance of intensity values of voxels in the three- dimensional images Ixyz over time t.

Preferably, the sequence of three dimensional lung surface meshes is calculated by steps including preprocessing the three-dimensional images for correcting intensity inhomogeneities. In some embodiments, a mask is used to focus registration on a specific part of the input image.

The present methods can be embodied on a non-transitory computer-readable medium with software instruction that when executed by one or more computers cause the one or more computers to execute a respective method as described herein.

Also an apparatus for quantifying lung kinematics is envisaged. For example, the apparatus comprises the non-transitory computer-readable medium. Alternatively, or in addition, the apparatus may comprise means for recording the sequence of three-dimensional images. For example, the sequence of three-dimensional images is recorded by Magnetic Resonance Imaging MRI. For the purpose of clarity and a concise description, features are described herein as part of the same or separate embodiments, however, it will be appreciated that the scope of the invention may include

embodiments having combinations of all or some of the features described. For example, while the present study and analysis is focused on a specific patient group, the proposed methods are applicable to a large variety of medical cases where the diaphragm function itself, the effect of other conditions on the respiratory muscles or general lung deformation needs to be investigated. Other areas of application for this kind of kinematic analysis are, for instance, other neuromuscular diseases affecting the diaphragm, congenital diaphragm hernia, chronic obstructive pulmonary disease, damage of the phrenic nerves, and adolescent idiopathic scoliosis. Of course, it is to be appreciated that any one of the above embodiments or processes may be combined with one or more other embodiments or processes to provide even further improvements in finding and matching designs and advantages. It is appreciated that this disclosure offers particular advantages to analysis of diaphragm movement, and in general can be applied also for other subsurfaces of the lungs.

The present disclosure may improve automated methods to quantify the contribution of the different muscles involved in the respiratory system e.g. in order to get deeper insights into the contribution of these muscles in various patient group scomp are d to healthy controls. In some embodiments, dynamic 3D MRI is used to capture a slow full exhalation breathing maneuver. In further embodiments, a registration-based approach is applied to obtain deformation fields describing the motion between frames of the image sequence. For example, lung segmentations can be propagated from a reference frame to all other frames, to facilitate quantification and analysis of the movement of the diaphragm and chest wall during the respiratory maneuvers. Contributions may include one or more of 1) the optimization of a 3D+t registration framework for lung motion estimation, 2) the automatic extraction of lung surface regions adjacent to the diaphragm and chest wall, 3) the computation of geometry and motion features for different lung surface parts and 4) the evaluation of the relationship between the proposed image features and PFT features as well as their ability to discriminate between patients with Pompe disease and healthy controls.

In the present disclosure we have shown that automatic analysis of 3D+t MRI sequences allows us to quantify respiratory dynamics in healthy and diseased subjects. By using a registration-based segmentation approach, lung surface extraction, and lung surface partitioning, we were able to obtain image-based features that quantify the contribution of the different muscles groups involved in the breathing maneuver. Experimental results show that several of the proposed features, especially features describing the shape change of the diaphragm (including diaphragm orientation and displacement differences within the diaphragm) enable a clear differentiation between Pompe patients and healthy controls. More importantly, the proposed image-based features have been shown to have the potential to be more sensitive to moderate diaphragm weakness than pulmonary function tests. The latter can only measure the combined contribution of all muscles involved with respiration, while image-based analysis can measure the contribution of the different muscle groups separately.

Finally, the above-discussion is intended to be merely illustrative of the present systems and/or methods and should not be construed as limiting the appended claims to any particular embodiment or group of embodiments. The specification and drawings are accordingly to be regarded in an illustrative manner and are not intended to limit the scope of the appended claims. In interpreting the appended claims, it should be understood that the word "comprising" does not exclude the presence of other elements or acts than those listed in a given claim; the word "a" or "an" preceding an element does not exclude the presence of a plurality of such elements; any reference signs in the claims do not limit their scope; several "means" may be represented by the same or different item(s) or implemented structure or function; any of the disclosed devices or portions thereof may be combined together or separated into further portions unless specifically stated otherwise. The mere fact that certain measures are recited in mutually different claims does not indicate that a combination of these measures cannot be used to advantage. In particular, all working combinations of the claims are considered inherently disclosed.

Claims

1. A computer-implemented method for quantifying lung kinematics, the method comprising

- receiving a sequence of three-dimensional images (lxyz) of a lung motion recorded as a function of time (t) during a breathing maneuver;

- calculating, based on the three-dimensional images (lxyz), a sequence of three dimensional lung surface meshes (Lxyz) representing the recorded lung motion;

- selecting a first lung surface mesh (LI) at a first recording time (ti) in the sequence of three dimensional lung surface meshes (Lxyz);

- identifying part of the first surface mesh (LI) as a first subsurface mesh (S i);

- selecting a second lung surface mesh (L_end) at a second recording time (tend) in the sequence of three dimensional lung surface meshes (Lxyz);

- identifying part of the second lung surface mesh (Lend) as a second subsurface mesh (Send); wherein the identification of the second subsurface mesh (Send) corresponds to the identification of the first

subsurface mesh (S i) propagated through time (ti...t_end) during the recorded lung motion; and

- calculating a kinematic parameter (A) based at least on the respective corresponding subsurface meshes (S i, Send) at different recording times (ti,tend).

2. The method according to claim 1, wherein identification of the first subsurface mesh (Si) comprises

- division of the first lung surface mesh (LI) into a costal surface mesh and a residual surface; and

- subdivision of the residual surface into a diaphragm surface and mediastinal surface.

3. The method according to claim 2, wherein the costal surface is identified based on a convexity of the respective subsurface and/or wherein the diaphragm surface is identified based on a concavity of the respective subsurface.

4. The method according to claim 3, wherein subdivision of the residual surface into a diaphragm surface and a mediastinal surface is based on a random walker algorithm, wherein the random walker algorithm starts with initial seed points on the residual surface assigned to the mediastinal surface or diaphragm surface, wherein the initial seed points are based on concave patches generated by dividing the first lung surface mesh (LI) using a watershed transform, wherein the initial seed points are automatically generated based on

- a criterion that the most superior patch belongs to the

mediastinal surface and the most posterior patch on the residual surface belongs to the diaphragm surface and/or

- a criterion that patches of the diaphragm surface are entirely visible when viewing the residual surface from the inferior side.

5. The method according to claim 4, wherein patches are expanded to adjacent patches based on a criterion that the adjacent patches are likewise oriented with surface normal angles smaller than a threshold angle, e.g. less than thirty degrees.

6. The method according to claim 5, wherein patches are expanded to adjacent patches based on an edge weight cost, wherein the edge weight cost between neighboring patches is calculated based on distance between vertices at an edge between the neighboring patches, a difference between surface normal of the neighboring patches, and a curvature of the neighboring patches being concave or convex.

7. The method according to claim 6, wherein a random walker algorithm is used to generate a probability map indicating a probability that a part of the residual surface belongs to the diaphragm surface and mediastinal surface, wherein subsequent seed points are generated based on the probability map and the random walker algorithm is executed again including the subsequent seed points, and wherein iteration of the random walker algorithm is terminated when less than a threshold percentage change is achieved.

8. The method according to any of the preceding claims, wherein the identifying of the second subsurface mesh (Send) comprises

- calculating a deformation field (T) which models the lung motion through time (t); and

- applying the deformation field to the first subsurface mesh (Si) to propagate identification of the first subsurface mesh (Si) from the first recording time (tl) to the second recording time (tend) to obtain

identification of the second subsurface mesh (Send) corresponding to the first subsurface mesh (Si).

9. The method according to any of the preceding claims, wherein a deformation field (T) is calculated based on a registration algorithm which models the lung motion according to a non-rigid transformation model, wherein the registration algorithm is based on an optimization problem where a cost function (C) is minimized with respect to parameters that describe a transformation (Τμ) between a first three-dimensional image (II) in the sequence of three-dimensional images (Ixyz) and a reference three- dimensional image.

10. The method according to any of the preceding claims, wherein transformation between the first three-dimensional image (11) and a second three-dimensional image (lend) is calculated by concatenation of

transformation between the first three-dimensional image (11) to the reference three-dimensional image and inverse transformation between the second three-dimensional image and the reference three-dimensional image, wherein the reference three-dimensional image is based on an average of the sequence of three-dimensional images (Ixyz).

11. The method according to any of the preceding claims, wherein the respective subsurface mesh is a diaphragm surface mesh corresponding to a moving part of the lung surface associated with diaphragm movement, wherein the kinematic parameter is based on

- a diaphragm contribution (A3) to displaced volume obtained by calculating a volume between the diaphragm surface at end-inspiration and the diaphragm surface at end-expiration and/or

- a maximum diaphragm orientation change (A7) obtained by calculating a maximum of angle differences between the diaphragm surface at end-inspiration and the diaphragm surface end-expiration.

12. The method according to any of the preceding claims, wherein the kinematic parameter is normalized by a total displaced lung volume obtained by calculating a difference between the maximum total volume enclosed by the lung surface mesh at end-inspiration and the minimum total volume enclosed by the lung surface mesh at end-expiration.

13. The method according to any of the preceding claims, wherein the kinematic parameter is calculated based on a comparison of the respective subsurface meshes at recording times (tl,tend) coinciding with end- inspiration and end-expiration of the lung motion.

14. A non-transitory computer-readable medium with software instruction that when executed by one or more computers cause the one or more computers to execute the method according to any of the preceding claims.

15. A system for quantifying lung kinematics, the system comprising

- the computer-readable medium according to claim 14;

- a computer for executing the software instructions stored on the computer-readable medium; and

- a display for displaying a representation of the kinematic parameter resulting from execution of the software instructions.