real-time EDITER

seq/a2re.py

@@ -30,6 +30,7 @@ def __init__(self):
         self.addParameter(key='shimming', string='线性匀场 [x, y, z]', val=[210.0, 210.0, 895.0], field='OTH')
         self.addParameter(key='gFactor', string='梯度效率 [x, y, z]', val=[1.0, 1.0, 1.0], field='OTH')
         self.addParameter(key='readPadding', string='读出边距 (μs)', val=1.0, field='OTH')
+        self.addParameter(key='emiCh', string='噪声通道', val=[2, 3], field='OTH')
 
     def sequenceInfo(self):
         print("============ Averaging Acquisition with Relaxation Enhancement ============")
@@ -182,23 +183,49 @@ def sequenceAnalysis(self):
             ksp = self.mapVals['ksp3d_ch%d' % ch] = np.mean(self.mapVals['mse3d_ch%d' % ch], 2)  # 对回波链取平均
             self.mapVals['img3d_ch%d' % ch] = ifftshift(ifftn(ifftshift(ksp)))  # 重建
 
-        img = self.mapVals['img3d_ch0']    
-        abs_img = np.abs(img)
+        ksp_editer = [self.editer(
+            self.mapVals['ksp3d_ch0'][layer],
+            [self.mapVals['ksp3d_ch%d' % ch][layer] for ch in self.emiCh]
+        ) for layer in range(self.nPoints[2])]
+        img_editer = ifftshift(ifftn(ifftshift(ksp_editer)))
+        self.mapVals['img3d_editer'] = img_editer
+
+        img = self.mapVals['img3d_ch0']
+        img_result = np.array([img, img_editer]).transpose(1, 0, 2, 3).reshape(self.nPoints[2], -1, self.nPoints[0])
+        print('图像尺寸: ', img_result.shape)
         self.output = self.out = [{
             'widget': 'image',
-            'data': np.concatenate((abs_img/np.max(abs_img)*2 - 1, np.angle(img) / np.pi)),
+            'data': np.abs(img_result),
             'xLabel': '相位编码',
             'yLabel': '频率编码',
-            'title': '幅值图与相位图',
+            'title': '幅值图',
             'row': 0,
             'col': 0
+        }, {
+            'widget': 'image',
+            'data': np.angle(img_result),
+            'xLabel': '相位编码',
+            'yLabel': '频率编码',
+            'title': '相位图',
+            'row': 0,
+            'col': 1
         }, {
             'widget': 'image',
             'data': np.log(np.abs(self.mapVals['ksp3d_ch0'])),
             'xLabel': 'mV',
             'yLabel': 'ms',
             'title': 'k-Space',
-            'row': 0,
+            'row': 1,
+            'col': 0
+        }, {
+            'widget': 'image',
+            'data': np.abs([
+                self.mapVals['emi3d_ch%d' % ch] for ch in self.emiCh
+            ]).reshape(-1, self.nPoints[1], self.nPoints[0]),
+            'xLabel': 'mV',
+            'yLabel': 'ms',
+            'title': '噪声',
+            'row': 1,
             'col': 1
         }]
         self.saveRawData()

seq/mriBlankSeq.py

@@ -1722,3 +1722,83 @@ def autoProcessing(self, sampled, k_space):
         image = self.runIFFT(k_sp_zp)
 
         return image
+
+    def editer(ksp: np.ndarray, emi: list, ksz_col=0, ksz_lin=0):
+        def get_noise_kernel(data: list, pe_rng: int, ksz_col: int, ksz_lin: int):
+            noise_mat = []
+            for dat in data:
+                tmp = dat[:, pe_rng]
+                if isinstance(pe_rng, int):
+                    tmp = [tmp]
+                dfp = np.pad(tmp, ((ksz_col, ksz_col), (ksz_lin, ksz_lin)))
+                for col_shift in range(-ksz_col, ksz_col + 1):
+                    for lin_shift in range(-ksz_lin, ksz_lin + 1):
+                        dftmp = np.roll(dfp, (col_shift, lin_shift), (0, 1))
+                        noise_mat.append(dftmp[
+                            ksz_col:dftmp.shape[0] - ksz_col, 
+                            ksz_lin:dftmp.shape[1] - ksz_lin
+                        ])
+            noise_mat = np.array(noise_mat)
+            # kernels
+            gmat = noise_mat.reshape(noise_mat.shape[0], -1).T
+            init_mat_sub = ksp[:, pe_rng]
+            kernel, _, _, _ = np.linalg.lstsq(gmat, init_mat_sub.flatten(), None)
+            return kernel, gmat, init_mat_sub
+
+        # image size
+        ncol = ksp.shape[0]
+        nlin = ksp.shape[1]
+        Nc = len(emi)
+
+        # Initial pass using single PE line ( Nw =1 )
+        ksz_col = 0  # Deltakx = 1
+        ksz_lin = 0  # Deltaky = 1
+
+        # kernels across pe lines
+        kern_pe = np.zeros((Nc * (2*ksz_col+1) * (2*ksz_lin+1), nlin), np.complex64)
+
+        for clin in range(nlin):
+            kern_pe[:, clin], _, _ = get_noise_kernel(emi, clin, ksz_col, ksz_lin)
+
+        # look at correlation between the kernels
+        kern_pe_normalized = kern_pe / np.linalg.norm(kern_pe, axis=0)
+
+        kcor = np.dot(kern_pe_normalized.T, kern_pe_normalized)
+
+        # threshold
+        kcor_thresh = np.abs(kcor) > 5e-1
+
+        # start with full set of lines
+        aval_lins = list(range(nlin))
+
+        # window stack
+        win_stack = [None] * (nlin + 1)
+
+        cwin = 0
+        while aval_lins:
+            clin = min(aval_lins)
+            pe_rng = list(range(clin, 1 + clin + np.max(np.where(kcor_thresh[clin, clin:]))))
+
+            win_stack[cwin] = pe_rng
+            aval_lins = sorted(list(set(aval_lins) - set(pe_rng)))
+            cwin += 1
+
+        # drop the empty entries
+        win_stack = win_stack[:cwin]
+
+        ksz_col = 7  # Deltakx
+        ksz_lin = 0  # Deltaky
+
+        # solution kspace
+        gksp = np.zeros((ncol, nlin), np.complex64)
+
+        for cwin in range(len(win_stack)):
+            pe_rng = win_stack[cwin]
+
+            kern, gmat, init_mat_sub = get_noise_kernel(emi, pe_rng, ksz_col, ksz_lin)
+
+            # put the solution back
+            tosub: np.ndarray = np.dot(gmat, kern)
+            gksp[:, pe_rng] = init_mat_sub - tosub.reshape(ncol, len(pe_rng))
+
+        return gksp