1 Introduction
DCEMRI renography involves injection of an MRI contrast agent into the blood stream for obtaining magnetic resonance images (MRI) of internal body structures such as kidney, liver and spleen as shown in Figure 1. The images are taken sequentially over the time, until the contrast agent is discharged. Delineation of the kidney region from these dynamic image sequence is required in order to quantify various aspects of renal function, such as filtration and perfusion or blood flow. However, Absolute quantification in DCEMRI is often obfuscated by the issues including region of interest (ROI) definition, which may be influenced by the adjacent regions giving rise to the so called Partial Volume Effect (PVE). PVE is the contamination of any single voxel with two or more signal intensities from different tissues or regions (refer Figure 1 for more explanation). The contamination is caused by the finite resolution of imaging systems that limits the higher frequencies i.e., fine details such as edges and structures. This results in producing inaccurate or blurred images especially in DCEMRI where the spatial information is sacrificed for better temporal resolution since time is an important factor in the assessment of kidney function. The two primary issues concerning PVE are: (i) the action of the Point Spread function on the image data, and (2) the action of finite voxel sampling (comb function) and the integration of signal within the volume of the voxel (rect function). Further details, albeit for the same effects in 2D can be found in Yip et al, Phys Med Biol 2011, listed under my papers on the uni web page. Therefore, DCEMRI is required to undergo several steps including: 1) suppression of PVE gutierrez2010partial , selection of kidney ROI zollner2009assessment ; rusinek2007performance and a final stage of modelling kidney function by analysis of the resulting timeintensity plots tofts2012precise .
(a)  (b) 
A templatebased approach based on tissue classification proposed in gutierrez2010partial
has demonstrated that PVE is a very significant issue in quantification of DCEMRI studies. From the literature, eigen based methods such as Independent component analysis (ICA) and its variants have also been used for obtaining the kidney ROI
zollner2007assessment. In addition, kmeans clustering
zollner2009assessment has also been employed. Although these methods have shown improvements in the processing of the data, a complete automated solution has not been provided. The major issue with the aforementioned approaches is the need for manual or semiautomatic delineation of the target region, and the need for prior knowledge of the point spread function of the particular MRI sequence used, which will vary between particular sequences and vendors. In addition, this can be labour intensive, timeconsuming and inefficient as the human expert has to examine the whole sequence of images to find the most suitable frame for delineation. Therefore, to address this problem in more holistic manner, we in this study, present a framework consisting of DMD, thresholding and blob analysis for automatic delineation of kidney ROI with the PVE compensation within the DCEMRI sequence.DMD was originally introduced in the area of computational fluid dynamics (CFD) Schmid2 , specifically for analysing the sequential image data generated by nonlinear complex fluid flows Schmid1 ; Schmid2 ; Schmid3 ; tirunagarianalysis . The DMD decomposes a given image sequence into several images, called dynamic modes. These modes essentially capture different large scale to small scale structures (sparse components) including a background structure (low rank model) DBLP:journals/corr/GrosekK14 ; randDMD ; tirunagari2016can . DMD has gained significant applications in various fields 6926317 ; 2015arXiv150804487M ; brunton2016extracting , including for detecting spoof samples from facial authentication video datasets tirunagari2015detection and for detecting spoofed fingervein images tirunagari2015windowed . The advantage of this method is its ability to identify regions of dominant motion in an image sequence in a completely datadriven manner without relying on any prior assumptions about the patterns of behaviour within the data. Therefore, it is thus potentially wellsuited to analyse a wide variation of blood flow and filtration patterns seen in renography pathology.
Our preliminary results show that DMD mode2 highlights the kidney region due to the presence of injected contrast agent. Therefore, as shown in Figure 2, this mode is then used for delineating the kidney region using thresholding and blob analysis.
Prior works directed towards tackling the issues related to Partial Volume Correction (PVE) and delineation of kidney ROI are often dealt separately using different methods. In particular, to our knowledge, there is no literature of using a single method that is able to tackle both of these issues simultaneously. Therefore, our first contribution is to introduce a single methodology that is based on DMD to tackle the aforementioned issues problem. Our second contribution is to Validate our technique for the first time on medical data with applications to DCEMRI; and finally, our third contribution is in improving the understanding of the application through our DMD framework.
The remainder of this paper is organised as follows: in Section 2, we consider the theory for DMD. Section 3 presents synthetically generated dataset as well as ten different DCEMRI datasets used in this study. In Section 4 we present our experiments and results and finally, conclusions are drawn in Section 5.
2 Methodology
Our methodological framework (Figure 3) consists of DMD, thresholding and selection of kidney region using connected component (blob) analysis.
In our previous studies Tirunagari2017 , we presented a novel automated, registrationfree movement correction approach based on windowed and reconstruction variants of Dynamic Mode Decomposition (WRDMD) to suppress unwanted complex organ motion in DCEMRI image sequences caused due to respiration.
2.1 Dynamic Mode Decomposition (DMD)
In a dynamic sequence of images P, let be the image whose size is . This image is converted to
column vector, resulting in the construction of a data matrix
P of size for images in the DCEMRI data.(1) 
The images in the DCEMRI data are collected over regularly spaced time intervals and hence each pair of consecutive images are linearly correlated. It can be justified that a linear mapping A exists between them forming a span of krylov subspace krylov1931numerical ; saad1981krylov ; ruhe1984rational :
(2) 
The krylov subspace can then be represented using two matrices and where and
.
(3) 
The mapping matrix is responsible for capturing the dynamics or the fluctuating intensities within the image sequence. The sizes of the matrices and are both each. Therefore, the size of unknown matrix would be . Unfortunately, solving for A is computationally very expensive due to it size. For instance, if an image has a size of i.e., and , the size of is then .
Since solving is computationally expensive for large image dimensions, we need an alternative solution. Our assumption that the images form a Krylov span, allows us to introduce from the standard Arnoldi iteration trefethen1997numerical ).
(4) 
Here, a companion matrix also known as a shifting matrix that simply shifts images through and approximates the last frame by linearly combining the previous images, i.e., . requires the storage of data matrix is significantly smaller than in dimensions.
(5) 
Thus, for the last frame , where is much fewer (dimension) than the dimensionality of (), one can write as a linear combination of the previous vectors. Consistent with Equations 3 and 4, we then have:
(6) 
In Schmid2
, the author describes a more robust solution, which is achieved by applying a singular value decomposition (SVD) on
. From Equation 6, SVD decomposition on subspace is calculated to obtain , and matrices that are left singular vectors, singular values and right singular vectors respectively. The inversions of these matrices are then multiplied with subspace to obtain the The fullrank matrix , determined on the subspace spanned by the orthogonal basis vectors of . described by:(7) 
Here, and are the conjugate transpose of and , respectively; and denote the inverse of the singular values . After obtaining the
matrix, the eigenvalue analysis is performed to obtain
eigenvectors and a diagonal matrix containing the corresponding eigenvalues.(8) 
It is known that the eigenvalues of approximate some of the eigenvalues of the full system . The associated eigenvectors of provide the coefficients for the linear combination that is necessary to express the dynamics within the image sequence basis. The dynamic modes are thus calculated as follows:
(9) 
Since the eigenvalues associated with are complex, it is not possible to establish the order of the dynamic modes directly. Nevertheless, we calculate the absolute value for the phaseangles associated with the complex eigenvalues and the modes with unique phaseangles are selected. Doing this will remove one of the conjugate pairs in the dynamic modes which have same phaseangles but with different signs and capture similar information sayadi2013dynamic . After discarding one of the conjugate pairs, the dynamic modes are then sorted in the ascending order of their phaseangles.
3 Dataset
In this section we briefly describe our synthetically generated data as well as the DCEMRI data.
(a)  (b) 
(c)  (d) 
3.1 Synthetic data
(1)  (2)  (3)  (4)  (5) 
(6)  (7)  (8)  (9)  (10) 
In order to validate our experiments with the groundtruth available, we artificially generated synthetic data corresponding to the DCEMRI events as that shown in Figure 4(c). That is, as the liver region denoted by in DCEMRI image (Figure 4(a)) covers both background region () and kidney region (), similarly the synthetically generated images follow the same structure as shown in Figure 4(c).
A series of images are then produced representing kidney, liver and background regions/functions as labelled in Figure 4(c) with respect to their temporal evolution shown in Figure 4
(b). The kidney region’s temporal evolution is thus modelled as a combination of Poisson and log distributions. The liver region is modelled as a sigmoid distribution and the background region as a random normal distribution. The generated regions are then convoluted using a Gaussian kernel operator (variance = 22 and block size of
). This operator thus simulates the point spread function’s (PSF’s) pixellevel mixing of the kidney regions as shown in Figure 4(d).The timeintensity plot of convolved kidney region thus shows the influence of the liver (sigmoid) and background functions as shown in Figure 6. Similarly the timeintensity plots of liver are influenced by background and kidney regions and vice versa.
3.2 DCEMRI data
The functional kidney DCEMRI datasets used in this study have been obtained from ten healthy volunteers as shown in Figure 5. These datasets are acquired after injection of of GdDTPA (Magnevist) contrast agent, on a Siemens Avanto MRI scanner, using a 32 channel body phased array coil. The MRI acquisition sequence consisted of a 3D spoilt gradient echo sequence utilising an Echo time, TE of 0.6ms, Repetition time, TR of 1.6ms, and a flip angle, FA = 17 degrees with a temporal resolution of collected for . The acquired DCEMRI datasets cover the abdominal region, enclosing left and right kidneys and abdominal aorta.
4 Experiments & Results
In this section, we present our experimental procedure as well as the results.
4.1 Evaluating the DMD framework over synthetically generated data
Synthetic data consisting of images is given as an input to DMD algorithm which produces DMD modes (In theory, for a image sequence with images, we obtain DMD modes). Recall that each dynamic mode captures a “principal dynamics” axes of the image sequence. Modes that show predominating liver, kidney and background functions are selected, as shown in Figure 6 (a). These modes are then thresholded to obtain the segmented versions of the respective functions as shown in Figure 6 (c). These segmented versions of the dynamic modes are then projected on to the convolved image sequence and mean pixel intensities of synthetic kidney, liver and background functions are computed. The timeintensity plots are shown in Figure 6 (c). It is clear that the mixing from the other functions and are totally suppressed and the timeintensity plots produced by the proposed framework recovers the original functions.
Kidney 


Liver 

Background 

(a) Dynamic modes  (b) Thresholded  (c) PVC 
4.2 Evaluating the DMD framework on DCEMRI
DCEMRI image sequences are given as an input to the DMD algorithm. The output of the DMD captures large (Figure 7(Top)) and small scale structures (Figure 7(bottom)) in the form of modes. Dynamic mode1 reveals the background (lowrank) model and the remaining modes capture the sparse representations. The contrast changes are captured in the top modes, particularly, mode2 capturing kidney region and mode3 and 4, the spleen and the liver regions respectively, as shown in Figure 7(top) on the dataset1. Noise and residuals are captured in the lower modes (Figure 7(bottom)).
Therefore, we have selected DMD mode2 across all the 10 datasets as shown in Figure 8 (a). These modes are then thresholded to obtain a binary of version (Figure 8 (b)). Later using connected component analysis on the left half of the image, the largest blob is selected as a kidney ROI as shown in Figure 8 (c). The mean of intensity values inside this automatically delineated kidney ROI are calculated to produce timeintensity plots.
(a) Dynamic mode2 across 10 datasets of DCEMRI data 
Thresholding across 10 datasets of DCEMRI data 
Delineated kidney region obtained via blob analysis across 10 datasets of DCEMRI data 
4.3 Comparison with human expert
As a baseline method a minimum bounding box region around the kidney region (3) is selected as shown in Figure 4(a) across all the datasets. To evaluate the performance of the DMD framework we compare it against the performance of the manual delineation from a human expert. From Figure 9 we see that the kidney function produced from the automatic delineation clearly compares favourably with the human expert; where as the minimum bounding box region due to the influence of the PVE deviates from the kidney function quantified by a human expert.
(1)  (2)  (3)  (4) 
(5)  (6)  (7)  (8) 
(9)  (10) 
The root mean square error (RMSE) between the quantification produced by the proposed framework and baseline method against the human expert is shown in Figure 10. It is clear that the quantification results obtained from the proposed framework has lower RMSE when compared to the baseline method.
5 Conclusions
This study shows the significance of the proposed framework based on DMD, thresholding and blob analysis as a viable automatic delineation algorithm to effectively quantify the kidney function by reducing the partial volume effect in the DCEMRI data.
Our proposed DMD framework is applied on the DCEMRI datasets collected from 10 healthy volunteers as well as the synthetic data mimicking the DCEMRI data. We found that DMD can capture functional structures like kidney, liver and spleen in the top ranked dynamic modes. Particularly dynamic mode2 capturing the kidney structure across all the datasets. This mode2 is then thresholded and kidney template is obtained using connected component analysis. The kidney function produced from our proposed framework compares favourably with the human expert; while the minimum bounding box region due to the influence of the PVE deviates. Our results demonstrate that our proposed framework is very promising in delineating and quantifying the kidneys in DCEMRI studies by removing the PVE.
Acknowledgements.
The funding for this work has been provided by the Department of Computer Science and the Centre for Vision, Speech and Signal Processing (CVSSP)  University of Surrey. We also would like to express our gratitude towards Kidney Research UK for funding the DCEMRI data acquisition as part of a reproducibility study.References

(1)
Berger, E., Sastuba, M., Vogt, D., Jung, B., Amor, H.B.: Dynamic mode decomposition for perturbation estimation in human robot interaction.
In: The 23rd IEEE International Symposium on Robot and Human Interactive Communication, pp. 593–600 (2014). DOI 10.1109/ROMAN.2014.6926317  (2) Brunton, B.W., Johnson, L.A., Ojemann, J.G., Kutz, J.N.: Extracting spatial–temporal coherent patterns in largescale neural recordings using dynamic mode decomposition. Journal of neuroscience methods 258, 1–15 (2016)
 (3) Erichson, N., Donovan, C.: Randomized lowrank dynamic mode decomposition for motion detection. Computer Vision and Image Understanding In Press (2016). DOI 10.1016/j.cviu.2016.02.005. N. Benjamin Erichson acknowledges support from the UK Engineering and Physical Sciences Research Council (EPSRC).
 (4) Grosek, J., Kutz, J.N.: Dynamic mode decomposition for realtime background/foreground separation in video. CoRR abs/1404.7592 (2014). URL http://arxiv.org/abs/1404.7592
 (5) Gutierrez, D.R., Wells, K., Diaz Montesdeoca, O., Moran Santana, A., Mendichovszky, I., Gordon, I.: Partial volume effects in dynamic contrast magnetic resonance renal studies. European journal of radiology 75(2), 221–229 (2010)
 (6) Krylov, A.: On the numerical solution of the equation by which in technical questions frequencies of small oscillations of material systems are determined. Izvestija AN SSSR (News of Academy of Sciences of the USSR), Otdel. mat. i estest. nauk 7(4), 491–539 (1931)
 (7) Mann, J., Kutz, J.N.: Dynamic Mode Decomposition for Financial Trading Strategies. ArXiv eprints (2015)
 (8) Ruhe, A.: Rational krylov sequence methods for eigenvalue computation. Linear Algebra and its Applications 58, 391–405 (1984)
 (9) Rusinek, H., Boykov, Y., Kaur, M., Wong, S., Bokacheva, L., Sajous, J.B., Huang, A.J., Heller, S., Lee, V.S.: Performance of an automated segmentation algorithm for 3d mr renography. Magnetic Resonance in Medicine 57(6), 1159–1167 (2007)
 (10) Saad, Y.: Krylov subspace methods for solving large unsymmetric linear systems. Mathematics of computation 37(155), 105–126 (1981)
 (11) Sayadi, T., Schmid, P., Nichols, J., Moin, P.: Dynamic mode decomposition of controlled hand ktype transitions. Annual Research Briefs (Center for Turbulence Research, 2013) p. 189 (2013)
 (12) Schmid, P., Li, L., Juniper, M., Pust, O.: Applications of the dynamic mode decomposition. Theoretical and Computational Fluid Dynamics 25(14), 249–259 (2011)
 (13) Schmid, P.J.: Dynamic mode decomposition of numerical and experimental data. Journal of Fluid Mechanics 656, 5–28 (2010)
 (14) Schmid, P.J., Meyer, K.E., Pust, O.: Dynamic mode decomposition and proper orthogonal decomposition of flow in a liddriven cylindrical cavity. In: 8th International Symposium on Particle Image Velocimetry, pp. 25–28 (2009)
 (15) Tirunagari, S., Poh, N., Bober, M., Windridge, D.: Windowed dmd as a microtexture descriptor for finger vein counterspoofing in biometrics. In: Information Forensics and Security (WIFS), 2015 IEEE International Workshop on. IEEE (2015)
 (16) Tirunagari, S., Poh, N., Bober, M., Windridge, D.: Can dmd obtain a scene background in color? In: Image, Vision and Computing (ICIVC), International Conference on, pp. 46–50. IEEE (2016)
 (17) Tirunagari, S., Poh, N., Wells, K., Bober, M., Gorden, I., Windridge, D.: Movement correction in dcemri through windowed and reconstruction dynamic mode decomposition. Machine Vision and Applications 28(3), 393–407 (2017). DOI 10.1007/s0013801708355. URL http://dx.doi.org/10.1007/s0013801708355
 (18) Tirunagari, S., Poh, N., Windridge, D., Iorliam, A., Suki, N., Ho, A.T.: Detection of face spoofing using visual dynamics. Information Forensics and Security, IEEE Transactions on 10(4), 762–777 (2015)
 (19) Tirunagari, S., Vuorinen, V., Kaario, O., Larmi, M.: Analysis of proper orthogonal decomposition and dynamic mode decomposition on les of subsonic jets. CSI Journal of Computing 1, 20–26 (2012)
 (20) Tofts, P.S., Cutajar, M., Mendichovszky, I.A., Peters, A.M., Gordon, I.: Precise measurement of renal filtration and vascular parameters using a twocompartment model for dynamic contrastenhanced mri of the kidney gives realistic normal values. European radiology 22(6), 1320–1330 (2012)
 (21) Trefethen, L.N., Bau III, D.: Numerical linear algebra, vol. 50. Siam (1997)
 (22) Zöllner, F.G., Kocinski, M., Lundervold, A., Rørvik, J.: Assessment of renal function from 3d dynamic contrast enhanced mr images using independent component analysis. In: Bildverarbeitung für die Medizin 2007, pp. 237–241. Springer (2007)
 (23) Zöllner, F.G., Sance, R., Rogelj, P., LedesmaCarbayo, M.J., Rørvik, J., Santos, A., Lundervold, A.: Assessment of 3d dcemri of the kidneys using nonrigid image registration and segmentation of voxel time courses. Computerized Medical Imaging and Graphics 33(3), 171–181 (2009)
Comments
There are no comments yet.