Sparse and Deep Learning Approaches for Biomedical Image Reconstruction

Summer School Lecture by Jong Chul Ye, 13th IEEE EMBS International Summer School on Biomedical Imaging, June 14-21, 2018 , Saint-Jacut de la Mer, France

Transcript

1. Jong Chul Ye Bio-Imaging & Signal Processing Lab. Dept. Bio & Brain Engineering Dept. Mathematical Sciences KAIST, Korea IEEE EMBS Summer School on Biomedical Imaging Sparse and Deep Learning Approaches for Biomedical Image Reconstruction: Part I: Compressed Sensing & Sparse Recovery This material can be downloaded from http://bispl.weebly.com
2. Bio Imaging Signal Processing Laboratory Prof. Jong Chul Ye jong.ye@kaist.ac.kr h#p://bisp.kaist.ac.kr *4.3.
   SFDPO
   $IBMMFOHF
   XJOOFS
   LU
   '0$644
   /P
   BMJBTJOH
   "DDFMFSBUFE
   'BTU
   JNBHJOH
   VTJOH
   $PNQSFTTFE
   TFOTJO H
   .3*
   "DDFMFSBUJPO
   /FVSP*NBHJOH
   TUBUJTUJDBM
   UPPM
   EFWFMPQNFOU
   .BDIJOF
   -FBSOJOH%FFQ
   -FBSOJOH
   GPS
   #JPNFEJDJOF
   /*3441.
   .3*
   
   4JHOBM
   1SPDFTTJOH
   .3
   BSUJGBDU
   SFNPWBM
   VTJOH
   "-0)"
   
   
   0QUJDBM
   *NBHJOH
   4VQFSSFTPMVUJPO
   JNBHJOH
   '"-$0/
   %JGGVTF
   PQUJDBM
   UPNPHSBQIZ
   
   %05
   YSBZ
   $5
   3FDPOTUSVDUJPO
   1&5
   3FDPOTUSVDUJPO
   $PSSFDUFE
   4DBUUFSFE
   %ZOBNJD
   %
   1&5
   SFDPOTUSVDUJPO
   The ima ge can not be The ima ge can not be The ima ge can not be The ima ge can not be (SPVOE
   USVUI
   TPVSDF
   PCKFDU
   '07
   EFUFDUPS
   #BDLQSPKFDUJPO
   SBOHF
   $IPSE
   MJOF
   Z
   Y
   *OUFSJPS
   5PNPHSBQIZ
   Bio-signal processing & Machine learning Applied mathematics, signal processing & machine learning for bio-medical imaging 
   ""1.
   -PX%PTF
   $5
   (SBOE
   $IBMMFOHF
   8JOOFS
   8BWF/FU
   4QBSTF41.
   
3. Course Overview 1.  Introduction 2.  Part I: Compressive Sensing •  Motion-compensated CS •  Learning-based CS •  Multiple measurement vector (MMV) CS 3.  Part II: Low-rank Hankel Matrix Approaches 4.  Part III: Deep learning approaches 5.  Open problems and outlook
4. Compressed Sensing Hankel Structured Matrix Comple7on Deep Learning 4 Roadmap: From CS to Deep Learning Coherent Theme of Sparse Recovery
5. Resolution Limits in Medical Imaging •  Diffraction Limit –  Resolution limit for optics, x-ray, PET, and etc. http://www.microscopyu.com/
6. Resolution Limits in Medical Imaging •  Nyquist Sampling limit –  Nyquist requires the minimun two times sample rate w.r.t. signal bandlimit –  Resolution limit for Fourier based imaging (MRI, etc.)
7. Law of Parsimony: Sparsity, Low-rank, LDM •  Occam’s Razor: law of parsimony by William Occam (14th century) Entities should not be multiplied beyond necessity. •  All things being equal, the simplest solution tends to be the best one.
8. Super-resolution Microscopy S. Hell et al, Science 2007.
9. Compressed Sensing MR Imaging •  Forward problem •  Sparse recovery problem Figure courtesy of Mathews Jacob CS
10. 9 View CT Reconstruction using Deep Learning
11. COMPRESSED SENSING
12. Compressed Sensing (CS) 12 measurements vector # non-zeros •  Incoherent projection •  Underdetermined system •  Sparse unknown vector Courtesy of Dr. Dror Baron b A x n ⇥ 1 k m ⇥ 1 m ' k log( n ) ⌧ n
13. Min-norm Solution is Not sparse
14. L0 Minimization •  Two fundamental questions –  Uniqueness ? –  Convex relaxation ? min x k x k0 subject to b = A x
15. Uniqueness and Spark
16. L1 Convex Relaxation min x k x k1 subject to b = A x k x k1 = n X i=1 | xi |
17. Intuition behind l1 R. Tibshirani, “Regression Shrinkage and Selection via the Lasso”, Jour. Royal Stat. Soc. B, 1996,
18. Sufficient Conditions for P1
19. RIP Conditions •  RIP has been developed for probabilistic argument .
20. RIP as Sufficient Condition
21. RIP: Robust Recovery
22. 22 Sparse Recovery Formulations min x kT x k1 subject to k b A x k  ✏ min ✓ k✓k1 subject to kb A ✓k  ✏ •  Analysis Formulation •  Synthesis Formulation
23. K. Scheffler et al, Eur Radiol (2003) 13:2409–2418 Sequence diagram of b-SSFP §  Demanding Acquisition Requirements An example of Cartesian TrueFISP parameter for cardiac imaging. The acquisition time is calculated based on 60bpm patients with 8 slice acquisition. Flip angle of 40-60° is commonly used. §  Balanced SSFP FLASH b-SSFP 23 / 38 Application : Cardiac MRI
24. Aliasing in accelerated Dynamic MRI t Down sampling along k-t space à Aliasing artifact
25. Classical Methods •  Parallel imaging (SENSE, GRAPPA) –  Prussmann et al, 1999, Griswold et al 2002 •  k-t space method –  UNFOLD –  UNFOLD-SMASH, TSENSE –  K-t BLAST/SENSE –  PARADISE, PARADIGM Blaimer et al, Top Magn Reson Imaging • V15, N 4, August 2004
26. How to Improve ? Sparsity !! t y KLT Cardiac MR fMRI
27. RIP in Dynamic MR k-t space sample Random phase encoding Temporal transform k t
28. k-t FOCUSS (Jung et al, PMB:2007, MRM:2009(a), MRM:2009(b)) k-t FOCUSS (our method) k-t BLAST/SENSE (J. Tsao, MRM, 2003) No-update Optimal from compressed sensing perspective Does not solve L1 minimization problem Not Optimal from CS perspective When there is no update, k-t FOCUSS is exactly same with k-t BLAST/SENSE
29. Zero-padding inverse FFT from measurements 11 x accelerated Sampling pattern
30. k-t BLAST 11 x accelerated Sampling pattern
31. k-t FOCUSS with temporal average 11 x accelerated Sampling pattern
32. MOTION-COMPENSATED COMPRESSED SENSING
33. Lessons from HDTV History MUSE MPEG The first HDTV (Japan) (obsolete) Modern Standard (MPEG1,2,4..) quincunx sampling ME/MC Lattice sampling theory à Low compression ME/MC + residual coding à High compression * MUSE: MUltiple Sub-nyquist Encoding systems Figures from Kovacevic et al, IEEE TIP, 1993
34. MUSE vs k-t MR MUSE k-t MR The first HDTV (Japan) (obsolete) UNFOLD k-t BLAST/SENSE quincunx sampling k-t space downsampling Lattice sampling theory à Low compression Lattice sampling theory à Low acceleration !! k time * MUSE: MUltiple Sub-nyquist Encoding systems Figures from Kovacevic et al, IEEE TIP, 1993
35. Reference is often Free •  Cardiac MR –  During volume acquisition, diastole phase can be acquired •  fMRI •  MR angiography, vocal tract, etc. •  example) RIGR, SPEAR, etc. 12 sec X 10 Finger Tapping 30 sec 30 sec 21 sec Rest Rest Dummy
36. Random Sampling with Reference When two reference frames can be obtained, k Reference frames Fully sampled center region t Random sampling When one reference frames can be obtained 0 x k t We want to distribute k-space samples to allow prediction and residual coding
37. ME/MC in MPEG Video t k Encoder Full sampled measurements t x y
38. ME/MC in MPEG Video t x y
39. ME/MC in MPEG Video t x y Residual encoding
40. ME/MC in Dynamic MRI t k Encoder Down sampled measurements t x y
41. Refinable Motion Estimation t k Encoder Down sampled measurements t x y
42. Refinable Motion Estimation t x y Initialization with k-t FOCUSS k t
43. Refinable Motion Estimation t x y Initialization with k-t FOCUSS k t
44. Refinable Motion Estimation t x y Apply ME/MC k t
45. Refinable Motion Estimation t x y Residual Encoding using k-t FOCUSS Residual Encoding ME/MC
46. Cine MRI •  Acquisition parameters - 1.5 T Philips scanner at Yonsei Unversity Medical Center in Korea - bSSFP - FOV = 345.00 x 270.00 mm2 - Matrix size : 256 x 220 ( Read-Out x PE ) - TR = 3.45ms - flip angle : 50° - heart rate : 66 bpm - # of cardiac phases : 25 frames •  Among fully sampled k-data, partial k-data were arbitrarily chosen.
47. Zero-padding inverse FFT from measurements 11 x accelerated Sampling pattern
48. k-t BLAST 11 x accelerated Sampling pattern
49. k-t FOCUSS 11 x accelerated Sampling pattern
50. Motion Compensated k-t FOCUSS 11 x accelerated Sampling pattern
51. Radial Acquisition
52. World-first In vivo Experiment (Joint work with K. Nayak at USC) •  Acquisition parameters - 3.0 T GE Signa EXCITE at the University of Southern California - bSSFP - FOV = 300 x 300 mm2 - Matrix size : 240 x 256 ( PE x Read-Out) - TR = 3.7ms -  views per segment = 10 - flip angle : 45° - heart rate : 65 bpm - # of ardiac phases : 25 frames -  acquisition time for 6 x accelerated data = 4 heart beats -  acquisition time for fully sampled data = 24 heart beats
53. In Vivo Experiment (6x accel.)
54. LEARNING-BASED COMPRESSED SENSING
55. Which One Is Optimal ? FT along t-axis PCA along t-axis WVT along x-axis + FT along t-axis x t x f
56. KLT (PCA): optimal transform Karhunen Loeve Transform is well known for an optimal sparsifying transform of general signals. KLT basis: Matrix of Eigen vectors for T X U V = Σ 0 0 L s s ⎛ ⎞ ⎜ ⎟ Σ = ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ T T U X V = Σ Considering SVD L: Large Singular values s: small Singular values Hence, KLT is SPARSE 0 L s >> → T XX X
57. PCA in k-t FOCUSS k-t space sample Random phase encoding Temporal transform t y SVD Significant eigenvector
58. 58 Two step Boosting for PCA basis
59. MR Angiography Golden Section Radial Trajetory kx ky kx ky kx ky kx ky … T=1 T=2 T=512 … ky kx ky kx ky kx … … Target time points: 39, 79, 118, 157, 197, 236, 275, 314, 354, 393, 432, 472 T=39 T=472 39 or 40 views are merged on each time frame Measurement size (512 read-out x 512 views) vs Total recon matrix (512 x 512 x 12 time frames) 12 x acceleration !!
60. KLT for CE-MRA x t k-t FOCUSS recon using FFT Summation along t-axis Thresholding
61. 2009 ISMRM Recon Challenge Award L2 minimization k-t FOCUSS using KLT
62. 62 T2 Molecular motion Fast (water) Slow (protein) (fat) Detect T2 changes due to variety of disease T2 Parameter Mapping (joint work D. Kim, NYU: Feng et al, MRM, 2011)
63. 63 Cardiac T2 Mapping 6 x accel. Conventional method 6 x accel. k-t FOCUSS
64. 64 Cardiac T2 Mapping 1.8 x accel. GRAPPA 6 x accel. k-t FOCUSS
65. 65 T2 Estimation
66. MULTIPLE MEASUREMENT VECTOR COMPRESSED SENSING
67. MMV problem - definition minimize kXk 0 subject to B = AX. k = kXk 0 : number of nonzero rows r : number of snapshots m : number of sensor elements Here, we assume that the columns of B are linearly independent. Jong Chul Ye (KAIST) The role of joint sparsity UIUC Talk 3 / 21 MMV: Multiple Measurement Vector Problems
68. Geophysical x x x x x x x x x x x x x x x x x x x Medical x x x x x x Industrial Electromagnetic Acoustic Ultrasonic Optical x x x x x x Electromagnetic Ultrasonic Optical Off-set VSP/ cross-well tomography GPR surface imaging induction imaging Ultrasound tomography optical microscopy photon imaging Ultrasound tomography optical microscopy induction imaging Inverse Scattering Problems Slide courtesy by Devaney
69. Diffuse Optical Tomography Near-infrared (NIR, ~650-950 nm) * Durduran et al., MICCAI, 2010 ⌦1 ⌦2
70. Elastic Source Imaging Elastography: Medical, geophysical applications https://iame.com/online/breast_elastography/ https://marcellusdrilling.com/2018/04/9-more-seismic- testing-devices-stolen-in-swpa-6-were-returned/
71. Joint Sparsity Multiple Measurement Vectors from Multiple illumination Patterns Finite # of snapshots are available
72. 72 Sparse Signal Perturbations in the optical parameters Angiogenesis in cancer Joint Sparsity
73. 73 Inverse Scattering Problems Ill-posed Non-linear Existing methods •  Born approximation à linearization error •  Iterative Born approximation à computationally expensive
74. Overcoming Nonlinearities using Joint Sparsity 74
75. l0 uniqueness results on MMV Definition Given a matrix A, let spark(A) denote the smallest number of linearly dependent columns of A. Theorem (Rao / Chen, Huo / Feng, Bresler / Davies, Eldar) If a matrix X satisfies AX = B, then kXk 0 < spark(A) + rank(B) 1 2 if and only if X is the unique solution to the l 0 minimization problem. With increasing number of snapshots, more non-zero elements can be recovered. Jong Chul Ye (KAIST) The role of joint sparsity UIUC Talk 5 / 21 L0 Uniqueness of MMV
76. Spark reduction 1 J. Kim, O. Lee, J. Ye, Compressive MUSIC (CS-MUSIC), IEEE Trans. on Inform. Theory, January 2012 Theorem Let r  m < n. Suppose that we are given a sensing matrix A 2 Rm⇥n and an observation matrix B 2 Rm⇥r . If the k nonzero rows of X are in general position (i.e., any collection of r nonzero rows are linearly independent) and A satisfies the RIP condition 0  L 2 k r +1 (A) < 1, then spark(Q⇤A) = k r + 1. Note that A 2 Rm⇥n satisfies RIP with 0  L 2 k r +1 (A) < 1 if and only if k < spark(A) + rank(B) 1 2 . :l 0 bound for MMV problem CS-MUSIC achieves l 0 bound with finite snapshot. Jong Chul Ye (KAIST) The role of joint sparsity UIUC Talk 6 / 21 Spark Reduction Principle
77. Generalized MUSIC criterion Theorem Assume that we have an MMV problem AX = B, where A, X and B as in the previous theorem. If Ik r ⇢ suppX with |Ik r | = k r and AIk r 2 Rm⇥ ( k r ) which consists of columns of A, whose indices are in Ik r , then for any j 2 {1, · · · , n} \ Ik r , a⇤ j h PR ( Q ) PR ( PR ( Q ) AIk r ) i aj = 0 () j 2 suppX. a⇤ j P? [ AIk r B ] aj = 0 () j 2 suppX. Augmented Signal Subspace: R[AIk r B] Jong Chul Ye (KAIST) The role of joint sparsity UIUC Talk 7 / 21 Generalized MUSIC Criterion
78. Geometric Interpretation Jong Chul Ye (KAIST) The role of joint sparsity UIUC Talk 8 / 21 Geometric Interpretation
79. Su cient condition for SS-OMP to find a partial support Case 1: r is a fixed number Theorem Let k ⇣ [ AIk r S ] ⌘ > 1 + k r . If we have m > 1 + 1 2k r ⌘k r 2k log (n k) r , then we can find k r correct indices of suppX by applying subspace S-OMP. 1. [Fletcher, Rangan] For SMV, when SNR ! 1, we need m > 2(1 + )k log (n k) for some > 0. 2. The sampling ratio is reciprocally proportional with respect to the number of multiple measurement vectors. 3. SNR! 1 is required, in the large system limit. Jong Chul Ye (KAIST) The role of joint sparsity UIUC Talk 9 / 21 (Similar to Reevs and Gaspar) Sampling Rate for MMV
80. Su cient condition for S-OMP to find a partial support Case 2: ↵ := limn!1 r(n)/k(n) > 0 Theorem Let k ⇣ [ AIk r S ] ⌘ > 1 + ↵ Then if we have m > (1 + )2 k 1 2⌘k r ↵ [2 F(↵)]2 , for some > 0 where F(↵) is an increasing function on [0, 1] such that F(1) = 1 and lim↵! 0+ F(↵) = 0. Then we can find k r correct indices of suppX by applying subspace S-OMP. 1. If ↵ ! 1 and SNR ! 1, then we only need to have m > (1 + )k for a small > 0. (cf. MUSIC) 2. In this case, the required number of sensor elements is at most 4(1 + )k when SNR! 1 (No log n factor). 3. The required SNR is finite. Jong Chul Ye (KAIST) The role of joint sparsity UIUC Talk 10 / 21 Sampling Rate for MMV
81. Other MMV Applications •  Parallel MRI + CS •  EEG/MEG •  Diffuse optical tomography •  Wave inverse scattering MMV for Bio-medical Imaging IFFT MR sca n
82. Optical Diffraction Tomography (ODT) Laser L1 L2 P BS1 GM L3 CL OL L4 L5 M1 L6 M2 BS2 Camera M3 P
83. ODT under Born/Rytov Approximation Sample z x(y) 3d k-space Scattered field 2 2-D Fourier transform and scaling
84. ODT using Joint Sparsity Lippman-Schwinger equation MMV for Support & Induced Current Recovery
85. ODT using Joint Sparsity Unknown Total Field Estimation Abnormality Estimation
86. Super-Resolution from Joint Sparse Recovery (Lim et al, OPEX, 2017)
87. Super-Resolution from Joint Sparse Recovery (Lim et al, OPEX, 2017)   = ; < ; < = = ; < ; < = D 5\WRY E 3URSRVHG
88. Super-Resolution from Joint Sparse Recovery (Lim et al, OPEX, 2017) = ; < ; < = = ; < ; < =   D 5\WRY E 3URSRVHG
89. 89 Diffuse Optical Tomography Applications (Lee, OPEX, 2013)
90. 90 Diffuse Optical Tomography Applications (Lee, OPEX, 2013)
91. Summary : Part I •  Sparsity principle is very important for biomedical image reconstruction︎ •  Compressed sensing has many extensions︎ •  Compressed sensing︎ •  Motion compensated CS︎ •  Learning-based CS︎ •  MMV CS︎ •  Downside: suffers from discretization, RIP, optimization bottleneck︎
92. References: Part I •  Ye, J. C., Tak, S., Han, Y., & Park, H. W. (2007). Projection reconstruction MR imaging using FOCUSS. Magnetic Resonance in Medicine, 57(4), 764-775. •  Jung, H., Ye, J. C., & Kim, E. Y. (2007). Improved k–t BLAST and k–t SENSE using FOCUSS.Physics in Medicine and Biology, 52(11), 320 •  Ye, J. C. (2007). Compressed sensing shape estimation of star-shaped objects in Fourier imaging. IEEE Signal Processing Letters, 14(10), 750-753. •  Jung, H., Sung, K., Nayak, K. S., Kim, E. Y., & Ye, J. C. (2009). K-t FOCUSS: A general compressed sensing framework for high resolution dynamic MRI. Magnetic Resonance in Medicine, 61(1), 103-116 •  Jin, K. H., Kim, Y., Yee, D. S., Lee, O. K., & Ye, J. C. (2009). Compressed sensing pulse-echo mode terahertz reflectance tomography. Optics Letters, 34(24), 3863-3865. •  Jung, H., & Ye, J. C. (2010). Motion estimated and compensated compressed sensing dynamic magnetic resonance imaging: What we can learn from video compression techniques. International Journal of Imaging Systems and Technology, 20(2), 81-98. •  Jung, H., Park, J., Yoo, J., & Ye, J. C. (2010). Radial k-t FOCUSS for high‐ resolution cardiac cine MRI. Magnetic Resonance in Medicine, 63(1), 68-7
93. References: Part I •  Choi, J., Kim, K. S., Kim, M. W., Seong, W., & Ye, J. C. (2011). Sparsity driven metal part reconstruction for artifact removal in dental CT. Journal of X-ray Science and Technology, 19(4), 457-475. •  Lee, O., Kim, J. M., Bresler, Y., & Ye, J. C. (2011). Compressive diffuse optical tomography: noniterative exact reconstruction using joint sparsity. IEEE Transactions on Medical Imaging, 30(5), 1129-1142. •  Lee, K., Tak, S., & Ye, J. C. (2011). A data-driven sparse GLM for fMRI analysis using sparse dictionary learning with MDL criterion. IEEE Transactions on Medical Imaging, 30(5), 1076-1089. •  Feng, L., Otazo, R., Jung, H., Jensen, J. H., Ye, J. C., Sodickson, D. K., & Kim, D. (2011). Accelerated cardiac T2 mapping using breath-hold multiecho fast spin- echo pulse sequence with kt FOCUSS. Magnetic Resonance in Medicine, 65(6), 1661-1669. •  Kim, J. M., Lee, O. K., & Ye, J. C. (2012). Compressive MUSIC: Revisiting the link between compressive sensing and array signal processing. IEEE Transactions on Information Theory, 58(1), 278-301.
94. References: Part I •  Kim, J. M., Lee, O. K., & Ye, J. C. (2012). Improving noise robustness in subspace-based joint sparse recovery. IEEE Transactions on Signal Processing, 60(11), 5799-5809. •  Min, J., Jang, J., Keum, D., Ryu, S. W., Choi, C., Jeong, K. H., & Ye, J. C. (2013). Fluorescent microscopy beyond diffraction limits using speckle illumination and joint support recovery. Scientific Reports, 3(2075) •  Lee, O., & Ye, J. C. (2013). Joint sparsity-driven non-iterative simultaneous reconstruction of absorption and scattering in diffuse optical tomography. Optics Express, 21(22), 26589-26604. •  Zong, X., Lee, J., Poplawsky, A. J., Kim, S. G., & Ye, J. C. (2014). Compressed sensing fMRI using gradient-recalled echo and EPI sequences. NeuroImage, 92, 312-321. •  J. Min C. Vonesch, H.Kirshner, L. Carlini, N. Olivier, S. Holden, S. Manley, J.C. Ye, M. Unser, (2014) FALCON: fast and unbiased reconstruction of high- density super-resolution microscopy data, Scientific Reports 4 , Article no 4577.
95. References: Part I •  Yoon, H., Kim, K. S., Kim, D., Bresler, Y., & Ye, J. C. (2014). Motion adaptive patch-based lowrank approach for compressed sensing cardiac cine MRI. IEEE Transactions on Medical Imaging, 33(11), 2069-2085. •  Min, J., Holden, S. J., Carlini, L., Unser, M., Manley, S., & Ye, J. C. (2014). 3D high-density localization microscopy using hybrid astigmatic/biplane imaging and sparse image reconstruction. Biomedical Optics Express, 5(11), 3935-3948. •  Han, P. K., Park, S. H., Kim, S. G., & Ye, J. C. (2015). Compressed sensing for fMRI: feasibility study on the acceleration of non-EPI fMRI at 9.4 T. BioMed research international, 2015. •  Ward, J. P., Lee, M., Ye, J. C., & Unser, M. (2015). Interior tomography using 1D generalized total variation. Part I: Mathematical foundation. SIAM Journal on Imaging Sciences, 8(1), 226-247 •  Lee, M., Han, Y., Ward, J. P., Unser, M., & Ye, J. C. (2015). Interior tomography using 1D generalized total variation. Part II: Multiscale implementation. SIAM Journal on Imaging Sciences, 8(4), 2452-2486.
96. References: Part I •  Kim, K., Ye, J. C., Worstell, W., Ouyang, J., Rakvongthai, Y., El Fakhri, G., & Li, Q. (2015). Sparse-view spectral CT reconstruction using spectral patch-based low-rank penalty. IEEE Transactions on Medical Imaging, 34(3), 748-760.. •  Kim, K., Son, Y. D., Bresler, Y., Cho, Z. H., Ra, J. B., & Ye, J. C. (2015). Dynamic PET reconstruction using temporal patch-based low rank penalty for ROI-based brain kinetic analysis. Physics in Medicine and Biology, 60(5), 2019. •  Lee, O. K., Kang, H., Ye, J. C., & Lim, M. (2015). A non-iterative method for the electrical impedance tomography based on joint sparse recovery. Inverse Problems, 31(7), 075002. •  Lim, J., Lee, K., Jin, K. H., Shin, S., Lee, S., Park, Y., & Ye, J. C. (2015). Comparative study of iterative reconstruction algorithms for missing cone problems in optical diffraction tomography. Optics Express, 23(13), 16933-16948. •  Lee, O., Tak, S., & Ye, J. C. (2015). A unified sparse recovery and inference framework for functional diffuse optical tomography using random effect model. IEEE Transactions on Medical Imaging, 34(7), 1602-1615.
97. References: Part I •  Ye, J. C., Kim, J. M., & Bresler, Y. (2015). Improving M-SBL for joint sparse recovery using a subspace penalty. IEEE Transactions on Signal Processing, 63(24), 6595-6605. •  Young-Beom Lee, Jeonghyeon Lee, Sungho Tak, Kangjoo Lee, Duk L. Na, Sangwon Seo, Yong Jeong, and Jong Chul Ye (2016). Sparse SPM: Group sparse-dictionary learning in SPM framework for resting-state functional connectivity MRI analysis, NeuroImage, 125, 1032–1045 •  Yoo, J., Jung, Y., Lim, M., Ye, J. C., & Wahab, A. (2017). A joint sparse recovery framework for accurate reconstruction of inclusions in elastic media. SIAM Journal on Imaging Sciences, 10(3), 1104-1138.
98. Jong Chul Ye Bio-Imaging & Signal Processing Lab. Dept. Bio & Brain Engineering Dept. Mathematical Sciences KAIST, Korea IEEE EMBS Summer School on Biomedical Imaging Sparse and Deep Learning Approaches for Biomedical Image Reconstruction: Part 2: Low-Rank Hankel Matrix Approaches
99. Compressed Sensing Hankel Structured Matrix Comple7on Deep Learning 99 Roadmap: From CS to Deep Learning Coherent Theme of Sparse Recovery
100. THEORY
101. Compressed Sensing (CS) 101 measurements vector # non-zeros •  Incoherent projection •  Underdetermined system •  Sparse unknown vector Courtesy of Dr. Dror Baron b A x n ⇥ 1 k m ⇥ 1 m ' k log( n ) ⌧ n
102. Limitations of CS •  Discrete domain theory •  Explicit form of sensing matrix •  RIP issue
103. Analytic Reconstruction (b) Time-reversal of a scattered wave (a) MR Imaging Beautiful analytic reconstruction results from fully sampled data
105. Finite Rate of Innovation Sampling Theory (Vetterlie, 2002) T -T 0 n1 -n1 0
106. FRI to Low-Rank Hankel matrix ︙ ︙ 1 2 3 4 5 6 7 8 9 -1 0 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 12 2 3 4 5 6 7 8 9 12 13 3 4 5 6 7 8 9 10 10 10 10 0 11 11 11 1 2 3 4 5 Finite length convolution Matrix Representation * ALOHA : Annihilating filter based LOw rank Hankel matrix Approach * Jin KH et al. IEEE TCI,2016 * Jin KH et al.,IEEE TIP, 2015 * Ye JC et al. IEEE TIT, 2016
107. Hankel Matrix
108. Low-Rank Hankel matrix minimization Missing elements can be found by low rank Hankel structured matrix compleAon Nuclear norm Projec;on on sampling posi;ons min m kH ( m ) k⇤ subject to P⌦(b) = P⌦( f ) RankH(f) = k * Jin KH et al IEEE TCI, 2016 * Jin KH et al.,IEEE TIP, 2015 * Ye JC et al., IEEE TIT, 2016 m
109. General TV Signals 109 Weighted Fourier data Piecewise smooth Splines, polynomials
110. 110 Existence of Annihilating Filter ˆ h(!) ⇤ ⇣ ˆ l(!) ˆ f(!) ⌘ = 0 Annihilating filter for weighted Fourier data General Low-Rank Hankel Matrix Completion
111. Visualization of concept
112. Extension to general signal models Stream of Diracs Stream of differentiated Diracs Non-uniform spline Piecewise smooth polynomial rank  X j dj rank = r rank  rq rank = r rank = X j dj rank = rq * Ye JC et al.,IEEE TIT 2016 With a proper weighAng, the Hankel matrix of the weighted k-space data à low ranked.
114. Performance Guarantees m c1µcsk log ↵ n min m kH(m)k⇤ subject to P⌦( m ) = P⌦(ˆ f ) min m kH(m)k⇤ subject to kP⌦( m ) P⌦(ˆ f ) k  kH(m) H(ˆ f)kF  c2n2 ↵ = ( 2 , on grid 4 , o↵ grid Exact Recovery Stable Recovery * Ye JC et al.,IEEE TIT 2016
115. Low-Rank Hankel matrix minimization Missing elements can be found by low rank Hankel structured matrix compleAon Nuclear norm Projec;on on sampling posi;ons min m kH ( m ) k⇤ subject to P⌦(b) = P⌦( f ) m
116. Optimization 116 •  Complexity depends on an estimated rank of patch !! ADMM
117. APPLICATIONS
118. Applications to Image Processing Inpainting & Impulse noise removal
119. Spectral Domain Sparsity 119 Smoothness, texture, pattern Sparse spectrum Smooth patch Edge patch Texture patch Concentrated at DC Elongated along Perpendicular to edge Distributed orthogonally w.r.t. texture
120. Low-rank Hankel matrix in Image
121. 2-D Hankel matrix 121 Hankel Matrix construction
122. 18. APR. 2015. 122 Why patch processing ? -  Spectrum changes for each patch -  Need to adapt the local Image statistics ALOHA
123. Rotation invariant sparsity 18. APR. 2015. 123 Hankel structured matrix is intrinsic low rank !! 0 deg. 90 deg. 0 90
124. Experimental results (x5) 18. APR. 2015. 124 *
125. Experimental results (x5) 18. APR. 2015. 125 *
126. 18. APR. 2015. 126 Text inlayed image reconstruction
127. 18. APR. 2015. 127 Line scratches
128. 18. APR. 2015. 128 Object removal
129. 129 *K.H. Jin, et. al, IEEE TIP (2015)
130. Lissajous scanning in MEMS lens scanning OCT microscopy 18. APR. 2015. 130 Park, Hyeon-Cheol, et al. "Forward imaging OCT endoscopic catheter based on MEMS lens scanning." Optics letters 37.13 (2012): 2673-2675. (*courtesy of H. Park and K. Jeong) * Lissajous scanning •  Lissajous scanning utilizes full speeds of both scanning axes! •  Lissajous scanner in MEMS controlled with only small voltage.
131. Temporal sampling in Lissajous scanning 18. APR. 2015. 131 Subsampled image in the regular grid can be interpolated by inpainting algorithms. Subsampled
132. Retrospective Experiments 18. APR. 2015. 132 * Data acquired by Lina Carlini (LEB, EPFL) *Endoplasmic Reticulum labeled with tdEOS fused to reticulon-4 in U2OS cells was taken on a modified epi-fluorescence microscope (Olympus IX71).
133. Real OCT experiment (x6) 18. APR. 2015. 133 * Data acquired by Hyeon-Cheol Park (Johns Hopkins Univ.) and Ki-Hun Jeong (KAIST) Korean coin 1.  140 x 140 px2 2.  Total data length : 150k
134. Magnified views (x6) 18. APR. 2015. 134
135. 135 Impulse Noise Removal (Jin et al, TIP, 2018)
136. 70% Salt and Pepper Noise Removal
137. 137 Multi-Channel Impulse Noise Removal
138. Recall: Sparsity in dynamic MRI t y KLT Cardiac MR fMRI
139. Sparsity of Dynamic MRI y (PE) Sparse In Wavelet domain x (RO) y (PE) x (RO) Sparsity #1 Sparsity #2 T2 images s-t spectrum y (PE) t y (PE) f Temporal Fourier tra nsform Hankel matrix with wavelet weigthing ky Hankel matrix on t 2D block Hankel matrix on ky -t t ky
140. k-t dynamic cases
141. Dynamic MRI – multi coil Six fold (x6) down sampling # of coils =4
142. Visualization of concept
143. Single coil static MRI
144. Parallel MRI Signal Model 144 i-th coil measurement Signals from parallel MRI i-th coil sensitivities Original signal × =
145. Multi-channel Generalized ALOHA in kx -ky 145 FT × × × = × Annihilation #2 : inter-coil relationships
146. Rank Bound for Parallel Imaging Jin et al, IEEE TCI, 2016 Sparsity of common image In transform domain Sparsity of sensitivity map In Fourier domain
147. 
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     DIBOOFM
     TJHOBM
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     DIBOOFM
     TJHOBM
     DIBOOFM
     TJHOBM
     E
     
     
     E
     
     
     E
     
     
     -PX
     3BOL
     .BUSJY
     
     $PNQMFUJPO
     -PX
     3BOL
     .BUSJY
     
     $PNQMFUJPO
     -PX
     3BOL
     .BUSJY
     
     $PNQMFUJPO
     -1 -1 -1 
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
148. Parallel MRI
149. What is MR parameter mapping? MR Parameter Mapping [1] Siemonsen, S et al., Radiology, 2008 [2] Rugg-Gunn, F. J., et al. NEUROLOGY, 2005 TE1 TE2 TE3 TE4 e.g. Multi-Echo Spin-Echo (ME-SE, T2 mapping) t Finding the quanAtaAve v alue of each ;ssue Cons Pros e.g. of T2 mapping for diagnosis for epilepsy Clinically valuable - As a quantitative diagnosis tool - Acute stroke, epilepsy, etc. Pros Cons TE1 TE2 TE3 TE4 e.g. Multi-echo images for T2 mapping Long scan time - Needs multiple scans - Variation of TI, TE, FA, etc.
150. Sparsity of Dynamic MRI y (PE) Sparse In Wavelet domain x (RO) y (PE) x (RO) Sparsity #1 Sparsity #2 T2 images s-t spectrum y (PE) t y (PE) f Temporal Fourier tra nsform Hankel matrix with wavelet weigthing ky Hankel matrix on t 2D block Hankel matrix on ky -t t ky
151. 47.5173×10-3
     52.9943×10-3
     47.6076×10-3
     47.9392×10-3
     17.6058×10-3
     60.6220×10-3
     ME-SE image Error ×5
     Reconstruction of x12.8 accelerated scan – ME-SE (4th echo) Result : in vivo acceleration study ( ME-SE, T2 ) k-t FOCUSS k-t SPARSE C-LORAKS Full k-t SLR ALOHA Patch Lee et al, MRM, 2015
152. 18.0169×10-2
     15.9958×10-2
     18.7474×10-2
     23.0070×10-2
     6.5596×10-2
     12.6578×10-2
     0 (ms) 50 100 150 T2 Map Error ×5
     Mapping from the reconstruction of x12.8 accelerated scan – T2 mapping Result : in vivo acceleration study ( ME-SE, T2 ) k-t FOCUSS k-t SPARSE C-LORAKS Full k-t SLR ALOHA Patch Lee et al, MRM, 2015
153. Result : Signal intensity curves ( SE-IR, T1 ) 400 800 1200 1600 0 Time (ms) -1 0 1 0SJHJOBM
     LU
     '0$644
     LU
     41"34&
     1BUDI
     LU
     4-3
     -03",4
     "-0)"
     Signal Intensity (a.u.) T1 relaxation
     LU
     '0$644
     LU
     41"34&
     1BUDI
     
     0SJHJOBM
     
     "-0)"
     -03",4
     LU
     4-3
     Lee et al, MRM, 2015
154. Summary of Results •  Goal : Acceleration of MR Parameter mapping by undersampling and reconstruction Full acquisition Scan Time Conventional scan Accelerated sc an 1min 2s 12min 50s Accelerated acq uisition Reconstructed i mages ALOHA Lee et al, MRM, 2015
155. MR artifacts •  What is MR artifacts? During acquisition, external interruptions (ex. fluctuation power supply of gradient, motion of object, etc.) distort signals.
     155 4QJLF
     OPJTF
     
     3FTQJSBUPSZ
     NPUJPO
     M. Graves, et. al., JMRI (2013)
156. Motivation 156
157. ü  )FSSJOHCPOF
     TQJLFT
     %
     LTQBDF
     Motivations: MR artifacts as sparse outliers 157
158. 
     ü  .PUJPO
     BSUJGBDU
     TQJLFT
     %
     LTQBDF
     QBSBMMFM
     UP
     SFBEPVU
     
     158 Motivations: MR artifacts as sparse outliers
159. 
     ;JQQFS
     BSUJGBDU
     TQJLFT
     %
     LTQBDF
     QFSQFOEJDVMBS
     UP
     SFBEPVU
     159 Motivations: MR artifacts as sparse outliers
160. Key Observation : Sparse outliers 160 ALOHA* Sparse MR artifact +
     4QBSTF
     PVUMJFS
     JT
     TUJMM
     TQBSTF
     JO
     XFJHIUFE
     )BOLFM
     NBUSJY
     † E. Candes, et. al, JACM (2011), R. Otazo, et. al, MRM (2015) *K.H. Jin, et. al, IEEE TIP (2015), K.H. Jin, et. al, arXiv (2015), J. C. Ye, et. al, arXiv (2015), J. Lee, et. al, MRM (2016), D. Lee, et. al, MRM (2016) •  ALOHA: Annihilating filter based LOw rank Hankel matrix Approach
161. Robust ALOHA 161 RPCA for weighted Hankel matrix ü  Extension of ALOHA for decomposition of sparse outliers (E) out of mixed signal* ü  Can be addressed ADMM† ü  K-space weighting signal Sparse outlier K-space weighting
     &
     $BOEFT
     FU
     BM
     +"$.
     
     3
     0UB[P
     FU
     BM
     .3.
     
     õ
     4
     #PZE
     FU
     BM
     'PVOEBUJPOT
     BOE
     5SFOET
     JO
     .BDIJOF
     -FBSOJOH
     
     ö
     ;
     8BO
     FU
     BM
     
     .BUIFNBUJDBM
     1SPHSBNNJOH
     $PNQVUBUJPO
     
     
162. 162 Algorithm Flowchart
163. Retrospective results )JHI
     JOUFOTJUZ
     4QJLF
     OPJTF
     -PX
     JOUFOTJUZ
     4QJLF
     OPJTF
     
     MPX
     GSFRVFODZ
     SFHJPO
     
     4QJLF
     OPJTF
     
     XJUI
     EPXO
     TBNQMJOH
     Y
     
164. Singular value distribution of Hankel matrix 164
165. In Vivo Motion artifact sudden motion (3 times)
166. Cardiac Motion artifact
167. Zipper artifact 167
168. 2-D herringbone
169. 2-D herringbone (in-vivo) 169 before after
170. 2-D herringbone (in-vivo)
171. Result without k-space weighting 171
172. EPI Ghost artifact Gx RO Gy PE Gz SS RF Ghost artifact image In EPI, Gradient is distorted by eddy currents and this causes phase shift Distorted gradient FT Even and odd echo mismatch causes ghost artifact! Phase shift
173. Conventional correction •  Navigator : pre-scan or reference scan Navigator-free PE RO Make phase dif ference map Navigator-based •  Pulse sequence compensation Calculate difference of phase between 1st -2nd line, 2nd -3rd line only possible to linear phase correction •  Without any modification lower performance compared to the reference-based approaches Gx RO Gy PE Gz SS RF Without PE gradient Xiang QS et al., MRM, 2007 Poser BA et al., MRM, 2013 - Using Parallel Imaging Information Kim YC et al., JMRI, 2007 Zhang et al., MRM, 2004 1) 2) 2) 1) t-1 t … … SENSE recon. SENSE recon. Phase disparity from EPI data itself Calculate - others
174. EPI model Image intensity Frequency offset Echo time Echo spacing (time between each echo) EPI data can be expressed as N : Total # of echoes n : Index of each line x : Read-out y : Phase-encoding Virtual k-space (even signals) Virtual k-space (odd signals) where Different!
175. Sparsity of difference The ghost generating phase term can be changed into a sine term Sparse How can we use this sparsity?
176. Sparsity of difference (Cont.) Low rank structured matrix completion algorithm EPI ghost correction P roblem k-space interpolation Problem using low rank structure
177. Reconstruction flow •  SE-EPI in-vivo data, 128x128 matrix size, 6/8 partial Fourier Image K-space Odd echo only Even echo only Reconstruction by odd echo Reconstruction by even echo SSOS result ALOHA
178. Result : GRE-EPI in-vivo Direct inverse FT Proposed Conventional -with reference Conventional -w/o reference Image Re-scaled Image (20%) Phase
179. Result : fMRI analysis fMRI analysis of GRE-EPI using SPM •  Pair hand squeezing stimulation ⇒ Motor cortex activation Proposed Conventional SPMmip [0, -1, 1.25] < < < SPM{T 56 } hand1 SPMresults:.\classical Height threshold T = 5.315495 {p<0.05 (FWE)} Extent threshold k = 0 voxels Design matrix 0.5 1 1.5 2 2.5 10 20 30 40 50 60 contrast 1 Statistics: p-values adjusted for search volume set-level p c cluster-level p FWE-corr q FDR-corr k E p uncorr peak-level p FWE-corr q FDR-corr T (Z ≡ ) p uncorr mm mm mm 0.000 15 0.000 0.000 110 0.000 0.000 0.000 13.22 Inf 0.000 -39 -28 50 0.000 0.000 8.75 6.92 0.000 -45 -13 61 0.000 0.019 6.82 5.80 0.000 -42 -25 65 0.000 0.000 146 0.000 0.000 0.000 13.13 Inf 0.000 45 -19 61 0.000 0.000 9.84 7.46 0.000 45 -22 46 0.000 0.001 7.86 6.43 0.000 51 -10 46 0.000 0.000 36 0.000 0.000 0.000 8.83 6.96 0.000 -51 -37 65 0.000 0.000 11 0.000 0.000 0.019 6.84 5.81 0.000 -6 53 -25 0.000 0.005 6 0.003 0.000 0.019 6.78 5.77 0.000 -39 53 -14 0.000 0.001 10 0.000 0.000 0.025 6.67 5.70 0.000 -51 47 -18 0.000 0.000 13 0.000 0.001 0.088 6.26 5.43 0.000 -3 -4 54 0.000 0.000 11 0.000 0.002 0.097 6.19 5.38 0.000 -39 -52 65 0.021 0.516 5.56 4.94 0.000 -42 -55 58 0.005 0.063 2 0.059 0.003 0.144 6.04 5.28 0.000 -33 59 5 0.005 0.063 2 0.059 0.004 0.154 6.00 5.25 0.000 -3 -94 1 0.005 0.063 2 0.059 0.005 0.188 5.92 5.20 0.000 -54 -31 50 0.005 0.063 2 0.059 0.009 0.285 5.78 5.09 0.000 6 -91 16 Proposed (multi-coil) Conventional Proposed (single-coil)
180. §  /BOPTDPQZ
     CBTFE
     PO
     MPDBMJ[BUJPO
     
     §  -PDBMJ[BUJPO
     QSFDJTJPO
     JT
     OPU
     EJGGSBDUJPO
     MJNJUFE
     §  4QBSTFMZ
     BDUJWBUFE
     QSPCFT
     
     MPDBMJ[BUJPO
     
     TVQFSS FTPMVUJPO
     JNBHF
     
     §  )PXFWFS
     TQBSTF
     BDUJWBUJPO
     TDIFNF
     IBT
     UPP
     TMPX
     U FNQPSBM
     SFTPMVUJPO
     GPS
     MJWF
     JNBHJOH
     
     §  5FOT
     PG
     TFDPOET
     PS
     TFWFSBM
     NJOVUFT
     
     §  )JHIEFOTJUZ
     JNBHJOH
     GPS
     GBTU
     MJWF
     JNBHJOH
     
     §  3FRVJSF
     B
     SPCVTU
     MPDBMJ[BUJPO
     BMHPSJUIN
     BOE
     TZTUFN
     
     6/ Low-density imaging High-density imaging Localization microscopy
181. Greedy approach Sparsity based approach Min, J.et al, Sci. Rep, 2014 Zhu, L.et al, Nat Methods, 2012 Holden, S.et al, Nat Methods, 2011 Better Localization Performance 8/ Existing high density algorithm
182. §  *NBHJOH
     NPEFM
     
     §  4VDDFTTGVM
     SFDPWFS
     JG
     TPMVUJPO
     JT
     TQBSTF
     CVU
     EJTDSFUF
     EPNBJO
     GPSNVMB
     Sub-pixel grid c g H b Camera image Convolution kernel background Point spread function Sparsity based approach g: measurement, f: sample, c:brightness, b: background, H: point spread function 9/
183. §  0VS
     QSFWJPVT
     BMHPSJUIN
     '"-$0/
     .JO
     
     FU
     
     BM
     
      
     §  GBTU
     MPDBMJ[BUJPO
     BMHPSJUIN
     CBTFE
     PO
     DPOUJOVPVT
     GPSNVMBUJPO
     §  4NBMM
     TIJGUFE
     14'
     DBO
     CF
     BQQSPYJNBUFE
     CZ
     VTJOH
     5BZMPS
     TFSJFT
     
     §  4QBSTJUZQSPNPUJOH
     MPDBMJ[BUJPO
     GPSNVMB
     PO
     DPOUJOVPVT
     TQBDF
     
     
     §  .JOJNJ[BUJPO
     CZ
     VTJOH
     GBTU
     TPMWFS
     DBMMFE
     "%..
     BMUFSOBUJOH
     EJSFD UJPO
     NFUIPE
     PG
     NVMUJQMJFST
     
     Continuous domain localization 10 /
184. Fast 2D live cell imaging results §  -JWF
     &3 FOEPQMBTNJD
     SFUJDVMVN
     1"-.
     EBUB
     §  5FNQPSBM
     
     SFTPMVUJPO
     
     TFD
     
     
     GSBNFT
     
     
     Min et al. Scientific reports (2014) 12 /
185. §  Hybrid system : Biplane + astigmatism imaging §  PSF mutual coherency test •  %.
     %JDISPJD
     NJSSPS
     •  &.
     &NJTTJPO
     GJMUFS
     
     •  #4
     #FBN
     4QMJUUFS
     •  $-
     $ZMJOESJDBM
     -FOT
     
     3D localization and hybrid system ON ON .JO
     FU
     BM
     #JPNFEJDBM
     0QUJDT
     &YQSFTT
      MFTT
     BYJBMMZ
     DPIFSFOU
     
     CFUUFS
     BYJBM
     MPDBMJ[BUJPO
     QSFDJTJPO
      
186. 3D Reconstruction of fixed cell §  'JYFE
     NJDSPUVCVMF
     4503.
     EBUB
     PG
     
     GSBNFT
     TJOHMF
     GSBNF XJEFGJFME -FBTUTRVBSF
     GJUUJOH
     '"-$0/% .JO
     FU
     BM
     #JPNFEJDBM
     0QUJDT
     &YQSFTT
       
187. 3D live cell imaging Min et al. Biomedical Optics Express (2014) §  -JWF
     &3 FOEPQMBTNJD
     SFUJDVMVN
     1"-.
     EBUB
     §  5IF
     GBTUFTU
     MJWF
     EFNPOTUSBUJPO
     TIPXJOH
     UFNQPSBM
     SFTPMVUJPO
     
     TFD
     
     
     
     15 /
188. 188 ALOHA principle
     ü  14'
     FTUJNBUJPO
     ü  %FDPOWPMVUJPO
     
     ü  (SJEGSFF
     MPDBMJ[BUJPO
     ALOHA for localization microscopy (Min et al, TIP, 2018)
189. 189/ .JOJNVN PSF estimation
190. 26 / Grid-free localization
191. 191/ True 3. Grid-free Localization 2. Deconvolution 1. PSF estimation Raw data Image Fourier ROI for PSF estimation Algorithm procedure
192. PSF variation along time
193. Reconstruction
194. Localization bias
195. Summary : Part II •  ALOHA: Measurement domain interpolation using the sparsity principle︎ •  ALOHA overcomes the many limitations of CS︎ •  Various biomedical applications of ALOHA ︎ •  Downside: computational complexity︎
196. References: Part II •  Jong Chul Ye, Jong Min Kim, Kyong Hwan Jin and Kiryung Lee, "Compressive sampling using annihilating filter-based low-rank interpolation", IEEE Trans. on Information Theory, vol. 63, no. 2, pp.777-801, Feb. 2017. •  Kyong Hwan Jin, Dongwook Lee, and Jong Chul Ye. "A general framework for compressed sensing and parallel MRI using annihilating filter based low-rank hankel matrix," IEEE Trans. on Computational Imaging, vol 2, no. 4, pp. 480 - 495, Dec. 2016. •  Kyong Hwan Jin, Ji-Yong Um, Dongwook Lee, Juyoung Lee, Sung-Hong Park and Jong Chul Ye, " MRI artifact correction using sparse + low-rank decomposition of annihilating filter-based Hankel matrix", Magnetic Resonance in Medicine (in press), 2016 •  Juyoung Lee, Kyong Hwan Jin, and Jong Chul Ye, "Reference-free single-pass EPI Nyquist ghost correction using annihilating filter-based low rank Hankel matrix (ALOHA)", Magnetic Resonance in Medicine, 10.1002/mrm.26077, Feb. 17, 2016.
197. References •  Dongwook Lee,, Kyong Hwan Jin, Eung Yeop Kim, Sung-Hong Park and Jong Chul Ye, "Acceleration of MR parameter mapping using annihilating filter-based low rank Hankel matrix (ALOHA)", Magnetic Resonance in Medicine, 10.1002/ mrm.26081, Jan. 1, 2016.​ •  Kyong Hwan Jin and Jong Chul Ye, "Annihilating filter based low rank Hankel matrix approach for image inpainting", IEEE Trans. Image Processing, 2015 Nov;24(11):3498-511. •  KH Jin, JC Ye, Sparse+ low rank decomposition of annihilating filter-based Hankel matrix for impulse noise removal, IEEE Trans. on Image Processing, 2018 •  Min, J., Carlini, L., Unser, M., Manley, S., & Ye, J. C. (2015, September). Fast live cell imaging at nanometer scale using annihilating filter-based low-rank Hankel matrix approach. In SPIE Optical Engineering+ Applications (pp. 95970V-95970V). International Society for Optics and Photonics. •  Jin, Kyong Hwan, Yo Seob Han, and Jong Chul Ye. "Compressive dynamic aperture B-mode ultrasound imaging using annihilating filter-based low-rank interpolation." Biomedical Imaging (ISBI), 2016 IEEE 13th International Symposium on. IEEE, 2016.
198. Jong Chul Ye Bio-Imaging & Signal Processing Lab. Dept. Bio & Brain Engineering Dept. Mathematical Sciences KAIST, Korea IEEE EMBS Summer School on Biomedical Imaging Sparse and Deep Learning Approaches for Biomedical Image Reconstruction: Part 3: Deep Learning Updated version of presentation material can be downloaded next week from http://bispl.weebly.com
199. Compressed Sensing Hankel Structured Matrix Comple7on Deep Learning 199 Roadmap: From CS to Deep Learning Coherent Theme of Sparse Recovery
200. Deep Learning Age •  Deep learning has been successfully used for classification, low-level computer vision, etc •  Even outperforms human observers
201. 201 Reinforcement Learning (AlphaGo)
202. Generative Adversarial Networks (GoodFellow, 2016) Figure adopted from •  BEGAN: Boundary Equilibrium Generative Adversarial Networks, 17’03 •  StackGAN : Text to Photo-realistic Image Synthesis with Stacked Generative Adversarial Networks, ‘16.12
204. •  Successful demonstraAon of deep learning for various image reconstrucAon problems –  Low-dose x-ray CT (Kang et al, Chen et al, Wolterink et al, Ye et al) –  Sparse view CT (Jin et al, Han et al, Adler et al) –  Interior tomography (Han et al) –  Stationary CT for baggage inspection (Han et al) –  CS-MRI (Hammernik et al, Schlemper et al, Yang et al, Lee et al, Zhu et al) –  US imaging (Yoon et al ) –  Diffuse optical tomography (Yoo et al) –  Elastic tomography (Yoo et al) –  Optical diffraction tomography (Kamilov et al) –  etc •  Advantages –  Very fast reconstruction time –  Significantly improved results Deep Learning for Inverse Problems
205. WavResNet results
206. WavResNet results
207. MBIR C D WavResNet results
208. 1st view 2nd view 3rd view 4th view 5th view 6th view 7th view 8th view 9th view 9-view Recon
209. FBP
210. TV
211. Deep Learning
212. WHY DEEP LEARNING WORKS FOR RECON ? DOES IT CREATE ANY ARTIFICIAL FEATURES ?
213. First CNN: LeNet (LeCun, 1998)
214. Convolutional Layer
215. Pooling Layer Figure from Leonardo Araujo Aantos
216. Too Simple to Analyze..? Convolution & pooling à stone age tools of signal processing What do they do ?
217. Dark Age of Applied Mathematics ?
218. CNN – BIOLOGICAL ORIGIN A LAYMAN’S EXPLANATION
219. 219 http://klab.smpp.northwestern.edu/wiki/images/4/43/NTM2.pdf Emergence of Hiearchical Features
220. 220 http://www.vicos.si/File:Lhop-hierarchy-second.jpg •  LeCun et al, Nature, 2015 Hierarchical representation
221. Visual Information Processing in Brain 221 Kravitz et al, Trends in Cognitive Sciences January 2013, Vol. 17, No. 1
222. Retina, V1 Layer 222 Receptive fields of two ganglion cells in retina à convolution Orientation column in V1 http://darioprandi.com/docs/talks/image-reconstruction-recognition/graphics/pinwheels.jpg Figure courtesy by distillery.com
223. Visual Pathway 223 Poggio et al NATURE|VOL 431 | 14 OCTOBER 2004
224. “The Jennifer Anniston Cell” 224 Quiroga et al, Nature, Vol 435, 24, June 2005
226. CNN – MATHEMATICAL UNDERSTANDING
227. Why Deep Learning works for recon ? Existing views 1: unfolded iteration •  Most prevailing views •  Direct connecAon to sparse recovery –  Cannot explain the filter channels Jin, arXiv:1611.03679
228. Why Deep Learning works for recon ? Existing views 2: generative model •  Image reconstruc;on as a distribu;on matching –  However, difficult to explain the role of black-box network Bora et al, Compressed Sensing using Generative Models, arXiv:1703.03208
229. •  What is the role of the nonlinearity such as rectified linear unit (ReLU) ? •  Why do we need a pooling and unpooling in some architectures ? •  Why do some networks need fully connected layers whereas the others do not ? •  What is the role of by-pass connection or residual network ? •  What is the role of the filter channels in convolutional layer ? Many Mysteries…
230. Skipped or Not ? Residual Network Clean image Standard Network Zhang, K., et al, IEEE TIP, 2017.
231. Even Confused about Where to Start Optimization Generalization Architecture
232. Our Proposal: Deep Learning == Deep Convolutional Framelets •  Ye et al, “Deep convolutional framelets: A general deep learning framework for inverse problems”, SIAM Journal Imaging Sciences, 11(2), 991-1048, 2018.
233. Matrix Representation of CNN Figure courtesy of Shoieb et al, 2016
234. Hankel Matrix: Lifting to Higher Dimensional Space
235. Why we are excited about Hankel matrix ? T -T 0 n1 -n1 0 * FRI Sampling theory (Ve#erlie et al) and compressed sensing
236. ︙ ︙ 1 2 3 4 5 6 7 8 9 -1 0 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 12 2 3 4 5 6 7 8 9 12 13 3 4 5 6 7 8 9 10 10 10 10 0 11 11 11 1 2 3 4 5 Finite length convolution Matrix Representation * ALOHA : Annihilating filter based LOw rank Hankel matrix Approach * Jin KH et al. IEEE TCI, 2016 * Jin KH et al.,IEEE TIP, 2015 * Ye JC et al. IEEE TIT, 2016 Annihilating filter-based low-rank Hankel matrix
237. Missing elements can be found by low rank Hankel structured matrix compleAon Nuclear norm Projec;on on sampling posi;ons min m kH ( m ) k⇤ subject to P⌦(b) = P⌦( f ) RankH(f) = k * Jin KH et al IEEE TCI, 2016 * Jin KH et al.,IEEE TIP, 2015 * Ye JC et al., IEEE TIT, 2016 m Annihilating filter-based low-rank Hankel matrix
238. Key Observation Data-Driven Hankel matrix decomposition => Deep Learning •  Ye et al, “Deep convolutional framelets: A general deep learning framework for inverse problems”, SIAM Journal Imaging Sciences, 11(2), 991-1048, 2018.
239. H d(f) = U⌃V T : Non-local basis : Local basis Convolution Framelets (Yin et al; 2017) > = I > = I H d(f)
240. H d(f) H d(f) = ˜ T ˜ T C C = T H d(f) C = T (f ~ ) Encoder: ˜ T = I ˜ = PR(V ) H d(f) = U⌃V T Unlifting: f = (˜C) ~ ⌧(˜ ) : Non-local basis : Local basis : Frame condition : rank condition convolution pooling un-pooling convolution : User-defined pooling : Learnable filters H pi (gi) = X k,l [Ci]kl e Bkl i Decoder: Deep Convolutional Framelets (Y, Han, Cha; 2018)
241. Single Resolution Network Architecture
242. Multi-Resolution Network Architecture
243. Network Width versus Depth H d(f) Filter length > Intrinsic complexity of signal = Hankel matrix rank Filter length < Intrinsic complexity of signal ? à recursive application of the extended Hankel matrix
244. Network Width versus Depth
245. Deep Pyramidal Residual Network Han et al, CVPR, 2017
246. SO FAR, THE THEORY IS LINEAR.. ROLE OF NONLINEARITY ?
247. Role of Nonliearity: ReLU à Conic coding D. D. Lee and H. S. Seung, Nature, 1999 ReLU: positive framelet coefficients Conic encoding à part by representation similar to visual cortex Han et al, 2018 arXiv:1805.03779
248. Conic fi [Ci]kl 0 H pi (gi) = X k,l [Ci]kl e Bkl i H pi (fi) ' Lifting Geometry of Neural Network gi Un-lifting Ci(fi) Ci(fi) 0 i ⇣i ⇣ e i ⌘
249. fi [Ci]kl 0 H pi (gi) = X k,l [Ci]kl e Bkl i H pi (fi) ' Lifting Geometry of Residual Neural Network Ci(fi) Ci(fi) 0 i ⇣i ⇣ e i ⌘ gi Un-lifting
250. https://www.youtube.com/watch?v=DdC0QN6f3G4 Relativity: Lifting to the Space-time à linear trajectory High Dimensional Lifting in Relativity Falling apple Apple at rest Universal law of gravity: 3-D Space à curved trajectory Falling apple Apple at rest
251. Deep CNN Lifting Un-lifting Conic Lifting Un-lifting Lifting Un-lifting
252. Successive Low-Rank Hankel Matrix Approximation Truncated channel à Low rank Hankel matrix approximation
253. Compressed Sensing Hankel Structured Matrix Comple7on Deep Learning 253 From CS to Deep Learning: Coherent Theme
254. APPLICATION-DRIVEN EVIDENCES
255. U-Net Architecture (Ronnenberger, 2015)
256. Problem of U-net Pooling does NOT satisfy the frame condition JC Ye et al, SIAM Journal Imaging Sciences, 2018 Y. Han et al, TMI, 2018. ext > ext = I + > 6= I
257. Improving U-net using Deep Conv Framelets •  Dual Frame U-net •  Tight Frame U-net JC Ye et al, SIAM Journal Imaging Sciences, 2018 Y. Han and J. C. Ye, TMI, 2018
258. U-Net versus Dual Frame U-Net
259. Tight-Frame U-Net JC Ye et al, SIAM Journal Imaging Sciences, 2018
260. Denoising: U-Net vs. Tight-Frame U-Net
261. Inpainting: U-Net vs. Tight-Frame U-Net
262. k-Space Deep Learning for Accelerated MRI Han et al, arXiv:1805.03779 ALOHA : k-space interpolation à k-space interpolation using deep learning ? Yes
263. k-Space Deep Learning for Accelerated MRI Han et al, arXiv:1805.03779 ~3dB gain
264. k-Space Deep Learning for Parallel MRI Han et al, arXiv:1805.03779
265. k-Space Deep Learning for Parallel MRI Cha et al, arXiv:1806.00806 ~13dB gain
266. k-Space Deep Learning for Parallel MRI Cha et al, arXiv:1806.00806
267. DEEP NETWORKS FOR CT
268. Low-Dose CT •  To reduce the radiation exposure, sparse-view CT, low-dose CT and interior tomography. Sparse-view CT (Down-sampled View) Low-dose CT (Reduced X-ray dose) Interior Tomography (Truncated FOV)
269. 269 Wavelet transform level 2 level 1 level 3 level 4 Wavelet recomposition + Residual learning : Low-resolution image bypass High SNR band CNN (Kang, et al, Medical Physics 44(10))
270. Routine dose Quarter dose (Kang, et al, Medical Physics 44(10) 2017)
271. Routine dose AAPM-Net results (Kang, et al, Medical Physics 44(10) 2017)
272. WavResNet results (Kang et al, TMI, 2018)
273. WavResNet results (Kang et al, TMI, 2018)
274. MBIR Our latest Result C D WavResNet results
275. Full dose Quarter dose
276. Full dose Quarter dose
277. Sparse-View CT •  To reduce the radiation exposure, sparse-view CT, low-dose CT and interior tomography. Sparse-view CT (Down-sampled View) Low-dose CT (Reduced X-ray dose) Interior Tomography (Truncated FOV)
278. Multi-resolution Analysis & Receptive Fields
279. Tight-Frame U-Net JC Ye et al, SIAM Journal Imaging Sciences, 2018 Han et al, TMI, 2018
280. 90 view recon U-Net vs. Tight-Frame U-Net •  JC Ye et al, SIAM Journal Imaging Sciences, 2018 •  Y. Han and J. C. Ye, TMI, 2018
283. •  Figures from internet 9 View CT for Baggage Screening
284. 9 View CT for Baggage Screening
285. 9 View Dual Energy CT for Baggage Screening Han et al, arXiv preprint arXiv:1712.10248, (2017). CT meeting 2018
286. 1st view 2nd view 3rd view 4th view 5th view 6th view 7th view 8th view 9th view
287. FBP
288. TV
289. Ours
290. ROI Reconstruction •  To reduce the radiation exposure, sparse-view CT, low-dose CT and interior tomography. Sparse-view CT (Down-sampled View) Low-dose CT (Reduced X-ray dose) Interior Tomography (Truncated FOV)
291. Deep Learning Interior Tomography Han et al, arXiv preprint arXiv:1712.10248, (2017): CT meeting 2018.
292. Ground Truth FBP
293. TV Chord Line Ours 8~10 dB gain
294. Ground Truth Chord Line TV Ours
295. Summary: Part III •  Deep learning for inverse problems︎ •  Significant performance gain︎ •  Deep convolutional framelets: new mathematical tools for understanding deep neural network for inverse problems︎ •  Depth vs. width︎ •  Role of ReLU︎ •  Residual net︎ •  Numerical results support the theory︎
296. Still Unresolved Problems.. •  Cascaded geometry of deep neural network •  Generalization capability •  Optimization landscape •  Training procedure •  Extension to classification problems
297. LAB CONTENTS
298. Dataset of Natural Images •  MNIST (http://yann.lecun.com/exdb/mnist/) –  Handwritten digits. •  SVHN (http://ufldl.stanford.edu/housenumbers/) –  House numbers from Google Street View. •  ImageNet (http://www.image-net.org/) –  The de-facto image dataset. •  LSUN (http://lsun.cs.princeton.edu/2016/) –  Large-scale scene understanding challenge. •  Pascal VOC ( http://host.robots.ox.ac.uk/pascal/VOC/) –  Standardized image dataset. •  MS COCO (http://cocodataset.org/#home) –  Common Objects in Context. •  CIFAR-10 / -100 (https://www.cs.utoronto.ca/~kriz/cifar.html) –  Tiny images data set. •  BSDS300 / 500 ( https://www2.eecs.berkeley.edu/Research/ Projects/CS/vision/grouping/resources.html) –  Contour detection and image Segmentation resources. https://en.wikipedia.org/wiki/List_of_datasets_for_machine_learning_research https://www.kaggle.com/datasets
299. Dataset for Medical Images •  HCP (https://www.humanconnectome.org/) –  Behavioral and 3T / 7T MR imaging dataset. •  MRI Data (http://mridata.org/) –  Raw k-space dataset acquired on a GE clinical 3T scanner. •  LUNA (https://luna16.grand-challenge.org/data/) –  Lung Nodule analysis dataset acquired on a CT scanner. •  Data Science Bowl (https://www.kaggle.com/c/data-science-bowl-2017) –  A thousand low-dose CT images. •  NIH Chest X-rays (https://nihcc.app.box.com/v/ChestXray-NIHCC) –  X-ray images with disease labels. http://www.cancerimagingarchive.net/ https://www.kaggle.com/datasets •  TCIA Collections (http://www.cancerimagingarchive.net/) •  De-identifies and hosts a large archive of medical images of cancer accessible for public download. •  The data are organized as “Collections”, typically patients related by a common disease, image modality (MRI, CT, etc).
300. Libraries for Deep learning •  TensorFlow (https://www.tensorflow.org/) –  Python •  Theano ( http://deeplearning.net/software/theano/) –  Python •  Keras (https://keras.io/) –  Python •  Caffe (http://caffe.berkeleyvision.org/) –  Python •  Torch (or PyTorch) (http://torch.ch/) –  C / C++ (or Python) •  Deeplearning4J ( https://deeplearning4j.org/) –  Java •  Microsoft Cognitive Toolkit (CNTK) ( https://www.microsoft.com/en-us/cognitive- toolkit/) –  Python / C / C++ •  MatConvNet ( http://www.vlfeat.org/matconvnet/) –  Matlab
301. Implementation on MatConvNet Download the MatConvNet Provide the pre-trained models
302. Step 1: Compile the toolbox 1.  Unzip the MatConvNet toolbox. 2.  Open ‘vl_compilenn.m’ in Matlab.
303. Step 1: Compile the toolbox (cont.) 3.  Check the options such as enableGpu and enableCudnn. 4.  Run the ‘vl_compilenn.m’. * To use GPU processing (false true), you must have CUDA installed. ( https://developer.nvidia.com/cuda-90-download-archive ) ** To use cuDNN library (false true), you must have cuDNN installed. ( https://developer.nvidia.com/cudnn )
304. Step 2: Prepare dataset As a classification example, MNIST consists of as follows, Images % struct-type Data % handwritten digit image labels % [1, …, 10] set % [1, 2, 3], 1, 2, and 3 indicate train, valid, and test set, respectively. Data Labels 6
305. Step 2: Prepare dataset (cont.) •  As a segmentation example, U-Net dataset consists of as follows, –  Images % struct-type Ø  data % Cell image Ø  labels % [1, 2], Mask image. 1 and 2 indicate back- and for-ground, respectively. Ø  set % [1, 2, 3] Data Labels
306. Step 3: Implementation of the network architecture •  Developers only need to program the network architecture code because MatConvNet supports the network training framework. Support famous network architectures, such as alexnet, vggnet, resnet, inceptionent, and so on.
307. Step 3: Implementation of the architecture (cont.) –  The implementation details of U-Net U-Net can be implemented, recursively. Stage 0 Stage 1 Stage 2 Stage 3 Stage 4
308. Step 3: Implementation of the architecture (cont.) 1.  Create objects of network and layers. Encoder Part Skip + Concat Part Decoder Part •  The structure of Stage 0 Network Part
309. Step 3: Implementation of the architecture (cont.) 2.  Connect each layers. •  The structure of Stage 0 Layer Name ( string-type ) Layer object ( object ) Input Name ( string-type ) Output Name ( string-type ) Parameters Name ( string-type ) All objects and names must be unique.
310. Step 3: Implementation of the architecture (cont.) 3.  Implement recurrently the each stages and add a loss function. Previous parts (3.1 and 3.2) become functional as ‘add_block_unet’.
311. Step 4: Network hyper-parameter set up •  MatConvNet supports the default hyper-parameters as follows, Refer the cnn_train.m ( or cnn_train_dag.m ) The supported hyper-parameters 1.  The size of mini-batch 2.  The number of epochs 3.  Learning rate 4.  Weight decay factor 5.  Solvers such as SGD, AdaDelta, AdaGrad, Adam, and RMS The kind of Optimization Solvers
312. Step 5: Run the training script 1. Training script 2. Training loss 3. Training loss graph •  Blue : train •  Orange : valid
313. Cell segmentation
314. train_cnn_cell_seg.m val_cnn_cell_seg.m
315. F5 ژח Run // Training code Training code
316. Ѿҗ : figure 1
317. Ѿҗ : figure 10
318. F5 ژח Run ߡౡ ௿ܼ // Validation code Validation code // ޷ܻ ੷੢೧֬਷ ֎௼ਕ௼ܳ ࠛ۞৬ प೯
319. F5 ژח Run ߡౡ ௿ܼ Validation code // ೟णػ ݃૑݄ ֎௼ਕ௼ܳ ࠛ۞৬ प೯
320. Ѿҗ : figure 20
321. References: Part III •  Ge Wang, Jong Chu Ye, Klaus Mueller, Jeffrey A Fessler, "Image Reconstruction Is a New Frontier of Machine Learning", IEEE Trans. on Medical Imaging, Vol. 37 no. 6, pp. 1289 - 1296, June 2018. •  Yoseob Han and Jong Chul Ye,"Framing U-Net via Deep Convolutional Framelets: Application to Sparse-view CT", Special Issue on Machine Learning for Image Reconstruction, IEEE Trans. on Medical Imaging, vol. 37, no. 6, pp. 1418-1429, June, 2018. •  Eunhee Kang, Won Chang, Jaejun Yoo, and Jong Chul Ye,"Deep Convolutional Framelet Denosing for Low-Dose CT via Wavelet Residual Network", Special Issue on Machine Learning for Image Reconstruction, IEEE Trans. on Medical Imaging, vol. 37, no.6, pp. 1358-1369, 2018. •  Dongwook Lee, Jaejun Yoo, Sungho Tak and Jong Chul Ye, "Deep Residual Learning for Accelerated MRI using Magnitude and Phase Networks", IEEE Trans on Biomedical Engineering (in press), Invited paper for Special Section on Deep Learning, 2018. •  Jong Chul Ye, Yoseob Han and Eunju Cha, "Deep convolutional framelets: a general deep learning framework for inverse problems", SIAM Journal on Imaging Sciences (in press), 2018. •  Yoseob Han, Jaejun Yoo, Hak Hee Kim, Hee Jung Shin, Kyunghyun Sung, and Jong Chul Ye, "Deep Learning with Domain Adaptation for Accelerated Projection-Reconstruction MR", Magnetic Resonance in Medicine (accepted), 2018. •  Eunhee Kang, Junhong Min and Jong Chul Ye, " A Deep Convolutional Neural Network using Directional Wavelets for Low-dose X-ray CT Reconstruction", Medical Physics 44.10 (2017).
322. Acknowledgements CT Team •  Yoseob Han •  Eunhee Kang •  Jawook Goo US Team •  Shujaat Khan •  Jaeyong Hur MR Team •  Dongwook Lee •  Juyoung Lee •  Eunju Cha •  Byung-hoon Kim Image Analysis Team •  Boa Kim •  Junyoung Kim Optics Team •  Sungjun Lim •  Junyoung Kim •  Jungsol Kim •  Taesung Kwon