PET重建技术 MLEM迭代法（C++）（一） 原理及成像

PET重建MLEM迭代法（一） 原理及成像

 一、名词解释：

1）系统矩阵：
        系统矩阵，也叫概率矩阵（通常用阵A表示），在重建中属于一个已知量或是一个固定参数。它是一个稀疏矩阵（大多数元素为零的矩阵）。在PET重建中所具有的系统矩阵为每个点所接收到电子的概率，从而判定该点应该有多少放射电子。这里A(i,j)表示第j个像素发射的光子对被第i个探测器接收的概率。
    系统矩阵的存储格式为三元（线）表：（将非零元素的行、列、值通过二维数组存储）

     2）测试数据：
     从机器上获得的每个点所接受到的电子的数量，一般为一列科学技术法，通过代码转换为十进制数后转存为向量Y；

二、成像原理：

    根据PET的工作原理：

     实际的接收到的电子的数据（Y）=电子接受（发送）的概率矩阵（A）*人体内的放射性元素数量（向量X）；
  
   即：Y=A*X；

     而将X根据各个点的放射性元素多少来分配明暗程度即可得到图像；
 
     因此求得X的过程即为求得图像过程；

    由于计算量的问题按照上式进行X的求解并不现实，因此我们采用迭代的方法使一个随机的X值通过多次的迭代接近真正的X从而得到图像X的近似图像；

三、迭代公式（MLEM）：

      
   我们可以通过读取矩阵A向量Y并进行计算得到X的图片的值
四、C++代码的实现：

 1）矩阵的存储：通过建立Class SparseMtx 来存储概率矩阵 
                    建立Class Vector 来存储向量 
                    建立class ELEM 来存储每个非零的点的信息（行、列、值）
为了方便，将三个存储用的class 存储在了一个.h头文件中：

#pragma once
#include "stdafx.h"

#include <string>
#include <stdio.h>
using namespace std;

#ifndef EPSILON
#define EPSILON       0.00000000001
#endif

class ELEM
{
public:
int row;
int col;
double val;
ELEM()
{
row = 0;
col = 0;
val = 0;
}
ELEM(int i, int j, double d)
{
row = i;
col = j;
val = d;
}
void operator=(ELEM t)
{
row = t.row;
col = t.col;
val = t.val;
}
bool operator<(ELEM t)
{
if (row < t.row)
return true;
else if (row == t.row && col < t.col)
{
return true;
}
else
return false;
}

bool operator >(ELEM t)
{
if (row > t.row)
return true;
else if (row == t.row && col > t.col)
{
return true;
}
else
return false;
}

bool operator == (ELEM t)
{
if (row == t.row && col == t.col)
{
return true;
}
else
{
return false;
}
}

bool operator!=(ELEM t)
{
return !(*this == t);
}

bool operator>=(ELEM t)
{
return (*this > t || *this == t);
}

bool operator<=(ELEM t)
{
return (*this < t || *this == t);
}
};
class Vector
{
public:
int iNum;
double* pdData;
// 构造函数
Vector(int iDimen)
{
iNum = iDimen;
pdData = new double[iNum];
}

// 析构函数
~Vector()
{
//delete[] pdData;
}
void getvector(int i, double val)
{
pdData[i] = val;
}
void Initilize()
{
memset(pdData, 0, iNum * sizeof(double));
}

void Initilize(double dVal)
{
int i;
for (i = 0; i < iNum; i++)
{
pdData[i] = dVal;
}
}

double Sum()
{
double d = 0.0;
for (int i = 0; i < iNum; i++)
{
d += pdData[i];
}
return d;
}

void ShowVector()
{
CString add, message;
int i;
for (i = 0; i < iNum; i++)
{
printf("%f ", pdData[i]);

add.Format("%f/t", pdData[i]);
message = message + add;
add = "";
}
MessageBox(NULL, message, "提示", MB_OK);
}
int SaveVector(CString pchFileName)
{
FILE *fp;
errno_t err;
if ( (err = fopen_s(&fp,pchFileName, "wb"))!=0)
{
printf("Fail to open writen file\n");
return (-1);
}
for (int i = 0; i < iNum; i++)
fprintf(fp,"%lf\r\n", pdData[i]);
//{
//	printf("Error data size\n ");
//	return (-1);
//}
fclose(fp);
return 0;
}
};

class SparseMtx
{
//private:
public:
int n_rows;
int n_cols;
int n_max_ELEM;
int n_actual_ELEM;
ELEM *items;

public:
// 构造函数
SparseMtx(int rows, int cols, int n_max);

// 析构函数
~SparseMtx()
{
//free(items);
}
// 稀疏矩阵初始化
void Empty();

// 矩阵转置
void Transpose();

4000

// 添加元素，添加了处理溢出的机制
bool AddElment(ELEM e);

// 对元素按行进行一次排序
void ElemRowSort();

// 缩减稀疏矩阵存储空间到最低值
void CutSparMtxMem();

// 显示函数 for debugging
void ShowElements();

// 元素交换
void MySwap(ELEM &s1, ELEM &s2);

// 递归快速排序，注意递归算法处理大的数据容易栈溢出
void QuickSort(ELEM* vec, int low, int high);
// 保存到文件
void SaveToFile(char* pchFileName);
};

 下面是.cpp：

#include "stdafx.h"
#include "SparseMtx.h"
#include<iostream>
#include<vector>
#include<stack>
#include<cstdlib>
#include<algorithm>
#include <memory.h>

static bool is_debugging = false;
#ifndef SPAR_MTX_LEN
#define SPAR_MTX_LEN 2000000
#endif

SparseMtx::SparseMtx(int rows, int cols, int n_max)
{
this->n_rows = rows;
this->n_cols = cols;
this->n_max_ELEM = n_max;
this->n_actual_ELEM = 0;
this->items = (ELEM*)malloc(n_max_ELEM * sizeof(ELEM));
if (NULL == this->items)
{
printf("稀疏矩阵构造函数分配内存空间失败\n");
}
}

void SparseMtx::MySwap(ELEM &s1, ELEM &s2)
{
ELEM t = s1;
s1 = s2;
s2 = t;
}

void SparseMtx::QuickSort(ELEM* arr, int left, int right)
{
// key位置和对应的值
int key_pos = left;
ELEM key_value = arr[key_pos];

// 循环搜索范围
int i = left;
int j = right;

// 一次快速排序
while (i < j)
{
// 从左向右，大于key就互换
while (arr[i]<key_value && i<right)
i++;
if (i<j)
{
MySwap(arr[i], arr[key_pos]);
key_pos = i;
}

// 从右向左，小于key就互换
while (arr[j]>key_value && j>left)
j--;
if (i<j)
{
MySwap(arr[j], arr[key_pos]);
key_pos = j;
}
}

// 递归调用，key的左右两个子集进行一次快速排序
if (key_pos>left)
{
QuickSort(arr, left, key_pos - 1);
}
if (key_pos<right)
{
QuickSort(arr, key_pos + 1, right);
}
}

/************************************************************************/
/*      为了提高速度，放弃了对重复元素判断
/************************************************************************/
bool SparseMtx::AddElment(ELEM e)//进行矩阵读入的函数，输入类型为为ELEM,可将三线表中行列值读入ELEM中 在将每个点放入矩阵中
{
// check，防止溢出
if (e.col <0 || e.col>n_cols
|| e.row <0 || e.row>n_rows)
{
printf("OverFlow when adding elements...\n");
return false;
}
else if (this->n_actual_ELEM >= this->n_max_ELEM)// 需要处理溢出
{
this->n_max_ELEM += SPAR_MTX_LEN;
printf("追加分配空间，新的n_max_ELEM = %d\n", n_max_ELEM);
this->items = (ELEM*)realloc(this->items, n_max_ELEM*sizeof(ELEM));
if (NULL == this->items)
{
printf("追加分配空间失败\n");
exit(-1);
}
this->items[this->n_actual_ELEM] = e;
this->n_actual_ELEM++;

return true;
}
else
{
this->items[this->n_actual_ELEM] = e;
this->n_actual_ELEM++;
return true;
}

}

void SparseMtx::ShowElements()
{
for (int i = 0; i<this->n_actual_ELEM; i++)
{
printf("(%d, %d)\t%.4f\n", this->items[i].row, this->items[i].col, this->items[i].val);
CString message;
message.Format("%d %d %lf", this->items[i].row, this->items[i].col, this->items[i].val);
MessageBox(NULL, message, "检测", MB_OK);
}
}

// 分成两部分：首先把所有点按行排列，再在每行中进行排序
void SparseMtx::ElemRowSort()
{
int i, j;
int iRowNum;
int iColNum;
int iTotElemNum;
ELEM e;
iRowNum = this->n_rows;
iColNum = this->n_cols;
iTotElemNum = this->n_actual_ELEM;

int* piRowElemNum = new int[iRowNum];
int* piStartPos = new int[iRowNum];
int* piStartPosCopy = new int[iRowNum];
double* pdMem = new double[iColNum];

ELEM* NewItem = new ELEM[iTotElemNum];

if (iTotElemNum > 0)
{
memset(piRowElemNum, 0, iRowNum * sizeof(int));
memset(piStartPos, 0, iRowNum * sizeof(int));
}
else
{
printf("数组内没有元素!\n");
return;
}

for (i = 0; i < iTotElemNum; i++)
{
piRowElemNum[this->items[i].row]++; //统计第行元素的个数
}
/**//*处理存放地址*/
piStartPos[0] = 0;
for (i = 1; i < iRowNum; i++)
{
piStartPos[i] = piStartPos[i - 1] + piRowElemNum[i - 1];
}

// 首先按行粗排序，使得同行的数据相连
memcpy(piStartPosCopy, piStartPos, iRowNum * sizeof(int));

for (i = 0; i < iTotElemNum; i++)
{
j = piStartPos[this->items[i].row];
NewItem[j] = this->items[i];
piStartPos[this->items[i].row]++;
}
// 对每一行下的所有元素进行排序

this->Empty();
for (i = 0; i < iRowNum; i++)
{
memset(pdMem, 0, iColNum * sizeof(double));
if (0 == piRowElemNum[i])
{
continue;
}
// 重排
for (j = piStartPosCopy[i]; j < piStartPosCopy[i] + piRowElemNum[i]; j++)
{
pdMem[NewItem[j].col] = NewItem[j].val;
}
for (j = 0; j < iColNum; j++)
{
if (pdMem[j] > EPSILON)
{
e.row = i;
e.col = j;
e.val = pdMem[j];
this->AddElment(e);
}
}
}

delete[] pdMem;
delete[] piRowElemNum;
delete[] piStartPos;
delete[] piStartPosCopy;
delete NewItem;

}

void SparseMtx::SaveToFile(char* pchFileName)
{
int i;
FILE* fp;
errno_t err;
err=fopen_s(&fp,pchFileName, "w+");
if (err!=0)
{
//prrintf("Cannot open file %s\n", pchFileName);
SparseMtx A(16354, 16354, 1);
CString message;
message.Format("Cannot open file %s!", pchFileName);
AfxMessageBox(message);
return;
}

for (i = 0; i< this->n_actual_ELEM; i++)
{
fprintf(fp, "%d\t%d\t%f\n", this->items[i].row, this->items[i].col, this->items[i].val);
}
fclose(fp);
}

void SparseMtx::Empty()
{
for (int i = 0; i < this->n_max_ELEM; i++)
{
this->items[i].col = -1;
this->items[i].row = -1;
this->items[i].val = -1;
}
this->n_actual_ELEM = 0;
}

void SparseMtx::CutSparMtxMem()
{
this->items = (ELEM*)realloc(items, this->n_actual_ELEM*sizeof(ELEM));
}

/**//*矩阵转置算法*/
void SparseMtx::Transpose()
{
int i, j;
int iRowNum;
int iColNum;
int iTotElemNum;
iRowNum = this->n_rows;
iColNum = this->n_cols;
iTotElemNum = this->n_actual_ELEM;

int* piColElemNum = new int[iColNum];
int* piStartPos = new int[iColNum];
ELEM* NewItem = new ELEM[iTotElemNum];

if (iTotElemNum > 0)
{
memset(piColElemNum, 0, iColNum * sizeof(int));
memset(piStartPos, 0, iColNum * sizeof(int));
}
else
{
printf("数组内没有元素!\n");
return;
}

for (i = 0; i < iTotElemNum; i++)
{
piColElemNum[this->items[i].col]++; //统计第i列的所有元素个数
}
/**//*处理存放地址*/
piStartPos[0] = 0;
for (i = 1; i < iColNum; i++)
{
piStartPos[i] = piStartPos[i - 1] + piColElemNum[i - 1];
}
/**//*存放*/
for (i = 0; i < iTotElemNum; i++)
{
j = piStartPos[this->items[i].col]; // 使得j能下一个地址
NewItem[j].row = this->items[i].col;
NewItem[j].col = this->items[i].row;
NewItem[j].val = this->items[i].val;
piStartPos[this->items[i].col]++;
}
for (i = 0; i < iTotElemNum; i++)
{
this->items[i] = NewItem[i];
}

i = this->n_rows;
this->n_rows = this->n_cols;
this->n_cols = i;

delete[] piColElemNum;
delete[] piStartPos;
delete NewItem;
}

其中Class SparseMtx 定义变量时需要输入矩阵的行列数 Vector 需要输入向量元素个数；
 
  2）文件的读入：

矩阵的读入：
 
 通过定义FILE 和使用 fopen函数 fscanf函数（加快读入速度）来进行文件的读入（读取Y时需要将科学计数法转化为十进制数字进行读入；

读入矩阵的代码：

double val;
FILE*fp;
errno_t err;
if ((err = fopen_s(&fp, a, "r")) != 0)
{
CString message;
message.Format("Can't open the file: %s", a);
AfxMessageBox(message);
}
else
{
int row,col;
while (fscanf(fp, "%d%d%lf", &row, &col, &val) != EOF)
{
//printf("%d %d %f /n");

ELEM e(row, col, val);
A.AddElment(e);
//			A.ShowElements();
}
}
//A.ElemRowSort();
fclose(fp);

读取向量的代码：

if (NULL == (fp = fopen(y, "r")))
{
CString message;
message.Format("Can't open the file: %s", y);
AfxMessageBox(message);
}
else
{
int k = 0;
while (fscanf(fp,"%s",&data)!=EOF)
{
if ((data[0]=='-'&&data[9]!='e')||(data[0]!='-'&&data[8]!='e'))
{
CString message;
message.Format("%s", data);

MessageBox(NULL, message, "测试", MB_OK);
}
//科学计数法转换：
if (data[0] == '\n')
break;
val = 0;
bool x = true;
int d;
double a=1.0;
if (data[0] == '-')
{
x = false;
for (int i = 0; i <= 13; i++)
data[i] = data[i + 1];
}
d = 100 * (data[10]-'0') + 10 * (data[11]-'0') + (data[12]-'0');
if (data[9] == '+')
{
for (int i = 0; i < d; i++)
a = a * 10;
}
else if(data[9]=='-')
{
for (int i = 0; i < d; i++)
a = a/10;
}
else MessageBox(NULL, "出错1", "测试", MB_OK);
for(int i = 0; i <= 7; i++)
{
if (i != 1)
{
double min = 1.0;
for (int l = 2; l <= i; l++)
min = min / 10;
val = val + (data[i]-'0')*min;
}
}
if (x == true)
val = val*a;
else if(x==false)
{
val = -1 * val*a;
}
else
{
MessageBox(NULL, "出错2", "测试", MB_OK);
}
//将转化完的科学计数法进行读入：
Y.getvector(k, val);
CString message;
message.Format("%lf", Y.pdData[k]);
//if (x== false)
//MessageBox(NULL, message, "测试", MB_OK);
k++;

}

}
fclose(fp);

 3）矩阵与向量计算：
 通过一个Class  SprsMtxOperator 将所有需要用到的 矩阵计算加入进去：

头文件：

#pragma once
#include "SparseMtx.h"
#include "stdafx.h"

#ifndef EPSILON
#define EPSILON       0.000000001
#endif
class SprsMtxOperator
{
public:
SprsMtxOperator();
~SprsMtxOperator();

// 矩阵与向量相乘
int MtxVecMultiple(const SparseMtx& Mtx, Vector& InVec, Vector& ResultVec);

// 向量与矩阵相乘
int VecMtxMultiple(Vector& InVec, const SparseMtx& Mtx, Vector& ResultVec);

// 向量逐项相乘
int VecMultiple(Vector& FirstVec, Vector& SecondVec, Vector& ResultVec);

// 向量逐项相除
int VecDiv(Vector& FirstVec, Vector& SecondVec, Vector& ResultVec);

// 向量的每一项乘以一个因子
int VecScale(Vector& FirstVec, double dScale, Vector& ResultVec);

// 计算稀疏矩阵乘积的列和
int MtxColSum(const SparseMtx& Mtx, Vector& ResultVec);

};

#pragma once

#include "stdafx.h"
#include <stdio.h>
#include "SprsMtxOperator.h"
using namespace std;
SprsMtxOperator::SprsMtxOperator(void)
{
}

SprsMtxOperator::~SprsMtxOperator(void)
{
}
// 矩阵乘以向量
int SprsMtxOperator::MtxVecMultiple(const SparseMtx& Mtx, Vector& InVec, Vector& ResultVec)
{
int i;
int iRowNum;
int iColNum;
int iTotElemNum;
int iVecNum;

iRowNum = Mtx.n_rows;
iColNum = Mtx.n_cols;
iTotElemNum = Mtx.n_actual_ELEM;
iVecNum = InVec.iNum;

if (iColNum != iVecNum || iRowNum != ResultVec.iNum)
{
printf("Dimensions of matrix and vector mismatch\n");
return (-1);
}

Vector NewVector(iRowNum);
NewVector.Initilize();

for (i = 0; i < iTotElemNum; i++)
{
NewVector.pdData[Mtx.items[i].row] += Mtx.items[i].val * InVec.pdData[Mtx.items[i].col];
}

for (i = 0; i < ResultVec.iNum; i++)
{
ResultVec.pdData[i] = NewVector.pdData[i];
}
return 0;
}

// 向量乘以矩阵,该乘法也可是使用矩阵转置和MtxVecMultiple函数实现，
// 但是考虑到系统矩阵太大，矩阵转置消耗内存，因此不使用
int SprsMtxOperator::VecMtxMultiple(Vector& InVec, const SparseMtx& Mtx, Vector& ResultVec)
{
int i;
int iRowNum;
int iColNum;
int iTotElemNum;
int iVecNum;

iRowNum = Mtx.n_rows;
iColNum = Mtx.n_cols;
iTotElemNum = Mtx.n_actual_ELEM;
iVecNum = InVec.iNum;

if (iRowNum != iVecNum || iColNum != ResultVec.iNum)
{
printf("Dimensions of matrix and vector mismatch\n");
return (-1);
}

Vector NewVector(iColNum);
NewVector.Initilize();

for (i = 0; i < iTotElemNum; i++)
{
NewVector.pdData[Mtx.items[i].col] += Mtx.items[i].val * InVec.pdData[Mtx.items[i].row];
}

for (i = 0; i < ResultVec.iNum; i++)
{
ResultVec.pdData[i] = NewVector.pdData[i];
}
return 0;
}

int SprsMtxOperator::VecMultiple(Vector& FirstVec, Vector& SecondVec, Vector& ResultVec)
{
int i, iFirstNum, iSecondNum;
iFirstNum = FirstVec.iNum;
iSecondNum = SecondVec.iNum;

if (ResultVec.iNum != iFirstNum || ResultVec.iNum != iSecondNum)
{
printf("Dimensions of two vectors mismatch\n");
return (-1);
}
Vector NewVec(iFirstNum);
NewVec.Initilize();
for (i = 0; i < iFirstNum; i++)
{
NewVec.pdData[i] = FirstVec.pdData[i] * SecondVec.pdData[i];
}
for (i = 0; i < ResultVec.iNum; i++)
{
ResultVec.pdData[i] = NewVec.pdData[i];
}
return 0;
}

int SprsMtxOperator::VecDiv(Vector& FirstVec, Vector& SecondVec, Vector& ResultVec)
{
int i, iFirstNum, iSecondNum;
iFirstNum = FirstVec.iNum;
iSecondNum = SecondVec.iNum;

if (iFirstNum != iSecondNum)
{
printf("Dimensions of two vectors mismatch\n");
return (-1);
}
Vector NewVec(iFirstNum);
NewVec.Initilize();
for (i = 0; i < iFirstNum; i++)
{
NewVec.pdData[i] = FirstVec.pdData[i] / (SecondVec.pdData[i] + EPSILON);
}
for (i = 0; i < ResultVec.iNum; i++)
{
ResultVec.pdData[i] = NewVec.pdData[i];
}
return 0;
}

// 向量与一个数相乘
int SprsMtxOperator::VecScale(Vector& InVec, double dScale, Vector& ResultVec)
{
int i, iNum;
iNum = InVec.iNum;
if (ResultVec.iNum != iNum)
{
printf("Dim of output does not equal to input\n");
return (-1);
}

Vector NewVec(iNum);
NewVec.Initilize();
for (i = 0; i < iNum; i++)
{
NewVec.pdData[i] = InVec.pdData[i] * dScale;
}
for (i = 0; i < iNum; i++)
{
ResultVec.pdData[i] = NewVec.pdData[i];
}
return 0;
}

// 计算稀疏矩阵乘积的列和
int SprsMtxOperator::MtxColSum(const SparseMtx& Mtx, Vector& ResultVec)
{
int i;
int iRowNum;
int iColNum;
int iTotElemNum;
iRowNum = Mtx.n_rows;
iColNum = Mtx.n_cols;
iTotElemNum = Mtx.n_actual_ELEM;
if (iColNum != ResultVec.iNum)
{
printf("Result vector has error length\n");
return (-1);
}
ResultVec.Initilize();
for (i = 0; i < iTotElemNum; i++)
{
ResultVec.pdData[Mtx.items[i].col] += Mtx.items[i].val;
}
return 0;
}

4）迭代公式的 代码编写：

通过代码将迭代公式完成  并可以进行迭代：

int i = 0;
Vector X1(16384);
Vector X2(16384);
Vector X3(16384);
Vector X4(16384);
Vector X5(16384);
Vector X6(16384);
GetBest *best;
best = new GetBest;

SparseMtx Aa(16384,16384,1);
for (int l = 0; l<A.n_actual_ELEM; l++)
Aa.AddElment(A.items[l]);
Aa.Transpose();
for (i = 0; i < t; i++)
{

SprsMtxOperator * handle;
handle = new SprsMtxOperator;
handle->MtxVecMultiple(A, X, X1);
handle->VecDiv(Y, X1, X2);
handle->MtxVecMultiple(Aa, X2, X3);
handle->VecMultiple(X, X3, X4);
handle->MtxColSum(A, X5);
handle->VecDiv(X4, X5, X6);
CString filename;
filename.Format("%d.txt", i);
X6.SaveVector(filename);
for (int l = 0; l < X6.iNum; l++)
X.pdData[l] = X6.pdData[l];
filename.Format("%d.jpg", i);
GetPic *Pic;
Pic = new GetPic;
Pic->GetSavePic(X, filename);
i
9758
f(best->getbest==false)
best->GetBestPicA(Y, X, A, i);

}

使用公式可以进行迭代  并通过循环  更改t的 数量进行迭代

5）将迭代后的X向量转换为jpg图片文本：
迭代后  X依然是一个向量 可以将其元素个数开方 得到a，将其看作一个axa的矩阵 通过CImage 和 SetPixel描点公式进行描点 将像素值换算成一个个点画到图片上
成像代码：

double  max = 0;
int point = 0;

for (int i = 0; i < X.iNum; i++)
if (X.pdData[i] >= max)
max = X.pdData[i];
COLORREF COLOR = RGB(40, 123, 90);
CImage image;
image.Create(128, 128, 24);
HDC DC;
DC = image.GetDC();

for (int i = 1; i <=128; i++)
{
for (int l = 1; l <=128; l++)
{
point = (X.pdData[(i - 1) * 128 + l - 1] / max) * 255;
SetPixel(DC, i - 1, l - 1, RGB(point, point, point));
}
}
image.ReleaseDC();
DC = NULL;

image.Save(Picname);

这样就得到了  我们迭代出来的图像！！！

那么得到的图像通过多次迭代可以得到近似图像，迭代次数越多越接近原始图像，但是我们却不能盲目的加多迭代次数，过多的迭代会导致图片质量下降，那么到底应该进行多少次迭代呢？
请关注下次更新：
 
PET重建技术 MLEM迭代法（C++）（二） 迭代次数判定算法

（未完待续）

本文版权归作者和CSDN共有，欢迎转载，但未经作者同意必须保留此段声明，且在文章页面明显位置给出原文连接，否则保留追究法律责任的权利。