振幅、相位,重建MRI图像

# -*- coding: utf-8 -*-
"""
Created on Wed Oct 16 12:42:27 2019

@author: cyy
"""

import numpy as np
import matplotlib.pyplot as plt
import os
import pydicom


dicomfile_ge = pydicom.read_file('D:/program/read_dicm/IM-0001-0082.dcm')
img_original = dicomfile_ge.pixel_array
img_original = img_original.astype(np.float64)
f = np.fft.fft2(img_original)    # 快速傅里叶变换算法得到频率分布
img_R=np.fft.ifft2(f)
a_fshift = np.abs(f)  # 取振幅
ph_fshift = np.angle(f)  # 取相位


plt.figure()
plt.subplot(151), plt.imshow(img_original, 'gray'), plt.title('original')
plt.subplot(152), plt.imshow(abs(f),'gray'), plt.title('k')
plt.subplot(153),plt.imshow(abs(img_R),'gray'),plt.title('img_R')
plt.subplot(154), plt.imshow(a_fshift, 'gray'), plt.title('amp')
plt.subplot(155), plt.imshow(ph_fshift, 'gray'), plt.title('phase')
plt.show()




# 逆变换--利用欧拉公式 cos(x) + isin(x) ⇔ a + ib
# 将振幅和相位整合成复数的实部和虚部,即整合成频域数据
s_real = a_fshift*np.cos(ph_fshift)  # 整合复数的实部
s_imag = a_fshift*np.sin(ph_fshift)  # 整合复数的虚部

s = np.zeros([256, 256], dtype=complex)  # 指定数据类型为复数
s.real = np.array(s_real)
s.imag = np.array(s_imag)
#s = s.real + s.imag*1j
 # 对整合好的频域数据进行逆变换
#fshift1 = np.fft.ifftshift(s) 
#img_back = np.fft.ifft2(fshift1)
img_back = np.fft.ifft2(s)
# 出来的是复数,取绝对值,转化成实数
img_back = np.abs(img_back)
plt.figure()
plt.subplot(141), plt.imshow(img_original, 'gray'), plt.title('original')
plt.subplot(142), plt.imshow(a_fshift, 'gray'), plt.title('amp')
plt.subplot(143), plt.imshow(ph_fshift, 'gray'), plt.title('phase')
plt.subplot(144), plt.imshow(img_back, 'gray'), plt.title('inverse fourier')
plt.show()
'''
#img_original = img_original.astype(np.float64)
f = np.fft.fft2(img_original)    # 快速傅里叶变换算法得到频率分布
fshift = np.fft.fftshift(f)     # 移频,默认结果中心点位置是在左上角,转移到中间位置

img_R=np.fft.ifft2(fshift)    
                             
a_fshift = np.log(np.abs(fshift))  # 取振幅
ph_fshift = np.angle(fshift)  # 取相位

plt.figure()

plt.subplot(151), plt.imshow(img_original, 'gray'), plt.title('original')
plt.subplot(152), plt.imshow(abs(fshift),'gray'), plt.title('k')
plt.subplot(153), plt.imshow(abs(img_R),'gray'), plt.title('img_R')
plt.subplot(154), plt.imshow(a_fshift, 'gray'), plt.title('amp')
plt.subplot(155), plt.imshow(ph_fshift, 'gray'), plt.title('phase')
plt.show()


# 逆变换--利用欧拉公式 cos(x) + isin(x) ⇔ a + ib
# 将振幅和相位整合成复数的实部和虚部,即整合成频域数据
s_real = a_fshift*np.cos(ph_fshift)  # 整合复数的实部
s_imag = a_fshift*np.sin(ph_fshift)  # 整合复数的虚部

s = np.zeros([256, 256], dtype=complex)  # 指定数据类型为复数
s.real = np.array(s_real)
s.imag = np.array(s_imag)
#s = s.real + s.imag*1j
 # 对整合好的频域数据进行逆变换
#fshift1 = np.fft.ifftshift(s) 
#img_back = np.fft.ifft2(fshift1)
img_back = np.fft.ifft2(s)
# 出来的是复数,取绝对值,转化成实数
img_back = np.abs(img_back)
plt.figure()
plt.subplot(141), plt.imshow(img_original, 'gray'), plt.title('original')
plt.subplot(142), plt.imshow(a_fshift, 'gray'), plt.title('amp')
plt.subplot(143), plt.imshow(ph_fshift, 'gray'), plt.title('phase')
plt.subplot(144), plt.imshow(img_back, 'gray'), plt.title('inverse fourier')
plt.show()




def make_transform_matrix(d, image):
    transform_matrix = np.zeros(image.shape)
    center_point = tuple(map(lambda x: (x-1)/2, image.shape))
    for i in range(transform_matrix.shape[0]):
        for j in range(transform_matrix.shape[1]):
            def cal_distance(pa, pb):
                from math import sqrt
                # d值决定中心圆大小,遍历并计算频率域中的每一个点与这个圆边界的距离
                dis = sqrt((pa[0]-pb[0])**2+(pa[1]-pb[1])**2)
                return dis
            dis = cal_distance(center_point, (i, j))
            if dis <= d:
                transform_matrix[i, j] = 1
            else:
                transform_matrix[i, j] = 0
    return transform_matrix

d_1 = make_transform_matrix(40, a_fshift)  # 40 是中心园的半径
img_d1 = np.abs(np.fft.ifft2(np.fft.ifftshift(fshift*d_1)))
plt.subplot(141), plt.imshow(dicomfile_ge.pixel_array, 'gray'), plt.title('original')
plt.subplot(142), plt.imshow(a_fshift, 'gray'), plt.title('k-space amp')
plt.subplot(143), plt.imshow(d_1, 'gray'), plt.title('filter')
plt.subplot(144), plt.imshow(img_d1, 'gray'), plt.title('result')
plt.show()

'''





参考博客:https://www.rithia.net/2018/11/15/20181115/

©️2020 CSDN 皮肤主题: 大白 设计师:CSDN官方博客 返回首页