// SnakeProcess.cpp: implementation of the SnakeProcess class.
//
//////////////////////////////////////////////////////////////////////
#include "stdafx.h"
#include "mammocad.h"
#include "SnakeProcess.h"
#include "ROI.h"
#include "imageprocess.h"
#include <fstream.h>
#include <list>
#include "Matrix.h"
#include <afxtempl.h>
#include <VECTOR>
#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) (a<b? a: b)
#define MAX_12_BIT 4095
typedef unsigned char PGbyte;
typedef char PGchar;
typedef unsigned short PGushort;
typedef short PGshort;
typedef unsigned int PGuint;
typedef int PGint;
typedef float PGfloat;
typedef double PGdouble;
typedef void PGvoid;
typedef unsigned short USHORT;
/* Prints error message: modified to work in mex environment */
void MinQuare(USHORT *Para_a, USHORT *Para_b, USHORT *Para_c,USHORT *Para_d,USHORT *pDataSrc)
{
int i,j,k,l;
int nRows,nCols;
nRows = 125;
nCols = 125;
double* pData_a = new double[nRows*nCols];
double* pData_b = new double[nRows*nCols];
double* pData_c = new double[nRows*nCols];
double* pData_d = new double[nRows*nCols];
int sizeWnd=5;
for(k=sizeWnd/2;k<nRows-sizeWnd/2;k++)
{
for(l=sizeWnd/2;l<nCols-sizeWnd/2;l++)
{
double a,b,c;
double sumII,sumJJ,sumIJ,sumI,sumJ,sum1;
sumII=sumJJ=sumIJ=sumI=sumJ=sum1=0;
double sumIPij,sumJPij,sumPij;
sumIPij=sumJPij=sumPij=0;
for(i=k-sizeWnd/2;i<k+sizeWnd/2+1;i++)
{
for(j=l-sizeWnd/2;j<l+sizeWnd/2+1;j++)
{
sumII += (i-k)*(i-k);
sumJJ += (j-l)*(j-l);
sumIJ += (i-k)*(j-l);
sumI += (i-k);
sumJ += (j-l);
sum1 += 1;
sumIPij += (double)(i-k) * (*(pDataSrc+125*i+j));
sumJPij += (double)(j-l) *(*(pDataSrc+125*i+j));
sumPij += (*(pDataSrc+125*i+j));
}
}
double value[9]={
sumII,sumIJ,sumI,
sumIJ,sumJJ,sumJ,
sumI ,sumJ ,sum1
};
CMatrix matrix(3,3,value);
BOOL flag = matrix.InvertSsgj();
if(flag)
{
double valueTmp[3]={
sumIPij,sumJPij,sumPij
};
CMatrix matrixTmp(3,1,valueTmp);
CMatrix para(3,1);
para = matrix*matrixTmp;
a=para.GetElement(0,0);
b=para.GetElement(1,0);
c=para.GetElement(2,0);
}
pData_a[k*nCols+l]=a;
pData_b[k*nCols+l]=b;
pData_c[k*nCols+l]=c;
double sin,cos,r;
r=sqrt((k-nRows/2.0)*(k-nRows/2.0)+(l-nCols/2.0)*(l-nCols/2.0));
sin = (nRows/2.0-k)/r;
cos = (l-nCols/2.0)/r;
pData_d[k*nCols+l]=a*sin-b*cos;
}//End for(l=2;l<nCols-2;l++)
}//End for(k=2;k<nRows-2;k++)
//进行标定
float min_a,min_b,min_c,min_d,max_a,max_b,max_c,max_d;
min_a=max_a=pData_a[sizeWnd/2*nCols+sizeWnd/2];
min_b=max_b=pData_b[sizeWnd/2*nCols+sizeWnd/2];
min_c=max_c=pData_c[sizeWnd/2*nCols+sizeWnd/2];
min_d=max_d=pData_d[sizeWnd/2*nCols+sizeWnd/2];
for(k=sizeWnd/2;k<nRows-sizeWnd/2;k++)
{
for(l=sizeWnd/2;l<nCols-sizeWnd/2;l++)
{
//Para_a
if(pData_a[k*nCols+l]<min_a)
{
min_a=pData_a[k*nCols+l];
}
else if(pData_a[k*nCols+l]>max_a)
{
max_a=pData_a[k*nCols+l];
}
//Para_b
if(pData_b[k*nCols+l]<min_b)
{
min_b=pData_b[k*nCols+l];
}
else if(pData_b[k*nCols+l]>max_b)
{
max_b=pData_b[k*nCols+l];
}
//Para_c
if(pData_c[k*nCols+l]<min_c)
{
min_c=pData_c[k*nCols+l];
}
else if(pData_c[k*nCols+l]>max_c)
{
max_c=pData_c[k*nCols+l];
}
//Para_d
if(pData_d[k*nCols+l]<min_d)
{
min_d=pData_d[k*nCols+l];
}
else if(pData_d[k*nCols+l]>max_d)
{
max_d=pData_d[k*nCols+l];
}
}//End for(l=2;l<nCols-2;l++)
}//End for(k=2;k<nRows-2;k++)
for(k=sizeWnd/2;k<nRows-sizeWnd/2;k++)
{
for(l=sizeWnd/2;l<nCols-sizeWnd/2;l++)
{
//Para_a
*(Para_a+125*k+l) = 4095*(pData_a[k*nCols+l]-min_a)/(max_a-min_a);
//Para_b
*(Para_b+125*k+l) = 4095*(pData_b[k*nCols+l]-min_b)/(max_b-min_b);
//Para_c
*(Para_c+125*k+l) = 4095*(pData_c[k*nCols+l]-min_c)/(max_c-min_c);
//Para_d
*(Para_d+125*k+l) = 4095*(pData_d[k*nCols+l]-min_d)/(max_d-min_d);
}//End for(l=2;l<nCols-2;l++)
}//End for(k=2;k<nRows-2;k++)
}
PGvoid pgError(char error_text[])//错误消息处理
{
char buf[255];
sprintf(buf, "GVFC.MEX: %s", error_text);
//mexErrMsgTxt(buf);
return;
}
/* Allocates a matrix of doubles with range [nrl..nrh][ncl..nch] */
PGdouble **pgDmatrix(int nrl, int nrh, int ncl, int nch)//分配二维数组,数组下标可以任意开始A[3..10][5..15]
{
int j;
long bufsize,bufptr;
PGdouble **m;
nrh++;
nch++;
bufsize = (nrh-nrl+1)*sizeof(PGdouble*)
+ (nrh-nrl+1)*(nch-ncl+1)*sizeof(PGdouble);//节省空间,sizeof(PGdouble*)=4
m=(PGdouble **) malloc(bufsize);
if (!m) pgError("allocation failure 1 in pgDmatrix()");
m -= nrl;//如果此处不这样写的话,下面m[j]就不是那个位置了!
bufptr = ((long) (m)) + (nrh-nrl+1)*sizeof(PGdouble*);//相当于到数据区的偏移量
for(j=nrl;j<=nrh;j++) {
m[j] = ((PGdouble *) (bufptr+(j-nrl)*(nch-ncl+1)*sizeof(PGdouble)));//此时相当于m+j
m[j] -= ncl;//此时将其地址向前移动,目的是调用时m[nrl][ncl]相当于m[0][0]
}
return m;
}
/* Frees a matrix allocated by dmatrix */
PGvoid pgFreeDmatrix(PGdouble **m, int nrl, int nrh, int ncl, int nch)//二维数组指针释放
{
free((char*) (m+nrl));
return;
}
double *dvector(long nl, long nh)
{
double *v;
v=(double *)malloc((size_t)((nl*nh)*sizeof(double)));
if (!v)
{
AfxMessageBox("内存分配失败!");
//exit(1);
}
return (v);
}
int *ivector(long nl, long nh)//申请int类型向量空间
{
int *v;
v = (int *)malloc((size_t) ((nh-nl+1+1)*sizeof(int)));
if( !v )
exit(1);
return v-nl+1;//?
}
USHORT *svector(long nl, long nh)//申请USHORT类型向量空间//?
{
USHORT *v;
v = (USHORT *)malloc((size_t) ((nl*nh)*sizeof(USHORT)));
if( !v )
exit(1);
return v;
}
void free_dvector(double *v, long nl, long nh)//释放double类型向量空间
{
free((char*) (v));
}
void free_ivector(int *v, long nl, long nh)//释放int类型向量空间
{
free((char*) (v+nl-1));
}
void free_svector(USHORT *v)//释放double类型向量空间//?
{
free(v);
}
void copy_svector(USHORT* img_buf, USHORT* EdgeMap_buf, int rows,int cols)
{
memcpy(EdgeMap_buf, img_buf, rows*cols*sizeof(USHORT));
}
SnakeMoxing.rar_snake
版权申诉
143 浏览量
2022-09-19
21:28:56
上传
评论
收藏 19KB RAR 举报
我虽横行却不霸道
- 粉丝: 75
- 资源: 1万+
最新资源
- TG-2024-05-23-204718255.mp4
- 候志强@181 5428 8938_20240420112107.amr
- spispispispispi
- 实验二:IP协议分析.zip
- 驱动代码驱动代码驱动代码驱动代码
- SVID_20240523_141155_1.mp4
- Code for the complete guide to tkinter tutorial
- 关于百货中心供应链管理系统.zip
- SimpleFolderIcon-master 修改Unity的Project下的文件夹图标
- A python Tkinter widget to display tile based maps
资源上传下载、课程学习等过程中有任何疑问或建议,欢迎提出宝贵意见哦~我们会及时处理!
点击此处反馈