www.pudn.com > snake.rar > SnakeProcess.cpp
// SnakeProcess.cpp: implementation of the SnakeProcess class.
//
//////////////////////////////////////////////////////////////////////
#include "stdafx.h"
#include "mammocad.h"
#include "SnakeProcess.h"
#include "ROI.h"
#include "imageprocess.h"
#include
#include
#include "Matrix.h"
#include
#include
#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
extern void gradient(USHORT*f , // image buffer of growth region defined by active contour
int rows, int cols);
SnakeProcess::SnakeProcess()
{
}
SnakeProcess::~SnakeProcess()
{
}
// Snake.cpp : Defines the class behaviors for the application.
//
/////////////////////////////////////////////////////////////////////////////
// The one and only CSnakeApp object
LONG glHeight,glWidth;
bool glFlag;
struct DIRECTION
{
bool direction[8];
DIRECTION()
{
for(int i=0;i<8;i++)
direction[i]=FALSE;
}
};
struct Ppoint
{
int x,y;
Ppoint(int a=-1,int b=-1)
{
x=a;
y=b;
}
bool operator==(Ppoint tempPoint)
{
return (x==tempPoint.x && y == tempPoint.y);
}
Ppoint& operator=(Ppoint& tempPoint)
{
x=tempPoint.x;
y=tempPoint.y;
return *this;
}
};
/* NOTE:
*
* You must have "mex" command working in your matlab. In matlab,
* type "mex gvfc.c" to compile the code. See usage for details of
* calling this function.
*
* Replace GVF(...) with GVFC in the snakedeform.m. The speed is
* significantly faster than the Matlab version.
*/
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "string.h"
#define PG_MAX(a,b) (a>b? a: b)
#define PG_MIN(a,b) (amax_a)
{
max_a=pData_a[k*nCols+l];
}
//Para_b
if(pData_b[k*nCols+l]max_b)
{
max_b=pData_b[k*nCols+l];
}
//Para_c
if(pData_c[k*nCols+l]max_c)
{
max_c=pData_c[k*nCols+l];
}
//Para_d
if(pData_d[k*nCols+l]max_d)
{
max_d=pData_d[k*nCols+l];
}
}//End for(l=2;l big)
big = temp;
if(big == 0.0)
nrerror(" Singular matrix in routine ludcmp! ");
vv[i] = 1.0 / big;
}
for(j=1; j<=n; j++)
{
for(i=1; i= big)
{
big = dum;
imax = i;
}//求最大值位置
}
if(j != imax)
{
for(k=1; k<=n; k++)
{
dum = a[imax][k];
a[imax][k] = a[j][k];
a[j][k] = dum;
}
*d = -(*d);
vv[imax] = vv[j];
}
indx[j] = imax;
if(a[j][j] == 0.0)
a[j][j] = TINY;
if(j != n)
{
dum = 1.0 / a[j][j];
for(i=j+1; i<=n; i++)
a[i][j] *= dum;
}
}
free_dvector(vv,1,n+1);
return;
}
//和ludcmp共同实现逆矩阵
void lubksb(double **a, int n, int *indx, double b[])
{
int i, ii, ip, j;
double sum;
ii = 0;
for(i=1; i<=n; i++)
{
ip = indx[i];
sum = b[ip];
b[ip] = b[i];
if( ii )
{
for(j=ii; j<=i-1; j++)
sum -= (a[i][j] * b[j]);
}
else if( sum )
{
ii = i;
}
b[i] = sum;
}
for(i=n; i>=1; i--)
{
sum = b[i];
for(j=i+1; j<=n; j++)
sum -= (a[i][j] * b[j]);
b[i] = sum / a[i][i];
}
return;
}
//输入：轮廓点x、y坐标数组， 轮廓点数目 array_length；
//输出：轮廓点梯度图 diff；
//功能：求轮廓点序列的差分图；
void ComputeDifference(int array_length,
double* diff,
double* x, double* y
)//计算每个点的梯度－梯度图
{
int i, N;
double dx, dy;
double x_shift[24000], y_shift[24000];
// shift array position by one bit.
N = array_length - 1;
for(i=0; i 0) // 如果index > 0 ，那么IDX中存的是插值索引 ，此时：threshold=max_dist
{
for(i=0; i threshold)
IDX[i] = 1;
else
IDX[i] = 0;
}
}
else // 如果index <= 0 ，那么IDX中存的是删除多余点索引 ，此时：threshold=min_dist
{
for(i=0; i max_value)
max_value = diff_array[i];
}
return( max_value );
}
//输入：插值标识数组：IDX 轮廓点x或者y方向数组：array_i 轮廓点的数目：array_length ；
//输出：插值后的轮廓点x或y方向数组：array_d ；
//功能：线性插值，用相邻元素取均值进行插值；
int interp1(int* IDX, double* array_i, int array_length,
double* array_d
)
{
int i, k;
int new_length;
k = array_length - 1;
new_length = 0;
for(i=0; irows-1||k<1||k>cols-1)
{
pgFreeDmatrix(fx, 1, rows, 1, cols);
return ;
}
y1 = fx[j][k];
y2 = fx[j+1][k];
y3 = fx[j+1][k+1];
y4 = fx[j][k+1];
t = x1 - (double)(j);
u = x2 - (double)(k);
//双线性插值
a1 = (1.0 - t) * (1.0 - u) * y1;
a2 = t * (1.0 - u) * y2;
a3 = t * u * y3;
a4 = (1.0 - t) * u * y4;
vfx[i] = a1 + a2 + a3 + a4;
}
pgFreeDmatrix(fx, 1, rows, 1, cols);
return;
}
//输入：轮廓点x、y坐标数组 contour_x、contour_y , 轮廓点的数目 contour_length ，
// 最大距离和最小距离max_dist 和min_dist ；
//输出：插值之后的轮廓x、y坐标数组 contour_x、contour_y
//返回：轮廓点的数目；
//功能：snake演变之后插值；
int SnakeInterp(double* contour_x, double* contour_y, int contour_length,
double max_dist, double min_dist
)
{
int i, k, m, N;
int Interp_length;
double max_d;
double *diff;
int *IDX;
double *xi, *yi;
double *xo, *yo;
// convert to double format
N = contour_length;
m = 30 * N;
//申请空间；
diff = dvector(1, m);
xo = dvector(1, m);
xi = dvector(1, m);
yo = dvector(1, m);
yi = dvector(1, m);
IDX = ivector(1, m);
for(i=0; i max_dist ) //当轮廓点序列中的最大差分值比max_dist大时，需要继续插值；
{
if(Interp_length>3000)
{
max_d=0;
break;
}
DefineIDX(N, IDX, diff, 1, max_dist); //定义插值标识；
Interp_length = interp1(IDX, xo, N, xi); //x方向线性插值两次；
Interp_length = interp1(IDX, yo, N, yi); //y方向线性插值两次；
N = Interp_length;
CopyDDArray(N, xi, xo);
CopyDDArray(N, yi, yo);
ComputeDifference(N, diff, xo, yo); //计算差分图；
max_d = FindMaxValue(N, diff); //找差分图中最大值；
k ++;
}
CopyDDArray(N, xo, contour_x);
CopyDDArray(N, yo, contour_y);
contour_x[N] = 0;
contour_y[N] = 0;
contour_length = N;
//释放空间；
free_dvector(diff,1,m);
free_dvector(xo,1,m);
free_dvector(xi,1,m);
free_dvector(yo,1,m);
free_dvector(yi,1,m);
free_ivector(IDX,1,m);
return( contour_length );
}
//输出：array_x和array_y，其他参数均为输入；
//功能：snake形变
void SnakeDeform(double* array_x, double* array_y, int array_length, //轮廓点x、y数组以及轮廓的长度；
double alpha_single, double beta_single, //固定snake形变参数；
double gamma_single, double kappa_single, //固定snake形变参数
double* OutputU_buf, double* OutputV_buf, //U、V场
int rows, int cols, int ITER //图像的行和列，以及迭代的次数；
)//snake形变
{
int i, j, k, m, N;
int pixel_count;
int *indx;
double dv;
double *d_buf1, *d_buf2;
double *alpha, *beta, *gamma, *kappa;
double *alpham1, *alphap1;
double *betam1, *betap1;
double *a, *b, *c, *d, *e;
double *Avec, *Gvec;
double *col;
double *vfx, *vfy;
double *tmpX, *tmpY;
PGdouble **Amatrix;
PGdouble **InvAmat;
N = array_length;
m = 4 * N;
indx = ivector(1, m);
col = dvector(1, m);
alpha = dvector(1, m);
beta = dvector(1, m);
gamma = dvector(1, m);
kappa = dvector(1, m);
alpham1 = dvector(1, m);
alphap1 = dvector(1, m);
betam1 = dvector(1, m);
betap1 = dvector(1, m);
a = dvector(1, m);
b = dvector(1, m);
c = dvector(1, m);
d = dvector(1, m);
e = dvector(1, m);
tmpX = dvector(1, m);
tmpY = dvector(1, m);
vfx = dvector(1, m);
vfy = dvector(1, m);
Avec = dvector(N, N);
Gvec = dvector(N, N);
Amatrix = pgDmatrix(1, N, 1, N);
InvAmat = pgDmatrix(1, N, 1, N);///
for(i=0; i 0)
{
min_dd = dd;
min_i = i;
}
}
i = x_array[k+1];
j = y_array[k+1];
x_array[k+1] = x_array[min_i];
y_array[k+1] = y_array[min_i];
x_array[min_i] = i;
y_array[min_i] = j;
}
return( contour_count);
} // record_contour_array()
/*
************************************************************************************
* *
* This function is used to compute active contour *
* *
************************************************************************************
*/
int Computing_ActiveContour(USHORT* img_buf, // image buffer of original image
USHORT* growth_buf, // image buffer recorded the original growth region (binary image)
int rows, int cols, // rows and columns of image buffer
int orig_contour_count, // count of contour pixels
int* x_array, int* y_array // x, y array of contour pixels
)//计算GVF场并计算活动轮廓
{
USHORT *r_buf1, *r_buf2;
USHORT *EdgeMap_buf;
int i, j, k, n;
int ITER;
int pixel_count;
int contour_count;
double mu; //固定参数 mu=0.1;
double u, v, mag;
double max_dist, min_dist; //snakeInterp参数，相邻两点最大的距离max_dist和最小距离min_dist；
double alpha, beta, gamma, kappa; //固定sake形变参数；
double *FMap_buf;
double *OutputU_buf, *OutputV_buf;
double *px_buf, *py_buf;
double *contour_x, *contour_y;
double *f_buf1, *f_buf2, *f_buf3, *f_buf4;
//申请空间；
contour_count = orig_contour_count;
EdgeMap_buf = svector(rows, cols);
FMap_buf = dvector(rows, cols);
OutputU_buf = dvector(rows, cols);
OutputV_buf = dvector(rows, cols);
px_buf = dvector(rows, cols);
py_buf = dvector(rows, cols);
contour_x = dvector(rows*cols/3, 1);//??
contour_y = dvector(rows*cols/3, 1);
copy_svector(img_buf, EdgeMap_buf, rows, cols);
// read initial contour data array
for(k=0; kwidthStep;
BYTE * pData=(BYTE*)temp->imageData;
for(i=0;i<125;i++)
for(j=0;j<125;j++)
{
int jyx;
jyx=(*(EdgeMap_buf+125*i+j))/16;
if(jyx>255)
pData[i*stepsPL+j]=255;
else
pData[i*stepsPL+j]=jyx;
}
/*
cvNamedWindow("jyx",1);
cvShowImage("jyx",temp); //显示梯度图；
AfxMessageBox("nihao",MB_OK);*/
// gradient(EdgeMap_buf,rows,cols); //差分梯度；
// remove pixels inside growth region
pixel_count = rows * cols;
r_buf1 = growth_buf;
r_buf2 = EdgeMap_buf;
while( pixel_count-- )
{
if(*r_buf1 == MAX_12_BIT)
*r_buf2 = 0;
r_buf1 ++;
r_buf2 ++;
}
pixel_count = rows * cols;
f_buf1 = FMap_buf;
r_buf1 = EdgeMap_buf;
while( pixel_count -- )
{
*f_buf1 ++ = (double)(*r_buf1 ++);
}
// compute the GVF of the edge map
mu = 0.1;
ITER = 80;
GVFC(rows, cols, FMap_buf, OutputU_buf, OutputV_buf, mu, ITER);
// normalize the GVF external force.
f_buf1 = OutputU_buf;
f_buf2 = OutputV_buf;
f_buf3 = px_buf;
f_buf4 = py_buf;
pixel_count = rows * cols;
double minva,maxva;
minva=1.0;
maxva=-1.0;
while( pixel_count -- )
{
u = *f_buf1;
v = *f_buf2;
mag = sqrt(u * u + v * v) + 1e-10;
*f_buf3 = u / mag;
*f_buf4 = v / mag;
minva=(*f_buf3)maxva?(*f_buf3):maxva;
f_buf1 ++;
f_buf2 ++;
f_buf3 ++;
f_buf4 ++;
}
//测试；
for(i=0;i<125;i++)
for(j=0;j<125;j++)
{
pData[i*stepsPL+j]=(BYTE)(((double)255*((*(px_buf+125*i+j))-minva)/(maxva-minva))+0.5);
// pData[i*stepsPL+j]=(BYTE)((double)255*(*(px_buf+125*i+j)));
}
/*
cvNamedWindow("jyx",1);
cvShowImage("jyx",temp);
AfxMessageBox("nihao",MB_OK);*/
// snake deformation
max_dist = 2.0;
min_dist = 0.5;
k = SnakeInterp(contour_x, contour_y, contour_count, max_dist, min_dist); //插值
contour_count = k;
// iteration for snake deformation and interpretation.
ITER = 5;
alpha = 0.05;
beta = 0.0;
gamma = 1.0;
kappa = 0.6;
int snake_iter =20; //snack_iter =9；迭代次数多些没关系；
for(k=1; k1&&contour_count<3000)
{
SnakeDeform(contour_x, contour_y, contour_count, alpha, beta, gamma, kappa,
px_buf, py_buf, rows, cols, ITER);
j = SnakeInterp(contour_x, contour_y, contour_count, max_dist, min_dist);
contour_count = j;
}
}
max_dist=1.0;
min_dist=0.5;
if(contour_count<3000)
{
j = SnakeInterp(contour_x, contour_y, contour_count, max_dist, min_dist);
contour_count = j;
}
// record the detected contour cordinate (x, y) into x_array and y_array
pixel_count = 0;
for(k=0; k124||j>124)
{
continue;
}
n = *(growth_buf+j*cols+i);//用到ActiveContourRegion函数返回的growth_buf值
if(n == 0)
{
x_array[pixel_count] = (int)(0.5 + contour_x[k]);
y_array[pixel_count] = (int)(0.5 + contour_y[k]);
pixel_count ++;
}
}
contour_count = pixel_count;
x_array[contour_count] = 0;
y_array[contour_count] = 0;
//释放空间；
free_svector(EdgeMap_buf);
free_dvector(FMap_buf,rows,cols);
free_dvector(OutputU_buf,rows,cols);
free_dvector(OutputV_buf,rows,cols);
free_dvector(px_buf,rows,cols);
free_dvector(py_buf,rows,cols);
free_dvector(contour_x,rows*cols/3, 1);
free_dvector(contour_y,rows*cols/3, 1);
free_svector(Para_a);
free_svector(Para_b);
free_svector(Para_c);
free_svector(OrigImg);
cvReleaseImage(&temp);
return( contour_count );
} // End of Computing_ActiveContour()
//输入：rows，scol；contour_count；x_array，y_array
//输出：growth_buf
//功能：计算包含肿块区域的矩形，该矩形的功能实际上是限制梯度图和GVF场的范围；
void ActiveContourRegion(USHORT* growth_buf, // image buffer of growth region defined by active contour
int rows, int cols, // rows and columns of image buffer
int contour_count, // number of pixels in boundary contour
int* x_array, int* y_array // contour (x, y) arrays
)//将growth_buf赋值
{
USHORT *LocalGrowth_buf;
USHORT *LocalMorph_buf;
USHORT *r_buf1;
int i, j, k, pixel_count;
int contour_pixels;
int cent_x, cent_y;
int max_x, max_y, min_x, min_y;
int scol, srow;
double *contour_x, *contour_y;
contour_x = dvector(rows*cols/3, 1);//??
contour_y = dvector(rows*cols/3, 1);
for(k=0; k max_x)
max_x = i;
if(j < min_y)
min_y = j;
if(j > max_y)
max_y = j;
}
r_buf1 ++;
}
}
//将包含肿块的最小矩形向外扩充20个像素点，得到的新矩形区域为受限肿块区域；
if(contour_pixels > 0)
{
cent_x /= contour_pixels; // (cent_x,cent_y)为肿块的几何中心；
cent_y /= contour_pixels;
min_x -= 20;
if(min_x < 0)
min_x = 0;
max_x += 20;
if(max_x >= scol)
max_x = scol - 1;
min_y -= 20;
if(min_y < 0)
min_y = 0;
max_y += 20;
if(max_y >= srow)
max_y = srow - 1;
}
// 将受限肿块区域之外的区域像素点灰度值设为4095；受限肿块区域灰度值设为0；
for(j=0; jwidthStep;
BYTE* pFlagData = (BYTE*)flagImg->imageData;
CvPoint candi[4] = {0, 0, flagImg->width-1, 0, 0, flagImg->height-1, flagImg->width-1, flagImg->height-1};
for(i=0; i<4; i++)
if(pFlagData[candi[i].y*steps_PL+candi[i].x] != 255)
return candi[i];
}
//根据边界点修改区域标志：0---背景，178---区域内像素，255---边界像素
void Cal_RegionFlag(IplImage* flagImg)//flagImg: 8U
{
using namespace std;
IplImage* tempImg = cvCloneImage(flagImg);
int i, j;
int steps_PL = flagImg->widthStep;
BYTE* pTempData = (BYTE*)tempImg->imageData;
BYTE* pFlagData = (BYTE*)flagImg->imageData;
CvPoint seed = FindSeed(flagImg);
bool** visited = new bool*[flagImg->height];
for(i=0; iheight; i++)
visited[i] = new bool[flagImg->width];
for(i=0; iheight; i++)
for(j=0; jwidth; j++)
{
visited[i][j] = false;
}
//生长
list stack;
stack.push_back(seed);
pFlagData[seed.y* steps_PL + seed.x] = 0;
visited[seed.y][seed.x] = true;
int dir[4][2] = {{1, 0}, {0, 1}, {0, -1}, {-1, 0}};
while(!stack.empty())
{
CvPoint curPoint = stack.back();//取出
stack.pop_back();//弹出
for(i=0; i<4; i++)
{
CvPoint nextPoint = cvPoint(curPoint.x + dir[i][0], curPoint.y + dir[i][1]);
if(nextPoint.x < 0 || nextPoint.y <0
|| nextPoint.x > flagImg->width-1 || nextPoint.y > flagImg->height-1)
continue;
if(visited[nextPoint.y][nextPoint.x])
continue;
if(pTempData[nextPoint.y* steps_PL + nextPoint.x] == 255)//遇到边界点
continue;
pFlagData[nextPoint.y* steps_PL + nextPoint.x] = 0;
stack.push_back(nextPoint);
visited[nextPoint.y][nextPoint.x] = true;
}
}
for(i=0; iheight; i++)
delete[] visited[i];
delete[] visited;
cvReleaseImage(&tempImg);
}
//实现平滑；
void Gaussion(IplImage *img)
{
int iTempW=3,iTempH=3;
int iTempMY=1,iTempMX=1;
float fCoef = 0.065;
float fpArray[]={1,2,1,2,4,2,1,2,1}; //平滑模板；
float fResult;
BYTE * m_pData=(BYTE *)img->imageData;
int steps_PL = img->widthStep * 8 / img->depth;
int m_height=img->height;
int m_width=img->width;
BYTE **pTemp=new BYTE*[m_height];
for (int i=0;i255)
m_pData[steps_PL*(m_height-1-i)+j]=255;
else
m_pData[steps_PL*(m_height-1-i)+j]=(BYTE)(fResult+0.5);
}
}
for(i=0;imaxg1)
{
maxg1=sum1;
pos1.y=centre.y+(5*i+2);
}
if(sum2>maxg2)
{
maxg2=sum2;
pos2.x=centre.x-(5*i+2);
pos2.y=centre.y+(5*i+2);
}
if(sum3>maxg3)
{
maxg3=sum3;
pos3.x=centre.x-(5*i+2);
}
if(sum4>maxg4)
{
maxg4=sum4;
pos4.x=centre.x-(5*i+2);
pos4.y=centre.y-(5*i+2);
}
if(sum5>maxg5)
{
maxg5=sum5;
pos5.y=centre.y-(5*i+2);
}
if(sum6>maxg6)
{
maxg6=sum6;
pos6.x=centre.x+(5*i+2);
pos6.y=centre.y-(5*i+2);
}
if(sum7>maxg7)
{
maxg7=sum7;
pos7.x=centre.x+(5*i+2);
}
if(sum8>maxg8)
{
maxg8=sum8;
pos8.x=centre.x+(5*i+2);
pos8.y=centre.y+(5*i+2);
}
}
//取八个方向上的半径的平均值为最后求得的半径；
radius=0;
radius=radius+sqrt((pos1.x-centre.x)*(pos1.x-centre.x)+(pos1.y-centre.y)*(pos1.y-centre.y))
+sqrt((pos2.x-centre.x)*(pos2.x-centre.x)+(pos2.y-centre.y)*(pos2.y-centre.y))
+sqrt((pos3.x-centre.x)*(pos3.x-centre.x)+(pos3.y-centre.y)*(pos3.y-centre.y))
+sqrt((pos4.x-centre.x)*(pos4.x-centre.x)+(pos4.y-centre.y)*(pos4.y-centre.y))
+sqrt((pos5.x-centre.x)*(pos5.x-centre.x)+(pos5.y-centre.y)*(pos5.y-centre.y))
+sqrt((pos6.x-centre.x)*(pos6.x-centre.x)+(pos6.y-centre.y)*(pos6.y-centre.y))
+sqrt((pos7.x-centre.x)*(pos7.x-centre.x)+(pos7.y-centre.y)*(pos7.y-centre.y))
+sqrt((pos8.x-centre.x)*(pos8.x-centre.x)+(pos8.y-centre.y)*(pos8.y-centre.y));
radius=radius/8+2;
assert(radius<55);
return radius;
}
//输入：轮廓半径； 输出：记录轮廓的缓冲区，轮廓点序列数组；
int GetContourPoint(int radius , USHORT * Contour_buf,int * x_array,int * y_array)
{
const int centre=63;
CPoint lefttop,rightbottom;
lefttop.x=lefttop.y=63-radius;
rightbottom.x=rightbottom.y=63+radius;
int i,j;
int dist;
int index=0;
//将到中心的距离为radius的点记为轮廓点；
for(i=lefttop.x;i<=rightbottom.x;i++)
for(j=lefttop.y;j<=rightbottom.y;j++)
{
dist=sqrt((i-centre)*(i-centre)+(j-centre)*(j-centre));
if(dist==radius)
{
*(Contour_buf+i*125+j)=4095;
x_array[index]=i;
y_array[index]=j;
index++;
}
}
return index;
}
/****计算轮廓点数和肿块的面积，若轮廓点数和肿块的面积大于特定的阈值，就要重新获得初始轮廓，
若重新定初始轮廓，返回true,否则返回false ***/
bool IsGetRoundContour(IplImage * RoiRegion)
{
const int cont_threshold=200; //定义轮廓点阈值为200；
const int area_threshold=1200; //定义肿块面积阈值为1200；
int steps_PL = RoiRegion->widthStep * 8 / RoiRegion->depth;
BYTE * pData = (BYTE*) RoiRegion->imageData;
int i,j;
//求肿块的轮廓点数和肿块的面积；
int contnum,areanum;
contnum=areanum=0;
for(i=0;i<125;i++)
for(j=0;j<125;j++)
{
if(pData[i*steps_PL+j]==255)
contnum++;
else if(pData[i*steps_PL+j]==178)
areanum++;
else
;
}
if(contnum>=cont_threshold&&areanum>=area_threshold) //判断；
return true;
else
return false;
}
//输入：img
//输出：resultImg
//功能：最小二乘法背景趋势去除；
void BackgrdCorrectMatlab(IplImage *img,IplImage * resultImg)
{
cvSetZero(resultImg);
//求最小二乘法拟和的平面参数
int i,j;
double a,b,c;//参数
int steps_PL = img->widthStep * 8 / img->depth;
BYTE * pData = (BYTE*) img->imageData;
BYTE* pDest = (BYTE*) resultImg->imageData;
int nRows,nCols;
nRows = nCols = img->width;
double sumII,sumJJ,sumIJ,sumI,sumJ,sum1;
sumII=sumJJ=sumIJ=sumI=sumJ=sum1=0;
double sumIPij,sumJPij,sumPij;
sumIPij=sumJPij=sumPij=0;
for(i=0;i max)
max = rawResult[i][j];
if(rawResult[i][j] < min)
min = rawResult[i][j];
}
}
//进行标定，保存成结果图像
for(i=0;idepth == 8);
assert(PImage->nChannels == 1);
int i, j;
int width = PImage->width;
int height = PImage->height;
int bytesPL = 8 * PImage->widthStep / PImage->depth;//4字节对齐
BYTE* pData = (BYTE*)PImage->imageData;
BYTE* pMask ;
if(mask != NULL)
pMask = (BYTE*)mask->imageData;
else
pMask = NULL;
for(i=0; i 0)
histogram[value] ++;
}
}
}
//实现Gama校正
IplImage * GamaCorrectMatlab(IplImage *image)
{
int i,j;
const double threshold=0.005;
BYTE* pData = (BYTE*) image->imageData;
int steps_PL = image->widthStep * 8 / image->depth;
int height,width;
height = image->height;
width = image->width;
double maxValue,minValue;
maxValue=pData[0];
minValue=pData[0];
int count = (int)(threshold * image->width * image->height);
int* histogram = new int[256];
CalHistogram8U(image, histogram, 256, NULL);
int sumLo = 0;
int sumHi = 0;
for(i=0; i<255; i++)
{
sumLo += histogram[i];
if(sumLo > count)
{
minValue = i;
break;
}
}
for(i=255; i>=0; i--)
{
sumHi += histogram[i];
if(sumHi > count)
{
maxValue = i;
break;
}
}
//进行标定，保存成结果图像
for(i=0;i 1)
rate = 1;
if(rate < 0)
rate = 0;
double temp = 255 * rate * rate ;
pData[i*steps_PL + j] = (BYTE) temp;
}
}
delete histogram;
return image;
}
double MaxInterstify(int m,double a[])
{
double m_pix1=0;
for(int i=0;i<=m;i++)
{
m_pix1+=a[i];
}
double m_pix2=1-m_pix1;
double h=0;
for(int j=45;j<=m;j++)
{
if(a[j] > 0.000001)
h-=a[j]*log(a[j]);
}
h=h/m_pix1;
double h2=0;
for(int t=m+1;t<256;t++)
{
if(a[t] > 0.000001)
h2-=a[t]*log(a[t]);
}
h2=h2/m_pix2;
h+=h2;
h+=logf(m_pix1);
h+=logf(m_pix2);
return h;
}
//最大熵分割法
void MaxEntropy(IplImage *image)
{
int m_height,m_width;
m_height=image->height;
m_width=image->width;
BYTE * m_pData=(BYTE *)image->imageData;
int steps_PL = image->widthStep * 8 / image->depth;
double m_Intensify[256]={0},sum_pix=0;
for(int i=0;imaxen)
{
maxen=entropy[i];
max_inter=i;
}
for(i=0;imax_inter)
m_pData[i*steps_PL + j]=255;
else
m_pData[i*steps_PL + j]=0;
}
}
//输入：原始图像 img
//输出和返回：经过背景趋势去除和最大熵方法得到的肿块区域图，返回的图像是一个二值图像，
// 肿块区域灰度值为255，背景区域灰度值为0；
IplImage * BackgroundTrendRemoveAndMaxEntropy(IplImage * img)
{
IplImage * resultimg;
resultimg=cvCreateImage(cvSize(125, 125), 8, 1);
BackgrdCorrectMatlab(img,resultimg);
GamaCorrectMatlab(resultimg); //Gama校正；
Gaussion(resultimg); //高斯滤波；
MaxEntropy(resultimg); //最大熵方法分割肿块；
return resultimg; //返回结果图像；
}
//输入：img
//输出：resultimg
//功能：实现背景趋势去除；
void BackgroundTrendRemove(IplImage * img,IplImage *resultimg)
{
BackgrdCorrectMatlab(img,resultimg);
GamaCorrectMatlab(resultimg); //Gama校正；
}
/*******基于四邻域取肿块区域轮廓，令轮廓点灰度值为255，背景区域灰度为0，
肿块内部点灰度为178 ******/
void GetRegionContour(IplImage * ImgRegion) //ImgRegion为8位区域图，
//肿块轮廓和背景的灰度为0，肿块区域内部灰度为178。
{
int stepsPL= ImgRegion->widthStep;
BYTE * pData=(BYTE*)ImgRegion->imageData;
int i,j;
BYTE top,bottom,left,right;
for(i=0;iheight;i++)
for(j=0;jwidth;j++)
{
if((pData[stepsPL*i+j]==0)&&i>1&&j>1&&i<(ImgRegion->height-1)&&j<(ImgRegion->width-1))
{ //四领域；
top=pData[stepsPL*(i-1)+j];
bottom=pData[stepsPL*(i+1)+j];
left=pData[stepsPL*i+(j-1)];
right=pData[stepsPL*i+(j+1)];
if(top==178||bottom==178||left==178||right==178)
{
pData[stepsPL*i+j]=255;
}
}
}
}
//找轮廓，并将轮廓灰度值值为255；
void FindContour(IplImage * roi_region)
{
int i,j;
int stepsPL =roi_region->widthStep;
BYTE * pData = (BYTE*)roi_region->imageData;
int lheight=roi_region->height;
int lwidth=roi_region->width;
CvMemStorage* storage = cvCreateMemStorage(0);
CvSeq* contours = NULL;
for(i=0;itotal);
cvCvtSeqToArray(contours, pointArray, CV_WHOLE_SEQ);
int total = contours->total;
for(i=0; iwidthStep;
BYTE * pData = (BYTE*)roi_region->imageData;
visited[seed.x][seed.y]=true; //入队时置1
CList queue(10) ;
//种子入队
queue.AddTail(seed);
long countInRegion=1;
while(!queue.IsEmpty()) //队列非空
{
CPoint curPoint;
//取队列头元素
curPoint=(CPoint)queue.RemoveHead();
//扩展种子点的8领域
for(i=-1;i<=1;i++)
{
for(j=-1;j<=1;j++)
{
if(curPoint.x+i<0 || curPoint.x+i>=125 ||
curPoint.y+j<0 || curPoint.y+j>=125 ) //越界处理
continue;
if(i==0 && j==0 )
continue;
lpSrc=pData[(curPoint.x+i)*stepsPL+(curPoint.y+j)];
if(!visited[curPoint.x+i][curPoint.y+j])
{
visited[curPoint.x+i][curPoint.y+j]=true;
if(lpSrc>0)
{
//入队
CPoint newSeed(curPoint.x+i,curPoint.y+j);
queue.AddTail(newSeed);
countInRegion++;
}
}
}
}
} //队列为空
return countInRegion;
}
//找最大的连通区域作为肿块区域
void FindMostLargeRegion(IplImage * roi_region)
{
CPoint resultCenter(0,0); //最大连通区域内部的一点；
LONG maxCountInRegion=0; //最大连通区域的点数；
LONG tempCountInRegion; //临时变量；
int stepsPL =roi_region->widthStep;
BYTE * pData = (BYTE*)roi_region->imageData;
int lHeight=roi_region->height;
int lWidth=roi_region->width;
int i,j;
USHORT lpSrc;
bool** visited=new bool*[lHeight]; //访问标志
for(i=0;i0&&(visited[i][j]==false)) //Not the background and not visited
{
CPoint startPoint(i,j);
tempCountInRegion=RegionGrowth(startPoint,visited,roi_region); //区域增长；
if (tempCountInRegion>maxCountInRegion)
{
maxCountInRegion=tempCountInRegion;
resultCenter.x=i;
resultCenter.y=j;
} //if
} //if
else
visited[i][j]=true;
} //for
for(i=0;iwidthStep / enhanceImg->depth;
BYTE * pRoiData = (BYTE*)enhanceImg->imageData;
int stepsPLRegion= roi.m_region->widthStep;
BYTE * pRegionData=(BYTE*)roi.m_region->imageData;
for(i=0;i<125;i++)
for(j=0;j<125;j++)
{
*(img_buf+125*i+j)=(USHORT)(pRoiData[i*stepsPLRoi+j]*16); //将enhanceImg拷贝到img_buf中；
*(Contour_buf+125*i+j)=0;
}
if(IsGetRoundContour(roi.m_region)) //判断是否重新取初始轮廓；
{
USHORT * Para_a, * Para_b, * Para_c, * OrigImg,*EdgeMap_buf;
int rows=125;
int cols=125;
Para_a = svector(rows, cols);
Para_b = svector(rows, cols);
Para_c = svector(rows, cols);
OrigImg = svector(rows, cols);
EdgeMap_buf= svector(rows, cols);
int radius;
copy_svector(img_buf, OrigImg, rows, cols);
MinQuare(Para_a,Para_b,Para_c,EdgeMap_buf,OrigImg); //最小二乘法就梯度图；
radius=FindMassRadius(EdgeMap_buf); //求新的圆形初始轮廓的半径；
GetContourPoint(radius,Contour_buf,x_array,y_array); //获得新的轮廓点数；
//释放空间；
free_svector(Para_a);
free_svector(Para_b);
free_svector(Para_c);
free_svector(OrigImg);
free_svector(EdgeMap_buf);
}
else
{
int index;
index=0;
for(i=0;i<125;i++)
for(j=0;j<125;j++)
{
if(pRegionData[i*stepsPLRegion+j]==255) //轮廓点的灰度值为最大值255；
{
*(Contour_buf+125*i+j)=4095;
x_array[index]=i; //初始化轮廓点数组；
y_array[index]=j;
index++;
}
else
{
*(Contour_buf+125*i+j)=0; //非轮廓点的灰度值令其为0；
}
}
}
Contour_count=record_contour_array(Contour_buf,125,125,x_array,y_array); //将曲线上的点按行走路线有序存储；
ActiveContourRegion(growth_buf,125,125,Contour_count,x_array,y_array); //求包含肿块区域的矩形区域，并记录在缓冲区growth_buf中；
Contour_count=Computing_ActiveContour(img_buf,growth_buf,125,125,Contour_count,x_array,y_array); //计算活动轮廓；
//令轮廓灰度值为255，其他点的灰度值为178；
for(i=0;i<125;i++)
for(j=0;j<125;j++)
{
pRegionData[i*stepsPLRegion+j]=178;
}
for(k=0;kwidthStep / enhanceImg->depth;
BYTE * pRoiData = (BYTE*)enhanceImg->imageData;
int stepsPLRegion= roi.m_region->widthStep;
BYTE * pRegionData=(BYTE*)roi.m_region->imageData;
int index;
index=0;
for(i=0;i<125;i++)
for(j=0;j<125;j++)
{
*(img_buf+125*i+j)=(USHORT)(pRoiData[i*stepsPLRoi+j]*16); //将enhanceImg拷贝到img_buf中；
if(pRegionData[i*stepsPLRegion+j]==255) //轮廓点的灰度值为最大值255；
{
*(Contour_buf+125*i+j)=4095;
x_array[index]=i;
y_array[index]=j; //初始化轮廓点数组；
index++;
}
else
{
*(Contour_buf+125*i+j)=0; //非轮廓点的灰度值令其为0；
}
}
Contour_count=record_contour_array(Contour_buf,125,125,x_array,y_array); //将曲线上的点按行走路线有序存储；
ActiveContourRegion(growth_buf,125,125,Contour_count,x_array,y_array); //求包含肿块区域的矩形区域，并记录在缓冲区growth_buf中；
Contour_count=Computing_ActiveContour(img_buf,growth_buf,125,125,Contour_count,x_array,y_array); //计算活动轮廓；
//令轮廓灰度值为255，其他点的灰度值为178；
for(i=0;i<125;i++)
for(j=0;j<125;j++)
{
pRegionData[i*stepsPLRegion+j]=178;
}
for(k=0;k