Sıkıştırılmış Çoklu-Kontrast Manyetik Rezonans Görüntülerinin Genişletilmiş Lagrange Metodu ile Gerikazanımı Compressed Multi-Contrast Magnetic Resonance Image Reconstruction using Augmented Lagrangian Method Alper Güngör1,∗ , Emre Kopanoğlu1,∗ , Tolga Çukur2 , H. Emre Güven1 Özetçe —Bu bildiride, çok-kanallı/çoklu-kontrast görüntüleme problemleri için birden fazla düzenlileştirme terimine sahip kısıtlı enküçültme problemine Genişletilmiş Lagrange Yöntemleri temelli bir çözüm sunulmaktadır. Önerilen algoritma verilen problemi çok kanal için birleşik hedef fonksiyonları ile ortak bilgiden faydalanarak benzerlerinden daha yüksek kalitede gerikazanımlar yapmaktadır. Grup seyreklik fonksiyonu ile renkli toplam değişinti fonksiyonlarının birlikte kullanımı önerilmiş ve benzerleriyle yakınsama hızı ve görüntü kalitesi bakımından kıyaslanmıştır. Son olarak, algoritma Manyetik Rezonans Görüntüleme (MRG) verileri üzerinde gösterilmiş ve benzer bir yöntem olan “Fast Composite Splitting Algorithm (FCSA-MT)” [1] ile pSNR (tepe Sinyal Gürültü Oranı) bakımından kıyaslanmıştır. Anahtar Kelimeler—Görüntü gerikazanımı, Sıkıştırılmış Algılama, Yön Değiştiren Çarpanlar Yöntemi, Seyreklik Abstract—In this paper, a Multi-Channel/Multi-Contrast image reconstruction algorithm is proposed. The method, which is based on the Augmented Lagrangian Method uses joint convex objective functions to utilize the mutual information in the data from multiple channels to improve reconstruction quality. For this purpose, color total variation and group sparsity are used. To evaluate the performance of the method, the algorithm is compared in terms of convergence speed and image quality using Magnetic Resonance Imaging data to FCSA-MT [1], an alternative approach on reconstructing multi-contrast MRI data. Keywords—Image Reconstruction, Compressed Sensing, Alternating Direction Method of Multipliers, Sparsity I. G İR İ Ş Bu bildiride, çok-kanallı görüntülerin birlikte gerikazanımında kısıtlı eniyileme teknikleri kullanılması ele alınmaktadır. Önerilen teknik, birden çok özellikten faydalanabilmesinin yanısıra, çok kanallı görüntülerde bulunan ortak özelliklerden de faydalanabilmektedir. Çeşitli MRG sıkıştırılmış algılama (SA) uygulamaları toplam değişinti ve dönüşüm uzayı seyrekliği gibi birden çok özelliğin bir arada artırılmasını sağlayan eniyileme problemleri çözümünü içermektedir [2], [3]. Yakın zamanlı çalışmalar, parçalı-sabit ve seyrek özelliklerinden grup seyreklik fonksiyonları kullanılarak daha fazla fayda sağlanabileceğini göstermiştir [1]. Görüntü, görüntü uzayında grup-seyrek ve gradyan uzayında grup-seyrek özelliklerden oluşacak şekilde modellenmiştir. Burada, Yön Değiştiren Çarpanlar Yöntemi (YDÇY) ile önerilen problemin hesapsal olarak verimli bir şekilde çözümü önerilmektedir. YDÇY teknikleri sinyal ve görüntü gerikazanımı problemlerine başarıyla uygulanmaktadır [4], [5]. Bu bildiride çok kanallı görüntü gerikazanımı probleminin YDÇY kullanılarak çözülmesi için bir algoritma önerilmektedir. Algoritmanın hızlı gerçeklenmesi sunulmakta, etkili kullanımı MRG verisi üzerinde gösterilmektedir. Aynı zamanda algoritma benzer bir yöntem olan Fast Composite Splitting Algorithm (FCSAMT) [1], kıyaslanmakta, süre ve görüntü kalitesi sonuçları verilmektedir. Görüntü kalitesi ölçütleri tepe Sinyal Gürültü Oranı (pSNR) ve karekök ortalama hatası (RMSE)’dir. II. G ENEL B ILGILER A. Veri modeli Birçok görüntüleme problemi görüntü vektörünün, x ∈ CN , gözlem vektörü y ∈ CM ile ilişkilendirilmesinde doğrusal operatörler kullanılarak modellenebilir. Gözlem matrisi B ise CM ×N kümesinin bir elemanıdır. Gürültü vektörünün, n ∈ CM , eklenmesiyle problem şu şekilde ifade edilebilmektedir: 1 ASELSAN Araştırma Merkezi, Ankara, Türkiye Elektrik-Elektronik Mühendisliği, İhsan Doğramacı Bilkent Üniversitesi, Ankara, Türkiye ∗ Bu yazarlar çalışmaya eşit katkıda bulunmuştur. alpergungor@aselsan.com.tr, ekopanoglu@aselsan.com.tr, cukur@ee.bilkent.edu.tr, heguven@aselsan.com.tr 2 978-1-5090-1679-2/16/$31.00 c 2016 IEEE y = Bx + n, (1) n genellikle normal dağılıma sahiptir. Bunun yanında, bazı görüntüleme problemleri çok-kanallı gerikazanım probleminin çözümüne ihtiyaç duymaktadır. Gözlem vektörlerinin kaynağı Yeşil, Kırmızı, Mavi kanallı bir kamera gibi farklı olsa da, bu kanalların birlikte gerikazanımı genellikle performansı arttırabilmektedir [8], [9]. i numaralı kanal için, düzenleyici denklem: y(i) = B(i) x(i) + n(i) , (2) olacaktır. Burada y(i) , B(i) , x(i) , ve n(i) değişkenleri, i numaralı kanal için denklem (1) ile aynı değişkeni ifade etmektedir. Bu bildiride MRG verisinin çok kanallı veri üzerinden gerikazanımı yapılmaktadır. MRG ile elde edilen Proton Yoğunluğu (PY), T1-ağırlıklı (T1a), ve T2-ağırlıklı (T2a) görüntülerine farklı kanaldan gelen veriler gibi davranılmıştır [1]. MRG gerikazanımı için, gözlem matrisi B(i) , kısmi Fourier gözlemleri ile ilişkili doğrusal operatördür. [7]’de, gözlem matrisi B’nin, aynı Fourier dönüşüm uzayı gibi, birimsel dönüşüm uzayında maskeleme olarak ifade edildiği durumlar için algoritmanın karmaşıklığının azaltılması için bir yöntem tarif edilmektedir. III. YÖNTEM Bu bölümde daha önce [7]’de önerilen hibrit YDÇY’nin çok-kanallı görüntülerin gerikazanımında kullanılması için genişletilmesi tarif edilmektedir. Beynin çok-kanallı –PY, T1a, T2a– MRG taramalarının alındığı durum incelenmiştir. A. Çok kanallı YDÇY B. SA gerikazanım yaklaşımları Gerikazanım işlemi kısıtlı bir eniyileme probleminin çözümünü içerir. Bu çalışmada, parametre seçiminin daha kolay olması sebebiyle kısıtlı biçim kullanılmıştır [7]. Eniyileme problemi şu biçimde ifade edilebilir: minimize α1 φ1 (x) + α2 φ2 (x) x , subject to kBx − yk2 ≤ ǫ (3) YDÇY’nin zorlayıcı olmayan koşullar altında yakınsaması garanti edilmektedir [4]. Hibrit IRWALM [7] az sayıda örneklenmiş veriden çok özellikli verinin gerikazanımı için kullanılan YDÇY tabanlı bir eniyileme yöntemidir (Algoritma 1). Algoritma öncelikle ayrılabilir her fonksiyon φi (·) için fonksiyonları ayırmaktadır ve her yinelemede Moreau yakınlık operatörlerini uygulamaktadır: µ Ψφ (v) = proxφ(·) (v) = arg min φ(x) + kx − vk22 . (4) x 2 Daha sonra sonuçlar minimum kareler adımı ile birleştirilmektedir. Algoritma 1: Hibrit IRWALM 1. Eşitle: k = 0, β0 = 1, (1) (2) (3) Seç: µ > 0, η ∈ (0, 1), z0 , z0 , z0 , (1) (2) (3) d0 , d0 , d0 , p, α1 , α2 2. tekrarla
(1) (1) (2) (2) (3) (3) 3. rk = zk + dk + zk + dk + BH zk + dk 4. xk+1 = I + B B (2) 5. zk+1 = Ψ α1 φ µ µ (3) rk

(2) (1) 1 xk+1 − dk 2 xk+1 − dk (2) 6. zk+1 = Ψ α2 φ
−1 (3) 7. zk+1 = ΨιE(ǫ,I,y) Bxk+1 − dk (1) (1) (1) 8. dk+1 = dk − xk+1 + zk+1 (2) (2) (2) 9. dk+1 = dk − xk+1 + zk+1 (3) (3) (3) 10. dk+1 = dk − Bxk+1 + zk+1 minimize x α1 φ1 (x) + · · · + αm φm (x) (i) (i) subject to kB x burada φi (x) ayrılabilir görüntü ön bilgisini, ǫ veri bağlılığı hata limitini ifade etmektedir. [2], [6], [7], [10] benzer gerikazanım problemini hem ℓp -norm (kxkpp ) hem de büyüklüğün Toplam Değişintisini (T V (|x|)) birlikte enküçülterek yapmaktadır. H Her kanalın bir kontrast ağırlığı ile ilintilendirildiği k kanal, ve m ayrılabilir hedef fonksiyonu için, problem şu şekilde ifade edilebilir:
11. k ← k + 1 12. ta ki: durma kriteri sağlanana kadar. Burada ΨιE(ǫ,I,y) (s), kBx − yk2 ≤ ǫ kısıtı ile ilintili Moreau yakınlık fonksiyonudur, ve ℓ2 -norm hiperküresine basit bir izdüşüm işlemi olarak uygulanabilir. Ayrıca, yine − y(i) k2 ≤ ǫi , i ∈ 1, · · · , k . (5) Her hedef fonksiyonu φi ya bir kanala özel ve herhangi bir xhi ile ilintili bir ifonksiyon, ya da tüm görüntü vektörü T x = x(1) · · · x(k) fonksiyonu olabilir. T T vektörünü içeren bir grup seyrekliği Bu problemin YDÇY kullanılarak çözülmesi için aşağıdaki YDÇY değişkenleri tanımlanır. i ∈ 1, ·i· · , k kanalları için, z(i) = h Bütün T T T (i) T T (i,m) T (i,1) (i,0) , ve Gi = [I · · · I (B ) ] . ···z z z G matrisini diyagonallerinde Gi bulunan bir blok-diyagonal matris olarak tanımlandığında, aşağıdaki problem çözülebilir: minimize x,z f1 (x) + f2 (z) subject to Gx + Qz − r = 0, (6) Q = −I, r = 0, f1 (x) = 0,   Pm f2 (z) = z(i,t) i=1,··· ,k + t=1 αt φt
Pk (i,0) . Bu tanımlar Gx = z, ve sonuç i=1 ιE(ǫi ,I,y(i) ) z (i,t) olarak: x(i) = zt=1,··· ,m , and B(i) x(i) = z(i,0) kısıtlarını sağlar. Bu tanımlar, grup seyrekliği, renk toplam değişintisi gibi grup fonksiyonlarının yanısıra, ℓ1 -norm ve toplam değişinti gibi ayrı hedef fonksiyonlarının da kullanımına imkan sağlamaktadır. YDÇY’nin adımları izlendiğinde, Algoritma 2’ye ulaşılabilir. Burada par-tekrar algoritmanın birbirinden bağımsız olarak tekrarlanabilen adımlarını ifade etmektedir. Burada bütün i değerleri için, 4,5, ve 6 numaralı adımlar paralel olarak hesaplanabilir. Ayrıca, [7]’de gösterildiği gibi, MRG’deki gibi kısmi-Fourier gözlemleri için 4. adım (7) kullanılarak hızlı hesaplanabilir.   −1  1 1 (i) H (i) (i) H (i) mI + (B ) B = I− . (B ) B m m+1 (7) Algoritma 2: Çok Kanallı Hibrit YDÇY Renk toplam değişintisi [8], şu şekilde tanımlanmıştır: v u k X uX t (∆ x(i) )[a, b]2 + (∆ x(i) )[a, b]2 . CT V (x) = 1 t 2 t 1. Bütün i ∈ 1, · · · , k ve t ∈ 0, · · · , m için, (i,t) Eşitle: n = 0, Seç: µ > 0, z0 (i,t) , d0 , αt 2. tekrarla a,b 3. par-tekrar i = 1, · · · , k (i) xn+1 4. = mI + (B(i) )H B Pm t=1 (i,t) (i,t) dn + zn + BH d(i,0) + z(i,0) n n 5. (i,0) zn+1 6. dn+1 = d(i,0) − B(i) xn+1 + zn+1 n = Ψι E ǫ,I,y(i) (i,0)
B(i) x(i) − d (i) (0)
(i,0) Renk toplam değişintisi için, Chambolle [11] tarafından önerilen enküçültme algoritması kullanılmaktadır ve bu algoritma yakınsayana kadar değil sabit sayıda bir iterasyon sayısında çalıştırılmaktadır. Böylece, çözülen problem aşağıdaki gibi ifade edilebilir:
 7. par-tekrar sonu minimize 8. par-tekrar t = 1 , · · · , m 9. 10.  (i,t) zn+1 i=1,··· ,k (i,t) (i,t) dn+1 = dn − = Ψφ (i) αt t µ  (i,t) xn+1 (t,i) − dn i=1,··· ,k (i,t) xn+1 + zn+1 , bütün t,i için x  IV. 12. n ← n + 1 13. ta ki durma kriteri sağlanana kadar... x ve Bx hesaplanması (9) ve (10) kullanılarak daha da hızlandırılabilir. q(i) n = m X + z(i,t) d(i,t) n n . (8) t=1 (i) B(i) xn+1  1 1 = q(i) + (B(i) )H · m n m+1    i − Bqn , m d(i,0) + z(i,0) n n i 1 h (i) (i)  (i,0) = . B qn + dn + z(i,0) n m+1 α1 kxk2,p + αCT V CTV (x) (i) (i) subject to kB x 11. par-tekrar sonu (i) xn+1 i=1 (13)
(i) −1 (9) (10) Böylece, algoritma yineleme başına 2 hızlı Fourier dönüşümü işlemi ile gerçeklenebilir. Dışbükey bir φi ile kullanıldığında algoritmanın yakınsama garantisi bulunacaktır. B. MRG’ye Uygulanması Çoklu-kontrast MRG için, daha önceki bir uygulamada Renkli Toplam Değişinti (CTV) ve bir dönüşüm uzayında grup seyreklik fonksiyonu kullanılmıştır [1]. Burada önerilen yöntemle algoritmanın daha hızlı gerçeklenmesi sağlanmaktadır. Altta yatan anatominin aynı oluşundan kaynaklı olarak, farklı kontrasta sahip görüntülerin bir dönüşüm uzayında benzer seyreklik terimlerine sahip olması beklenmektedir. Grup seyrekliği, farklı görüntülerde aynı konumlardaki seyrekliği birlikte öne çıkarır. Yukarıdaki notasyonu kullanarak, burada “grup” tanımını farklı görüntülerde aynı konumlarda bulunan pikseller, yani şu şekilde yapabiliriz: h iT xt = x(1) [t] · · · x(k) [t] . (11) − y(i) k2 ≤ ǫi , i ∈ 1, · · · , k . (14) S ONUÇLAR Algoritma 2, Chambolle izdüşümü işlemi mex fonksiyonundan çağırılarak MATLAB’da gerçeklenmiştir. Deneyler Intel Xeon E5-2650 v2 işlemcili ve 64 GB RAM’e sahip bir iş istasyonunda yapılmıştır. Başlangıçta zi,t 0 vektörleri sıfırladoldurulmuş görüntü gerikazanımlarına eşitlenmiştir, bütün t = 1, · · · , m, zi,0 = y(i) için zi,t = (B(i) )H y(i) ve 0 0 i,t bütün Lagrange değişkenleri d0 sıfır vektörlerine eşitlenmiştir. Simülasyonlar çalıştırılmadan önce bütün görüntüler [0, 255] arasına normalize edilmiştir ve µ hızlı ve doğru yakınsama için 0.1/256 olarak seçilmiştir. Her yinelemede µ, [4]’ün önerdiği üzere daha iyi bir yakınsama için 1.02 ile çarpılmıştır. Durma kriteri “maksimum yineleme sayısına ulaşılana kadar (i) (i) veya kxk+1 − xk k22 < 10−3 ” olarak seçilmiştir. İlk olarak, Monte Carlo benzetimleri FCSA-MT [1], ayrık Hibrit YDÇY [7], ve önerilen yöntem üzerinde farklı örnekleme şablonları kullanılarak koşturulmuştur. İki farklı gerçekdeğerli görüntü veri kümesi kullanılmıştır: %12.5 ve %6.25’i örneklenen bir beyin silüeti (BS) benzetimi ve %25 ve %6.25’e örneklenen SRI atlas (SR) görüntüsü [1]. Bütün deneyler için üç kanal da eşit ölçüde alt örneklenmiştir. α1 ve αCT V için daha önce FCSA-MT’de seçilen sırasıyla 0.0173 ve 0.0606 değerleri kullanılmıştır [1]. ǫ için ise, hiç bir RF uyarılması olmadan yapılan bir MRG veri edinmesi simüle edilerek gürültü verisi toplanmış, ve ǫ toplanan verinin ℓ2 -normuna eşitlenmiştir. Uygun bir karşılaştırma için yayınlanan FCSAMT kodu kullanılmış ve maksimum iterasyon sayısı FCSAMT’nin çalışma süresine eşit olacak şekilde seçilmiştir. Bu sayı önerilen yöntem için 200 olmuştur. (12) Tablo 1’de farklı algoritma ve parametreler için pSNR değerleri gösterilmektedir. Tabloda görüldüğü üzere, en yüksek başarım önerilen algoritma tarafından elde edilmiştir. Ayrık kanal gerikazanımı yüksek alt-örneklemeden dolayı kötü performansa sahip olmaktadır. Şekil 1 (a) ve (b)’de farklı alt örnekleme şablonları için pSNR değerinin ortalaması ve varyansı %12.5 ve %6.25’lik BS verisi için gösterilmektedir. pSNR’daki yüksek varyans, başarımdaki belirsizliği arttırdığı için istenmemektedir. Bu durumda da, önerilen algoritma en düşük varyans ve en yüksek ortalamaya sahip olduğu için en iyi başarımı elde etmiştir. Burada tüm işlemler eleman üzerindendir ve hızlı bir şekilde gerçeklenebilir. Şekil 2’de %12.5’lik alt örnekleme için BS gerikazanım sonuçları gösterilmektedir. FCSA-MT ile ayrık YDÇY’nin kıyası, birlikte gerikazanımın çoklu-kontrast görüntü için O zaman ilgili yakınlık fonksiyonu şu şekilde tanımlanabilir: S1 (xt , α) = xt max {0, kxt k2 − α} . kxt k2 Ratio 12.5% 6.25% 25% 6.25% Önerilen 31.42 21.43 47.52 33.98 FCSA-MT 22.71 14.73 44.65 12.58 ayrık YDÇY 18.83 14.67 37.73 17.58 (b) Şekil 1: Monte-Carlo benzetim sonuçlarına göre pSNR değerlerinin ortalama ve varyans grafikleri görüntü kalitesi üzerindeki etkiyi göstermektedir. Bu hızlandırma faktöründe, görüntü kalitesi ve çözünürlük düşmekte ve görüntüdeki bozulmalar doku sınırlarını gölgelemektedir: ok ile gösterilmiş örnek bölgelerde doku ayrımları daha az belirgindir. Dahası, diğer yöntemlerde düşük çözünürlükten kaynaklı hatalar bulunmaktadır. Önerilen metod yalnızca farklı doku tipleri arasındaki kontrastı korumakla kalmayıp, arka planda bulunan gürültünün en aza indirgenmesinde de pSNR açısından (tablo 1) değerlendirildiğinde daha iyi bir performans elde etmiştir. V. TARTI ŞMA Bu çalışma YDÇY tabanlı hesapsal olarak verimli bir çoklu-kontrast görüntü gerikazanım problemi için bir algoritmayı sunmaktadır. Algoritma sıkıştırılmış çok-karşıklıklı MRG verisi üzerinde gösterilmiş ve hem ayrık YDÇY [7] hem de FCSA-MT [1] ile kıyaslanmıştır. YDÇY tabanlı algoritmanın özellikle çok az veri ile hem algoritmayı hızlandırdığı hem de pSNR açısından iyileşme sağladığı gösterilmiştir. Çoklukontrast görüntülere ihtiyaç duyulduğu durumlarda bilinen yöntemlerden çok daha kısa sürede görüntülemenin tamlamlanmasının önü açılmıştır. Bu yöntem çoklu-kontrast MRG için gösterilmiş olmasına rağmen, herhangi bir çok-kanallı/çoklukontrast gerikazanım problemi için kullanılabilir. K AYNAKÇA [1] [2] Huang, J., Chen, C., Axel, L. (2014). Fast multi-contrast MRI reconstruction. Magnetic Resonance Imaging, 32(10), 1344–1352. Lustig, M., Donoho, D., Pauly, J.M. (2007). Sparse MRI: The application of compressed sensing for rapid MR imaging. Magnetic Resonance in Medicine, 58(6),1189–1195. Önerilen (a) FCSA-MT ayrık YDÇY Data BS BS SR SR Referans Tablo I: dB cinsinden ortalama pSNR değerleri Şekil 2: Soldan Sağa Kanallar: Proton Yoğunluğu (PY), T1ağırlıklı (T1a), T2-ağırlıklı (T2a) görüntüleri, tüm verinin %12.5’i bulunan BS görüntüleri [3] [4] [5] [6] [7] [8] [9] [10] [11] Cukur, T., Lustig, M., Saritas, E. U., Nishimura, D.G. (2011). Signal Compensation and Compressed Sensing for Magnetization-Prepared MR Angiography. IEEE Transactions on Medical Imaging, 30(5), 1017–1027. Boyd, S., Parikh, N., Chu, E., Peleato, B., Eckstein, J. (2011). Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers. Found. Trends Mach. Learn., 3(1), 1–122. Afonso, M. V., Bioucas-Dias, J. M., Figueiredo, M. A. T. (2011). An Augmented Lagrangian Approach to the Constrained Optimization Formulation of Imaging Inverse Problems. Trans. Img. Proc., 20(3), 681– 695. Guven, H. E., Cetin, M. (2014). An Augmented Lagrangian method for sparse SAR imaging. Paper presented at Synthetic Aperture Radar (EUSAR 2014), 10th European Conference on. Guven, H. E., Gungor, A., Cetin, M. (2015, 27-30 Sept. 2015). An Augmented Lagrangian Method for image reconstruction with multiple features. Paper presented at the Image Processing (ICIP), 2015 IEEE International Conference on. Bresson, X., Chan, T. F. (2008). Fast Dual Minimization of The Vectorial Total Variation Norm and Applications to Color Image Processing. Inverse Problems and Imaging, 2(4), 455–484. Aujol, J.-F. (2009). Some First-Order Algorithms for Total Variation Based Image Restoration. Journal of Mathematical Imaging and Vision, 34(3), 307–327. Cetin, M., Karl, W. C. (2001). Feature-enhanced synthetic aperture radar image formation based on nonquadratic regularization. Image Processing, IEEE Transactions on, 10(4), 623–631. Chambolle, A. (2005). Total variation minimization and a class of binary MRF models. Paper presented at the Proceedings of the 5th international conference on Energy Minimization Methods in Computer Vision and Pattern Recognition, St. Augustine, FL.