• Previous Article
    Inverse problems for a half-order time-fractional diffusion equation in arbitrary dimension by Carleman estimates
  • IPI Home
  • This Issue
  • Next Article
    An inverse problem for a fractional diffusion equation with fractional power type nonlinearities
doi: 10.3934/ipi.2021059
Online First

Online First articles are published articles within a journal that have not yet been assigned to a formal issue. This means they do not yet have a volume number, issue number, or page numbers assigned to them, however, they can still be found and cited using their DOI (Digital Object Identifier). Online First publication benefits the research community by making new scientific discoveries known as quickly as possible.

Readers can access Online First articles via the “Online First” tab for the selected journal.

Joint reconstruction and low-rank decomposition for dynamic inverse problems

1. 

Department of Computer Science, University College London, London WC1E 6BT, United Kingdom

2. 

Center for Industrial Mathematics, University of Bremen, 28359 Bremen, Germany

3. 

Department of Mathematical Sciences, University of Oulu, 90014 Oulu, Finland

*Corresponding author: Pascal Fernsel

Received  June 2020 Revised  June 2021 Early access October 2021

A primary interest in dynamic inverse problems is to identify the underlying temporal behaviour of the system from outside measurements. In this work, we consider the case, where the target can be represented by a decomposition of spatial and temporal basis functions and hence can be efficiently represented by a low-rank decomposition. We then propose a joint reconstruction and low-rank decomposition method based on the Nonnegative Matrix Factorisation to obtain the unknown from highly undersampled dynamic measurement data. The proposed framework allows for flexible incorporation of separate regularisers for spatial and temporal features. For the special case of a stationary operator, we can effectively use the decomposition to reduce the computational complexity and obtain a substantial speed-up. The proposed methods are evaluated for three simulated phantoms and we compare the obtained results to a separate low-rank reconstruction and subsequent decomposition approach based on the widely used principal component analysis.

Citation: Simon Arridge, Pascal Fernsel, Andreas Hauptmann. Joint reconstruction and low-rank decomposition for dynamic inverse problems. Inverse Problems & Imaging, doi: 10.3934/ipi.2021059
References:
[1]

S. Arridge, P. Fernsel and A. Hauptmann, Joint reconstruction and low-rank decomposition for dynamic, Available online on GitLab: Inverse Problems - Support Code and Reconstruction Videos, 2021. https://gitlab.informatik.uni-bremen.de/s_p32gf3/joint_reconstruction_lowrank_decomp_dynamicip. Google Scholar

[2]

D. Böhning and B. G. Lindsay, Monotonicity of quadratic-approximation algorithms, Ann. Inst. Statist. Math., 40 (1988), 641-663.  doi: 10.1007/BF00049423.  Google Scholar

[3]

C. Boutsidis and E. Gallopoulos, SVD based initialization: A head start for nonnegative matrix factorization, Pattern Recognition, 41 (2008), 1350-1362.  doi: 10.1016/j.patcog.2007.09.010.  Google Scholar

[4]

D. BrunetE. R. Vrscay and Z. Wang, On the mathematical properties of the structural similarity index, IEEE Trans. Image Process., 21 (2012), 1488-1499.  doi: 10.1109/TIP.2011.2173206.  Google Scholar

[5]

T. A. Bubba, M. März, Z. Purisha, M. Lassas and S. Siltanen, Shearlet-based regularization in sparse dynamic tomography, In Wavelets and Sparsity XVII, International Society for Optics and Photonics, 10394 (2017), 236–245. doi: 10.1117/12.2273380.  Google Scholar

[6]

M. Burger, H. Dirks, L. Frerking, A. Hauptmann, T. Helin and S. Siltanen, A variational reconstruction method for undersampled dynamic X-ray tomography based on physical motion models, Inverse Problems, 33 (2017), 24pp. doi: 10.1088/1361-6420/aa99cf.  Google Scholar

[7]

M. BurgerH. Dirks and C.-B. Schönlieb, A variational model for joint motion estimation and image reconstruction, SIAM J. Imaging Sci., 11 (2018), 94-128.  doi: 10.1137/16M1084183.  Google Scholar

[8]

J. CaiX. JiaH. GaoS. B. JiangZ. Shen and H. Zhao, Cine cone beam ct reconstruction using low-rank matrix factorization: Algorithm and a proof-of-principle study, IEEE Transactions on Medical Imaging, 33 (2014), 1581-1591.  doi: 10.1109/TMI.2014.2319055.  Google Scholar

[9]

E. J. CandèsX. LiY. Ma and J. Wright, Robust principal component analysis?, J. ACM, 58 (2011), 1-37.  doi: 10.1145/1970392.1970395.  Google Scholar

[10]

B. ChenJ. Abascal and M. Soleimani, Extended joint sparsity reconstruction for spatial and temporal ERT imaging, Sensors, 18 (2018), 4014.  doi: 10.3390/s18114014.  Google Scholar

[11]

C. Chen and O. Öktem, Indirect image registration with large diffeomorphic deformations, SIAM J. Imaging Sci., 11 (2018), 575-617.  doi: 10.1137/17M1134627.  Google Scholar

[12]

A. Cichocki, R. Zdunek, A. H. Phan and S. Amari, Nonnegative Matrix and Tensor Factorizations: Applications to Exploratory Multi–Way Data Analysis and Blind Source Separation, Wiley Publishing, 2009. Google Scholar

[13]

C. De Mol, Blind Deconvolution and Nonnegative Matrix Factorization, Oberwolfach Reports 51/2012, Mathematisches Forschungsinstitut Oberwolfach, 2012. Google Scholar

[14]

M. Defrise, C. Vanhove and X. Liu, An algorithm for total variation regularization in high-dimensional linear problems, Inverse Problems, 27 (2011), 16pp. doi: 10.1088/0266-5611/27/6/065002.  Google Scholar

[15]

C. Ding, X. He and H. D. Simon, On the equivalence of nonnegative matrix factorization and spectral clustering, In Proceedings of the 2005 SIAM International Conference on Data Mining, 5 2005,606–610. doi: 10.1137/1.9781611972757.70.  Google Scholar

[16]

N. DjurabekovaA. GoldbergA. HauptmannD. HawkesG. LongF. Lucka and M. Betcke, Application of proximal alternating linearized minimization (PALM) and inertial PALM to dynamic 3D CT, 15th International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine, 11072 (2019), 30-34.  doi: 10.1117/12.2534827.  Google Scholar

[17]

D. Driggs, J. Tang, J. Liang, M. Davies and C.-B. Schönlieb, Spring: A fast stochastic proximal alternating method for non-smooth non-convex optimization, preprint, arXiv: 2002.12266. Google Scholar

[18]

L. FengR. GrimmK. T. BlockH. ChandaranaS. KimJ. XuL. AxelD. K. Sodickson and R. Otazo, Golden-angle radial sparse parallel MRI: Combination of compressed sensing, parallel imaging, and golden-angle radial sampling for fast and flexible dynamic volumetric MRI, Magnetic Resonance in Medicine, 72 (2014), 707-717.  doi: 10.1002/mrm.24980.  Google Scholar

[19]

P. Fernsel and P. Maass, A survey on surrogate approaches to non-negative matrix factorization, Vietnam J. Math., 46 (2018), 987-1021.  doi: 10.1007/s10013-018-0315-x.  Google Scholar

[20]

C. FévotteN. Bertin and J.-L. Durrieu, Nonnegative matrix factorization with the itakura-saito-divergence: With application to music analysis, Neural Computation, 21 (2009), 793-830.   Google Scholar

[21]

H. GaoJ. CaiZ. Shen and H. Zhao, Robust principal component analysis-based four-dimensional computed tomography, Phys. Med. Biol., 56 (2011), 3181-3198.  doi: 10.1088/0031-9155/56/11/002.  Google Scholar

[22]

H. GaoY. ZhangL. Ren and F.-F. Yin, Principal component reconstruction (PCR) for cine CBCT with motion learning from 2d fluoroscopy, Medical Physics, 45 (2018), 167-177.  doi: 10.1002/mp.12671.  Google Scholar

[23]

T. Goldstein and S. Osher, The split Bregman method for l1-regularized problems, SIAM J. Imaging Sci., 2 (2009), 323-343.  doi: 10.1137/080725891.  Google Scholar

[24] G. H. Golub and C. F. Van Loan, Matrix Computations, 4$^nd$ edition, Johns Hopkins Studies in the Mathematical Sciences. Johns Hopkins University Press, Baltimore, MD, 2013.   Google Scholar
[25]

B. Gris, C. Chen and O. Öktem, Image reconstruction through metamorphosis, Inverse Problems, 36 (2020), 27pp. doi: 10.1088/1361-6420/ab5832.  Google Scholar

[26]

J. HakkarainenZ. PurishaA. Solonen and S. Siltanen, Undersampled dynamic X-ray tomography with dimension reduction Kalman filter, IEEE Transactions on Computational Imaging, 5 (2019), 492-501.  doi: 10.1109/TCI.2019.2896527.  Google Scholar

[27]

A. Hauptmann, O. Öktem and C. Schönlieb, Image reconstruction in dynamic inverse problems with temporal models, preprint, arXiv: 2007.10238. Google Scholar

[28]

K. S. Kim and J. C. Ye, Low-dose limited view 4d ct reconstruction using patch-based low-rank regularization, In 2013 IEEE Nuclear Science Symposium and Medical Imaging Conference (2013 NSS/MIC), (2013), 1–4. doi: 10.1109/NSSMIC.2013.6829178.  Google Scholar

[29]

B. KlingenbergJ. Curry and A. Dougherty, Non-negative matrix factorization: Ill-posedness and a geometric algorithm, Pattern Recognition, 42 (2009), 918-928.  doi: 10.1016/j.patcog.2008.08.026.  Google Scholar

[30]

D. Kressner and A. Uschmajew, On low-rank approximability of solutions to high-dimensional operator equations and eigenvalue problems, Linear Algebra Appl., 493 (2016), 556-572.  doi: 10.1016/j.laa.2015.12.016.  Google Scholar

[31]

K. Lange, Optimization, Springer Texts in Statistics, Springer-Verlag, New York, 2004. doi: 10.1007/978-1-4757-4182-7.  Google Scholar

[32]

K. LangeD. R. Hunter and I. Yang, Optimization transfer using surrogate objective functions, J. Comput. Graph. Statist., 9 (2000), 1-59.  doi: 10.2307/1390605.  Google Scholar

[33]

L. Lecharlier and C. De Mol, Regularized blind deconvolution with poisson data, Journal of Physics: Conference Series, 464 (2013), 012003.  doi: 10.1088/1742-6596/464/1/012003.  Google Scholar

[34]

D. D. Lee and H. S. Seung, Learning the parts of objects by non-negative matrix factorization, Nature, 401 (1999), 788-791.  doi: 10.1038/44565.  Google Scholar

[35]

D. D. Lee and H. S. Seung, Algorithms for non-negative matrix factorization, In Advances in Neural Information Processing Systems 13 - NIPS 2000, (2001), 535–541. Google Scholar

[36]

J. LeuschnerM. SchmidtP. FernselD. LachmundT. Boskamp and P. Maass, Supervised non-negative matrix factorization methods for maldi imaging applications, Bioinformatics, 35 (2019), 1940-1947.  doi: 10.1093/bioinformatics/bty909.  Google Scholar

[37]

S. G. LingalaY. HuE. V. R. DiBella and M. Jacob, Accelerated dynamic MRI exploiting sparsity and low-rank structure: K-t SLR, IEEE Transactions on Medical Imaging, 30 (2011), 1042-1054.  doi: 10.1109/TMI.2010.2100850.  Google Scholar

[38]

F. LuckaN. HuynhM. BetckeE. ZhangP. BeardB. Cox and S. Arridge, Enhancing compressed sensing 4D photoacoustic tomography by simultaneous motion estimation, SIAM J. Imaging Sci., 11 (2018), 2224-2253.  doi: 10.1137/18M1170066.  Google Scholar

[39]

M. Lustig, J. M. Santos, D. L. Donoho and J. M. Pauly, kt SPARSE: High frame rate dynamic MRI exploiting spatio-temporal sparsity, In Proceedings of the 13th annual meeting of ISMRM, Seattle, 2420 (2006). Google Scholar

[40]

E. NiemiM. LassasA. KallonenL. HarhanenK. Hämäläinen and S. Siltanen, Dynamic multi-source X-ray tomography using a spacetime level set method, J. Comput. Phys., 291 (2015), 218-237.  doi: 10.1016/j.jcp.2015.03.016.  Google Scholar

[41]

J. P. OliveiraJ. M. Bioucas-Dias and M. A. T. Figueiredo, Review: Adaptive total variation image deblurring: A majorization-minimization approach, Signal Processing, 89 (2009), 1683-1693.  doi: 10.1016/j.sigpro.2009.03.018.  Google Scholar

[42]

P. Paatero and U. Tapper, Positive matrix factorization: A non-negative factor model with optimal utilization of error estimates of data values, Environmetrics, 5 (1994), 111-126.  doi: 10.1002/env.3170050203.  Google Scholar

[43]

K. Pearson, On lines and planes of closest fit to systems of points in space, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 2 (1901), 559-572.  doi: 10.1080/14786440109462720.  Google Scholar

[44]

U. Schmitt and A. K. Louis, Efficient algorithms for the regularization of dynamic inverse problems: Ⅰ. Theory, Inverse Problems, 18 (2002), 645-658.  doi: 10.1088/0266-5611/18/3/308.  Google Scholar

[45]

U. SchmittA. K. LouisC. Wolters and M. Vauhkonen, Efficient algorithms for the regularization of dynamic inverse problems: Ⅱ. Applications, Inverse Problems, 18 (2002), 659-676.  doi: 10.1088/0266-5611/18/3/309.  Google Scholar

[46]

J. A. SteedenG. T. KowalikO. TannM. HughesK. H. Mortensen and V. Muthurangu, Real-time assessment of right and left ventricular volumes and function in children using high spatiotemporal resolution spiral bSSFP with compressed sensing, Journal of Cardiovascular Magnetic Resonance, 20 (2018), 79.  doi: 10.1186/s12968-018-0500-9.  Google Scholar

[47]

M. Tao and X. Yuan, Recovering low-rank and sparse components of matrices from incomplete and noisy observations, SIAM J. Optim., 21 (2011), 57-81.  doi: 10.1137/100781894.  Google Scholar

[48]

B. R. Trémoulhéac, Low-rank and Sparse Reconstruction in Dynamic Magnetic Resonance Imaging Via Proximal Splitting Methods, PhD thesis, University College London, 2014. Google Scholar

[49]

B. R. TrémoulhéacN. DikaiosD. Atkinson and S. Arridge, Dynamic MR image reconstruction-separation from undersampled (k, t)-space via low-rank plus sparse prior, IEEE Transactions on Medical Imaging, 33 (2014), 1689-1701.  doi: 10.1109/TMI.2014.2321190.  Google Scholar

[50]

A. Uschmajew, On Low-Rank Approximation in Tensor Product Hilbert Spaces, PhD thesis, Technische Universität Berlin, 2013. Google Scholar

[51]

S. WundrakJ. PaulJ. UlriciE. Hell and V. Rasche, A small surrogate for the golden angle in time-resolved radial MRI based on generalized Fibonacci sequences, IEEE Transactions on Medical Imaging, 34 (2014), 1262-1269.  doi: 10.1109/TMI.2014.2382572.  Google Scholar

[52]

X. YuS. ChenZ. HuM. LiuY. ChenP. Shi and H. Liu, Sparse/low rank constrained reconstruction for dynamic pet imaging, PLOS ONE, 10 (2015), 1-18.  doi: 10.1371/journal.pone.0142019.  Google Scholar

[53]

X. M. Yuan and J. F. Yang, Sparse and low-rank matrix decomposition via alternating direction methods, Pac. J. Optim., 9 (2013), 167-180.   Google Scholar

show all references

References:
[1]

S. Arridge, P. Fernsel and A. Hauptmann, Joint reconstruction and low-rank decomposition for dynamic, Available online on GitLab: Inverse Problems - Support Code and Reconstruction Videos, 2021. https://gitlab.informatik.uni-bremen.de/s_p32gf3/joint_reconstruction_lowrank_decomp_dynamicip. Google Scholar

[2]

D. Böhning and B. G. Lindsay, Monotonicity of quadratic-approximation algorithms, Ann. Inst. Statist. Math., 40 (1988), 641-663.  doi: 10.1007/BF00049423.  Google Scholar

[3]

C. Boutsidis and E. Gallopoulos, SVD based initialization: A head start for nonnegative matrix factorization, Pattern Recognition, 41 (2008), 1350-1362.  doi: 10.1016/j.patcog.2007.09.010.  Google Scholar

[4]

D. BrunetE. R. Vrscay and Z. Wang, On the mathematical properties of the structural similarity index, IEEE Trans. Image Process., 21 (2012), 1488-1499.  doi: 10.1109/TIP.2011.2173206.  Google Scholar

[5]

T. A. Bubba, M. März, Z. Purisha, M. Lassas and S. Siltanen, Shearlet-based regularization in sparse dynamic tomography, In Wavelets and Sparsity XVII, International Society for Optics and Photonics, 10394 (2017), 236–245. doi: 10.1117/12.2273380.  Google Scholar

[6]

M. Burger, H. Dirks, L. Frerking, A. Hauptmann, T. Helin and S. Siltanen, A variational reconstruction method for undersampled dynamic X-ray tomography based on physical motion models, Inverse Problems, 33 (2017), 24pp. doi: 10.1088/1361-6420/aa99cf.  Google Scholar

[7]

M. BurgerH. Dirks and C.-B. Schönlieb, A variational model for joint motion estimation and image reconstruction, SIAM J. Imaging Sci., 11 (2018), 94-128.  doi: 10.1137/16M1084183.  Google Scholar

[8]

J. CaiX. JiaH. GaoS. B. JiangZ. Shen and H. Zhao, Cine cone beam ct reconstruction using low-rank matrix factorization: Algorithm and a proof-of-principle study, IEEE Transactions on Medical Imaging, 33 (2014), 1581-1591.  doi: 10.1109/TMI.2014.2319055.  Google Scholar

[9]

E. J. CandèsX. LiY. Ma and J. Wright, Robust principal component analysis?, J. ACM, 58 (2011), 1-37.  doi: 10.1145/1970392.1970395.  Google Scholar

[10]

B. ChenJ. Abascal and M. Soleimani, Extended joint sparsity reconstruction for spatial and temporal ERT imaging, Sensors, 18 (2018), 4014.  doi: 10.3390/s18114014.  Google Scholar

[11]

C. Chen and O. Öktem, Indirect image registration with large diffeomorphic deformations, SIAM J. Imaging Sci., 11 (2018), 575-617.  doi: 10.1137/17M1134627.  Google Scholar

[12]

A. Cichocki, R. Zdunek, A. H. Phan and S. Amari, Nonnegative Matrix and Tensor Factorizations: Applications to Exploratory Multi–Way Data Analysis and Blind Source Separation, Wiley Publishing, 2009. Google Scholar

[13]

C. De Mol, Blind Deconvolution and Nonnegative Matrix Factorization, Oberwolfach Reports 51/2012, Mathematisches Forschungsinstitut Oberwolfach, 2012. Google Scholar

[14]

M. Defrise, C. Vanhove and X. Liu, An algorithm for total variation regularization in high-dimensional linear problems, Inverse Problems, 27 (2011), 16pp. doi: 10.1088/0266-5611/27/6/065002.  Google Scholar

[15]

C. Ding, X. He and H. D. Simon, On the equivalence of nonnegative matrix factorization and spectral clustering, In Proceedings of the 2005 SIAM International Conference on Data Mining, 5 2005,606–610. doi: 10.1137/1.9781611972757.70.  Google Scholar

[16]

N. DjurabekovaA. GoldbergA. HauptmannD. HawkesG. LongF. Lucka and M. Betcke, Application of proximal alternating linearized minimization (PALM) and inertial PALM to dynamic 3D CT, 15th International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine, 11072 (2019), 30-34.  doi: 10.1117/12.2534827.  Google Scholar

[17]

D. Driggs, J. Tang, J. Liang, M. Davies and C.-B. Schönlieb, Spring: A fast stochastic proximal alternating method for non-smooth non-convex optimization, preprint, arXiv: 2002.12266. Google Scholar

[18]

L. FengR. GrimmK. T. BlockH. ChandaranaS. KimJ. XuL. AxelD. K. Sodickson and R. Otazo, Golden-angle radial sparse parallel MRI: Combination of compressed sensing, parallel imaging, and golden-angle radial sampling for fast and flexible dynamic volumetric MRI, Magnetic Resonance in Medicine, 72 (2014), 707-717.  doi: 10.1002/mrm.24980.  Google Scholar

[19]

P. Fernsel and P. Maass, A survey on surrogate approaches to non-negative matrix factorization, Vietnam J. Math., 46 (2018), 987-1021.  doi: 10.1007/s10013-018-0315-x.  Google Scholar

[20]

C. FévotteN. Bertin and J.-L. Durrieu, Nonnegative matrix factorization with the itakura-saito-divergence: With application to music analysis, Neural Computation, 21 (2009), 793-830.   Google Scholar

[21]

H. GaoJ. CaiZ. Shen and H. Zhao, Robust principal component analysis-based four-dimensional computed tomography, Phys. Med. Biol., 56 (2011), 3181-3198.  doi: 10.1088/0031-9155/56/11/002.  Google Scholar

[22]

H. GaoY. ZhangL. Ren and F.-F. Yin, Principal component reconstruction (PCR) for cine CBCT with motion learning from 2d fluoroscopy, Medical Physics, 45 (2018), 167-177.  doi: 10.1002/mp.12671.  Google Scholar

[23]

T. Goldstein and S. Osher, The split Bregman method for l1-regularized problems, SIAM J. Imaging Sci., 2 (2009), 323-343.  doi: 10.1137/080725891.  Google Scholar

[24] G. H. Golub and C. F. Van Loan, Matrix Computations, 4$^nd$ edition, Johns Hopkins Studies in the Mathematical Sciences. Johns Hopkins University Press, Baltimore, MD, 2013.   Google Scholar
[25]

B. Gris, C. Chen and O. Öktem, Image reconstruction through metamorphosis, Inverse Problems, 36 (2020), 27pp. doi: 10.1088/1361-6420/ab5832.  Google Scholar

[26]

J. HakkarainenZ. PurishaA. Solonen and S. Siltanen, Undersampled dynamic X-ray tomography with dimension reduction Kalman filter, IEEE Transactions on Computational Imaging, 5 (2019), 492-501.  doi: 10.1109/TCI.2019.2896527.  Google Scholar

[27]

A. Hauptmann, O. Öktem and C. Schönlieb, Image reconstruction in dynamic inverse problems with temporal models, preprint, arXiv: 2007.10238. Google Scholar

[28]

K. S. Kim and J. C. Ye, Low-dose limited view 4d ct reconstruction using patch-based low-rank regularization, In 2013 IEEE Nuclear Science Symposium and Medical Imaging Conference (2013 NSS/MIC), (2013), 1–4. doi: 10.1109/NSSMIC.2013.6829178.  Google Scholar

[29]

B. KlingenbergJ. Curry and A. Dougherty, Non-negative matrix factorization: Ill-posedness and a geometric algorithm, Pattern Recognition, 42 (2009), 918-928.  doi: 10.1016/j.patcog.2008.08.026.  Google Scholar

[30]

D. Kressner and A. Uschmajew, On low-rank approximability of solutions to high-dimensional operator equations and eigenvalue problems, Linear Algebra Appl., 493 (2016), 556-572.  doi: 10.1016/j.laa.2015.12.016.  Google Scholar

[31]

K. Lange, Optimization, Springer Texts in Statistics, Springer-Verlag, New York, 2004. doi: 10.1007/978-1-4757-4182-7.  Google Scholar

[32]

K. LangeD. R. Hunter and I. Yang, Optimization transfer using surrogate objective functions, J. Comput. Graph. Statist., 9 (2000), 1-59.  doi: 10.2307/1390605.  Google Scholar

[33]

L. Lecharlier and C. De Mol, Regularized blind deconvolution with poisson data, Journal of Physics: Conference Series, 464 (2013), 012003.  doi: 10.1088/1742-6596/464/1/012003.  Google Scholar

[34]

D. D. Lee and H. S. Seung, Learning the parts of objects by non-negative matrix factorization, Nature, 401 (1999), 788-791.  doi: 10.1038/44565.  Google Scholar

[35]

D. D. Lee and H. S. Seung, Algorithms for non-negative matrix factorization, In Advances in Neural Information Processing Systems 13 - NIPS 2000, (2001), 535–541. Google Scholar

[36]

J. LeuschnerM. SchmidtP. FernselD. LachmundT. Boskamp and P. Maass, Supervised non-negative matrix factorization methods for maldi imaging applications, Bioinformatics, 35 (2019), 1940-1947.  doi: 10.1093/bioinformatics/bty909.  Google Scholar

[37]

S. G. LingalaY. HuE. V. R. DiBella and M. Jacob, Accelerated dynamic MRI exploiting sparsity and low-rank structure: K-t SLR, IEEE Transactions on Medical Imaging, 30 (2011), 1042-1054.  doi: 10.1109/TMI.2010.2100850.  Google Scholar

[38]

F. LuckaN. HuynhM. BetckeE. ZhangP. BeardB. Cox and S. Arridge, Enhancing compressed sensing 4D photoacoustic tomography by simultaneous motion estimation, SIAM J. Imaging Sci., 11 (2018), 2224-2253.  doi: 10.1137/18M1170066.  Google Scholar

[39]

M. Lustig, J. M. Santos, D. L. Donoho and J. M. Pauly, kt SPARSE: High frame rate dynamic MRI exploiting spatio-temporal sparsity, In Proceedings of the 13th annual meeting of ISMRM, Seattle, 2420 (2006). Google Scholar

[40]

E. NiemiM. LassasA. KallonenL. HarhanenK. Hämäläinen and S. Siltanen, Dynamic multi-source X-ray tomography using a spacetime level set method, J. Comput. Phys., 291 (2015), 218-237.  doi: 10.1016/j.jcp.2015.03.016.  Google Scholar

[41]

J. P. OliveiraJ. M. Bioucas-Dias and M. A. T. Figueiredo, Review: Adaptive total variation image deblurring: A majorization-minimization approach, Signal Processing, 89 (2009), 1683-1693.  doi: 10.1016/j.sigpro.2009.03.018.  Google Scholar

[42]

P. Paatero and U. Tapper, Positive matrix factorization: A non-negative factor model with optimal utilization of error estimates of data values, Environmetrics, 5 (1994), 111-126.  doi: 10.1002/env.3170050203.  Google Scholar

[43]

K. Pearson, On lines and planes of closest fit to systems of points in space, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 2 (1901), 559-572.  doi: 10.1080/14786440109462720.  Google Scholar

[44]

U. Schmitt and A. K. Louis, Efficient algorithms for the regularization of dynamic inverse problems: Ⅰ. Theory, Inverse Problems, 18 (2002), 645-658.  doi: 10.1088/0266-5611/18/3/308.  Google Scholar

[45]

U. SchmittA. K. LouisC. Wolters and M. Vauhkonen, Efficient algorithms for the regularization of dynamic inverse problems: Ⅱ. Applications, Inverse Problems, 18 (2002), 659-676.  doi: 10.1088/0266-5611/18/3/309.  Google Scholar

[46]

J. A. SteedenG. T. KowalikO. TannM. HughesK. H. Mortensen and V. Muthurangu, Real-time assessment of right and left ventricular volumes and function in children using high spatiotemporal resolution spiral bSSFP with compressed sensing, Journal of Cardiovascular Magnetic Resonance, 20 (2018), 79.  doi: 10.1186/s12968-018-0500-9.  Google Scholar

[47]

M. Tao and X. Yuan, Recovering low-rank and sparse components of matrices from incomplete and noisy observations, SIAM J. Optim., 21 (2011), 57-81.  doi: 10.1137/100781894.  Google Scholar

[48]

B. R. Trémoulhéac, Low-rank and Sparse Reconstruction in Dynamic Magnetic Resonance Imaging Via Proximal Splitting Methods, PhD thesis, University College London, 2014. Google Scholar

[49]

B. R. TrémoulhéacN. DikaiosD. Atkinson and S. Arridge, Dynamic MR image reconstruction-separation from undersampled (k, t)-space via low-rank plus sparse prior, IEEE Transactions on Medical Imaging, 33 (2014), 1689-1701.  doi: 10.1109/TMI.2014.2321190.  Google Scholar

[50]

A. Uschmajew, On Low-Rank Approximation in Tensor Product Hilbert Spaces, PhD thesis, Technische Universität Berlin, 2013. Google Scholar

[51]

S. WundrakJ. PaulJ. UlriciE. Hell and V. Rasche, A small surrogate for the golden angle in time-resolved radial MRI based on generalized Fibonacci sequences, IEEE Transactions on Medical Imaging, 34 (2014), 1262-1269.  doi: 10.1109/TMI.2014.2382572.  Google Scholar

[52]

X. YuS. ChenZ. HuM. LiuY. ChenP. Shi and H. Liu, Sparse/low rank constrained reconstruction for dynamic pet imaging, PLOS ONE, 10 (2015), 1-18.  doi: 10.1371/journal.pone.0142019.  Google Scholar

[53]

X. M. Yuan and J. F. Yang, Sparse and low-rank matrix decomposition via alternating direction methods, Pac. J. Optim., 9 (2013), 167-180.   Google Scholar

Figure 1.  Illustration of a phantom that can be represented by the decomposition in (2). The phantom consists of $ K = 3 $ components: the background and two dynamic components with periodically changing intensity (left and right plot). As such, this phantom can be efficiently represented by a low-rank decomposition considered in this study
Figure 2.  Structure of the NMF in the context of the dynamic Shepp-Logan phantom as shown in Figure 1. Here, the nonnegative spatial and temporal basis functions can be naturally represented by the matrices $ B $ and $ C $
Figure 3.  Results for the dynamic Shepp-Logan phantom considered in Section 3.1.1 with $ \vert \mathcal{I}_t \vert = 6 $ angles per time step and 1% Gaussian noise. Shown are the leading extracted features for the BC model and for BC-X
Figure 4.  Results for the dynamic Shepp-Logan phantom considered in Section 3.1.1 with $ \vert \mathcal{I}_t \vert = 6 $ angles per time step and 1% Gaussian noise. Shown are the leading extracted features for the $\texttt{gradTV_PCA}$ model and for $\texttt{gradTV_NMF}$
Figure 5.  Results for the dynamic Shepp-Logan phantom considered in Section 3.1.1 with $ \vert \mathcal{I}_t \vert = 6 $ angles per time step and 3% Gaussian noise. Shown are the leading extracted features for the BC model and for $\texttt{gradTV_PCA}$
Figure 6.  Results for the dynamic Shepp-Logan phantom considered in Section 3.1.1 with $ \vert \mathcal{I}_t \vert = 3 $ angles per time step and 1% Gaussian noise. Shown are the leading extracted features for the BC model and for BC-X
Figure 7.  Mean PSNR and SSIM values of the reconstructions of the dynamic Shepp-Logan phantom considered in Section 3.1.1 with 1% Gaussian noise for different numbers of projection angles
Figure 8.  Mean PSNR and SSIM values of the reconstructions of the dynamic Shepp-Logan phantom considered in Section 3.1.1 with 3% Gaussian noise for different numbers of projection angles
Figure 9.  Needed time in seconds for the reconstruction and feature extraction of the dynamic Shepp-Logan phantom considered in Section 3.1.1 with 1% Gaussian noise
Figure 10.  Results for the dynamic Shepp-Logan phantom considered in Section 3.1.1 with a stationary operator and 1% Gaussian noise. Shown are the leading extracted features with the sBC model for $ \vert \mathcal{I}_t \vert = 6 $ angles per time step and $ \vert \mathcal{I}_t \vert = 30. $
Figure 11.  Needed time in seconds, mean PSNR and mean SSIM values of the reconstructions of the dynamic Shepp-Logan phantom with 1% Gaussian noise for the stationary case sBC and different numbers of projection angles
Figure 12.  Illustration of the vessel phantom dataset consisting of $ T = 100 $ phantoms of dimension $ 264\times264, $ where the intensity of the blue highlighted area changes according to blue curve on the left
Figure 13.  Results for the vessel phantom with $ \vert \mathcal{I}_t \vert = 12 $ angles per time step and 1% Gaussian noise. Shown are the leading extracted features for the BC model and for BC-X
Figure 14.  Results for the vessel phantom with $ \vert \mathcal{I}_t \vert = 12 $ angles per time step and 1% Gaussian noise. Shown are the leading extracted features for the $\texttt{gradTV_PCA}$ model and for $\texttt{gradTV_NMF}$
Figure 15.  Results for the vessel phantom with $ \vert \mathcal{I}_t \vert = 12 $ angles per time step and 3% Gaussian noise. Shown are the leading extracted features for the BC model and for $\texttt{gradTV_PCA}$
Figure 16.  Selected time steps of the ground truth based on the dynamic Shepp-Logan phantom used in Section 3.1.3 for the numerical experiments. The dynamic parts of the dataset consist of the left ellipse, which periodically change its size, and the right ellipse, which periodically change its intensity
Figure 17.  Results for the dynamic Shepp-Logan phantom considered in Section 3.1.3 with $ \vert \mathcal{I}_t \vert = 12 $ angles per time step and 1% Gaussian noise. Shown are the extracted features for BC and BC-X with $ K = 4 $ and $ K = 3 $ respectively
Figure 18.  Results for the dynamic Shepp-Logan phantom considered in Section 3.1.3 with $ \vert \mathcal{I}_t \vert = 12 $ angles per time step and 1% Gaussian noise. Shown are the 5 leading extracted features for $\texttt{gradTV_PCA}$ and all features for $\texttt{gradTV_NMF}$
Figure 19.  Results for the dynamic Shepp-Logan phantom considered in Section 3.1.3 with $ \vert \mathcal{I}_t \vert = 12 $ angles per time step and 1% Gaussian noise. Shown are the extracted features for BC with $ K = 3 $ and $\texttt{gradTV_NMF}$ with $ K = 5. $
Figure 20.  Mean PSNR and SSIM values of the reconstructions of the dynamic Shepp-Logan phantom considered in Section 3.1.3 with 1% Gaussian noise for different numbers of projection angles
Figure 21.  Illustration of two iteration steps of the MM principle for a cost function $ \mathcal{F} $ with bounded curvature and a surrogate function $ \mathcal{Q}_\mathcal{F}, $ which is strictly convex in the first argument
Figure 22.  Plots of $ \Vert B_{\bullet, k} C_{k, \bullet} \Vert_{\infty} $ with $ K = 10 $ in descending order for the dynamic Shepp-Logan phantom (Figure 22a) and the vessel phantom (Figure 22b) considered in Section 3.1.1 and 3.1.2 with 1% Gaussian noise and the parameters given in Table 3 and 4. In the case of $\texttt{gradTV_PCA}$, the ten leading features with respect to the singular values are considered
Figure 23.  Plots of $ \Vert B_{\bullet, k} C_{k, \bullet} \Vert_{\infty} $ with $ K = 10 $ in descending order for the dynamic Shepp-Logan phantom considered in Section 3.1.3 with 1% Gaussian noise and the parameters given in Table 5. In the case of $\texttt{gradTV_PCA}$, the ten leading features with respect to the singular values are considered
Table0 
Algorithm 1 $\texttt{gradTV}$
1: Initialise: $ X $
2: Input: $ \rho_{\text{grad}}, \rho_{\text{thr}}, \rho_{\text{TV}} >0 $
3: repeat
4:   $ X_{\bullet, t} \gets X_{\bullet, t} - \rho_{\text{grad}} (A_t^\intercal A_tX_{\bullet, t} - A_t^\intercal Y_{\bullet, t}) \ \ \ \text{for all}\ t $
5:    $ (U,\Sigma, V) \gets {\text{SVD}}(X) $
6:   $ \Sigma \gets {\text{SOFTTHRESH}}_{\rho_{\text{thr}}}(\Sigma) $
7:   $ X \gets U\Sigma V^\intercal $
8:   $ X\gets \max(X, 0) $
9: until STOPPINGCRITERION satisfied
10: $ X \gets {\text{TVDENOISER}}_{\rho_{\text{TV}}}(X) $
11: return $ X $
Algorithm 1 $\texttt{gradTV}$
1: Initialise: $ X $
2: Input: $ \rho_{\text{grad}}, \rho_{\text{thr}}, \rho_{\text{TV}} >0 $
3: repeat
4:   $ X_{\bullet, t} \gets X_{\bullet, t} - \rho_{\text{grad}} (A_t^\intercal A_tX_{\bullet, t} - A_t^\intercal Y_{\bullet, t}) \ \ \ \text{for all}\ t $
5:    $ (U,\Sigma, V) \gets {\text{SVD}}(X) $
6:   $ \Sigma \gets {\text{SOFTTHRESH}}_{\rho_{\text{thr}}}(\Sigma) $
7:   $ X \gets U\Sigma V^\intercal $
8:   $ X\gets \max(X, 0) $
9: until STOPPINGCRITERION satisfied
10: $ X \gets {\text{TVDENOISER}}_{\rho_{\text{TV}}}(X) $
11: return $ X $
Table 1.  Summary and short explanation of the considered algorithms in the experimental section
Algorithm Description
BC Joint reconstruction and feature extraction with the NMF model BC without constructing $X$, see algorithm in Theorem 2.3
BC-X Joint reconstruction and feature extraction with NMF model BC-X and explicit construction of $X$, see algorithm in Theorem 2.2
sBC Joint reconstruction and feature extraction method with NMF model sBC for stationary operator, see algorithm in Corollary 1
$\texttt{gradTV}$ Low-rank based reconstruction method for $X$, see Algorithm 1
$\texttt{gradTV_PCA}$ Separated reconstruction and feature extraction with Algorithm 1 and subsequent PCA computation
$\texttt{gradTV_NMF}$ Separated reconstruction and feature extraction with Algorithm 1 and subsequent NMF computation based on the model in (6)
Algorithm Description
BC Joint reconstruction and feature extraction with the NMF model BC without constructing $X$, see algorithm in Theorem 2.3
BC-X Joint reconstruction and feature extraction with NMF model BC-X and explicit construction of $X$, see algorithm in Theorem 2.2
sBC Joint reconstruction and feature extraction method with NMF model sBC for stationary operator, see algorithm in Corollary 1
$\texttt{gradTV}$ Low-rank based reconstruction method for $X$, see Algorithm 1
$\texttt{gradTV_PCA}$ Separated reconstruction and feature extraction with Algorithm 1 and subsequent PCA computation
$\texttt{gradTV_NMF}$ Separated reconstruction and feature extraction with Algorithm 1 and subsequent NMF computation based on the model in (6)
Table 2.  Mean PSNR and SSIM values of the reconstruction results for the vessel phantom for different noise levels and numbers of projection angles. Values in brackets indicate that the dynamic part of the dataset in the corresponding experiment could not be reconstructed sufficiently well
BC BC-X $\texttt{gradTV}$
Noise $\vert \mathcal{I}_t \vert$ PSNR SSIM PSNR SSIM PSNR SSIM
1% 7 (34.130) (0.9016) 32.969 0.8382 31.463 0.8414
1% 12 35.050 0.9068 33.919 0.8496 34.309 0.8839
3% 12 30.148 0.7484 28.119 0.5708 29.375 0.6698
BC BC-X $\texttt{gradTV}$
Noise $\vert \mathcal{I}_t \vert$ PSNR SSIM PSNR SSIM PSNR SSIM
1% 7 (34.130) (0.9016) 32.969 0.8382 31.463 0.8414
1% 12 35.050 0.9068 33.919 0.8496 34.309 0.8839
3% 12 30.148 0.7484 28.119 0.5708 29.375 0.6698
Table 3.  Parameter choice of the experiments in Section 3.1.1 for the dynamic Shepp-Logan phantom for 1% and 3% Gaussian noise
BC BC-X $\texttt{gradTV}$
Parameter 1% noise 3% noise 1% noise 3% noise 1% noise 3% noise
$\alpha$ - - 70 70 - -
$\mu_C$ 0.1 0.1 0.1 0.1 - -
$\tau$ 10 50 6 20 - -
$\rho_{\text{grad}}$ - - - - $1\cdot 10^{-3}$ $8\cdot 10^{-4}$
$\rho_{\text{thr}}$ - - - - $7\cdot 10^{-4}$ $1\cdot 10^{-3}$
$\rho_{\text{TV}}$ - - - - $1\cdot 10^{-2}$ $2.5\cdot 10^{-2}$
$\tilde{\mu}_{C}$ - - - - $0.1$ $0.1$
BC BC-X $\texttt{gradTV}$
Parameter 1% noise 3% noise 1% noise 3% noise 1% noise 3% noise
$\alpha$ - - 70 70 - -
$\mu_C$ 0.1 0.1 0.1 0.1 - -
$\tau$ 10 50 6 20 - -
$\rho_{\text{grad}}$ - - - - $1\cdot 10^{-3}$ $8\cdot 10^{-4}$
$\rho_{\text{thr}}$ - - - - $7\cdot 10^{-4}$ $1\cdot 10^{-3}$
$\rho_{\text{TV}}$ - - - - $1\cdot 10^{-2}$ $2.5\cdot 10^{-2}$
$\tilde{\mu}_{C}$ - - - - $0.1$ $0.1$
Table 4.  Parameter choice of the experiments in Section 3.1.2 for the vessel phantom for 1% and 3% Gaussian noise
BC BC-X $\texttt{gradTV}$
Parameter 1% noise 3% noise 1% noise 3% noise 1% noise 3% noise
$\alpha$ - - $3\cdot 10^{2}$ $3\cdot 10^{2}$ - -
$\mu_C$ 1 1 1 1 - -
$\tau$ $1.3\cdot 10^{2}$ $4.3\cdot 10^{2}$ $90$ $3\cdot 10^{2}$ - -
$\rho_{\text{grad}}$ - - - - $2\cdot 10^{-4}$ $8\cdot 10^{-5}$
$\rho_{\text{thr}}$ - - - - $2\cdot 10^{-4}$ $2.5\cdot 10^{-4}$
$\rho_{\text{TV}}$ - - - - $2\cdot 10^{-2}$ $4\cdot 10^{-2}$
$\tilde{\mu}_{C}$ - - - - $0.1$ $0.1$
BC BC-X $\texttt{gradTV}$
Parameter 1% noise 3% noise 1% noise 3% noise 1% noise 3% noise
$\alpha$ - - $3\cdot 10^{2}$ $3\cdot 10^{2}$ - -
$\mu_C$ 1 1 1 1 - -
$\tau$ $1.3\cdot 10^{2}$ $4.3\cdot 10^{2}$ $90$ $3\cdot 10^{2}$ - -
$\rho_{\text{grad}}$ - - - - $2\cdot 10^{-4}$ $8\cdot 10^{-5}$
$\rho_{\text{thr}}$ - - - - $2\cdot 10^{-4}$ $2.5\cdot 10^{-4}$
$\rho_{\text{TV}}$ - - - - $2\cdot 10^{-2}$ $4\cdot 10^{-2}$
$\tilde{\mu}_{C}$ - - - - $0.1$ $0.1$
Table 5.  Parameter choice of the experiments in Section 3.1.3 for the dynamic Shepp-Logan phantom for 1% Gaussian noise
Parameter BC BC-X gradTV
$\alpha$ - $70$ -
$\mu_C$ $0.1$ $0.1$ -
$\tau$ $10$ $10$ -
$\rho_{\text{grad}}$ - - $1\cdot 10^{-3}$
$\rho_{\text{thr}}$ - - $2\cdot 10^{-4}$
$\rho_{\text{TV}}$ - - $1\cdot 10^{-2}$
$\tilde{\mu}_{C}$ - - $0.1$
Parameter BC BC-X gradTV
$\alpha$ - $70$ -
$\mu_C$ $0.1$ $0.1$ -
$\tau$ $10$ $10$ -
$\rho_{\text{grad}}$ - - $1\cdot 10^{-3}$
$\rho_{\text{thr}}$ - - $2\cdot 10^{-4}$
$\rho_{\text{TV}}$ - - $1\cdot 10^{-2}$
$\tilde{\mu}_{C}$ - - $0.1$
[1]

Tao Wu, Yu Lei, Jiao Shi, Maoguo Gong. An evolutionary multiobjective method for low-rank and sparse matrix decomposition. Big Data & Information Analytics, 2017, 2 (1) : 23-37. doi: 10.3934/bdia.2017006

[2]

Dan Zhu, Rosemary A. Renaut, Hongwei Li, Tianyou Liu. Fast non-convex low-rank matrix decomposition for separation of potential field data using minimal memory. Inverse Problems & Imaging, 2021, 15 (1) : 159-183. doi: 10.3934/ipi.2020076

[3]

Xianchao Xiu, Lingchen Kong. Rank-one and sparse matrix decomposition for dynamic MRI. Numerical Algebra, Control & Optimization, 2015, 5 (2) : 127-134. doi: 10.3934/naco.2015.5.127

[4]

Yangyang Xu, Ruru Hao, Wotao Yin, Zhixun Su. Parallel matrix factorization for low-rank tensor completion. Inverse Problems & Imaging, 2015, 9 (2) : 601-624. doi: 10.3934/ipi.2015.9.601

[5]

Yun Cai, Song Li. Convergence and stability of iteratively reweighted least squares for low-rank matrix recovery. Inverse Problems & Imaging, 2017, 11 (4) : 643-661. doi: 10.3934/ipi.2017030

[6]

Zhouchen Lin. A review on low-rank models in data analysis. Big Data & Information Analytics, 2016, 1 (2&3) : 139-161. doi: 10.3934/bdia.2016001

[7]

Huseyin Coskun. Nonlinear decomposition principle and fundamental matrix solutions for dynamic compartmental systems. Discrete & Continuous Dynamical Systems - B, 2019, 24 (12) : 6553-6605. doi: 10.3934/dcdsb.2019155

[8]

Daijun Jiang, Hui Feng, Jun Zou. Overlapping domain decomposition methods for linear inverse problems. Inverse Problems & Imaging, 2015, 9 (1) : 163-188. doi: 10.3934/ipi.2015.9.163

[9]

Simon Foucart, Richard G. Lynch. Recovering low-rank matrices from binary measurements. Inverse Problems & Imaging, 2019, 13 (4) : 703-720. doi: 10.3934/ipi.2019032

[10]

Junfeng Yang. Dynamic power price problem: An inverse variational inequality approach. Journal of Industrial & Management Optimization, 2008, 4 (4) : 673-684. doi: 10.3934/jimo.2008.4.673

[11]

Ke Wei, Jian-Feng Cai, Tony F. Chan, Shingyu Leung. Guarantees of riemannian optimization for low rank matrix completion. Inverse Problems & Imaging, 2020, 14 (2) : 233-265. doi: 10.3934/ipi.2020011

[12]

Jonathan H. Tu, Clarence W. Rowley, Dirk M. Luchtenburg, Steven L. Brunton, J. Nathan Kutz. On dynamic mode decomposition: Theory and applications. Journal of Computational Dynamics, 2014, 1 (2) : 391-421. doi: 10.3934/jcd.2014.1.391

[13]

Steven L. Brunton, Joshua L. Proctor, Jonathan H. Tu, J. Nathan Kutz. Compressed sensing and dynamic mode decomposition. Journal of Computational Dynamics, 2015, 2 (2) : 165-191. doi: 10.3934/jcd.2015002

[14]

Hao Zhang, Scott T. M. Dawson, Clarence W. Rowley, Eric A. Deem, Louis N. Cattafesta. Evaluating the accuracy of the dynamic mode decomposition. Journal of Computational Dynamics, 2020, 7 (1) : 35-56. doi: 10.3934/jcd.2020002

[15]

Simone Cacace, Maurizio Falcone. A dynamic domain decomposition for the eikonal-diffusion equation. Discrete & Continuous Dynamical Systems - S, 2016, 9 (1) : 109-123. doi: 10.3934/dcdss.2016.9.109

[16]

Alexandr Mikhaylov, Victor Mikhaylov. Dynamic inverse problem for Jacobi matrices. Inverse Problems & Imaging, 2019, 13 (3) : 431-447. doi: 10.3934/ipi.2019021

[17]

Michael Herty, Giuseppe Visconti. Kinetic methods for inverse problems. Kinetic & Related Models, 2019, 12 (5) : 1109-1130. doi: 10.3934/krm.2019042

[18]

Wei Wan, Weihong Guo, Jun Liu, Haiyang Huang. Non-local blind hyperspectral image super-resolution via 4d sparse tensor factorization and low-rank. Inverse Problems & Imaging, 2020, 14 (2) : 339-361. doi: 10.3934/ipi.2020015

[19]

Dongho Kim, Eun-Jae Park. Adaptive Crank-Nicolson methods with dynamic finite-element spaces for parabolic problems. Discrete & Continuous Dynamical Systems - B, 2008, 10 (4) : 873-886. doi: 10.3934/dcdsb.2008.10.873

[20]

Bernadette N. Hahn. Dynamic linear inverse problems with moderate movements of the object: Ill-posedness and regularization. Inverse Problems & Imaging, 2015, 9 (2) : 395-413. doi: 10.3934/ipi.2015.9.395

2020 Impact Factor: 1.639

Article outline

Figures and Tables

[Back to Top]