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#endifclass 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;}elsereturn false;}bool operator >(ELEM t){if (row > t.row)return true;else if (row == t.row && col > t.col){return true;}elsereturn 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();// 添加元素，添加了处理溢出的机制bool AddElment(ELEM e);// 对元素按行进行一次排序void ElemRowSort();// 缩减稀疏矩阵存储空间到最低值void CutSparMtxMem();// 显示函数 for debuggingvoid 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#endifSparseMtx::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#endifclass 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);if(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共有，欢迎转载，但未经作者同意必须保留此段声明，且在文章页面明显位置给出原文连接，否则保留追究法律责任的权利。