mc算法实现

1. 什么是marching cubes算法具体怎么讲的

Marching Cubes算法（医学图像三维绘制中的面绘制）2007-08-16 00:50建议要看的资料
[1] Lorensen W E, Cline H E .Marching cubes: a high-resoulution 3D suface construction algorithm [J], Computer Graphics,1987, 21(4):163~169
[2]集成化医学影像算法平台理论与实践田捷，赵明昌，何晖光 清华大学出版社2004年10月
[3]Polygonising a scalar field Also known as: "3D Contouring", "Marching Cubes", "Surface Reconstruction"
http://local.wasp.uwa.e.au/~pbourke/geometry/polygonise/Marching Cubes；
[4]www.3dmed.net

Marching Cubes算法工作原理
Marching Cubes算法是三维数据场等值面生成的经典算法，是体素单元内等值面抽取技术的代表。
等值面是空间中所有具有某个相同值的点的集合。它可以表示为， ，c是常数。则称F(f)为体数据f中的等值面。
在MC算法中，假定原始数据是离散的三维空间规则数据场。用于医疗诊断的断层扫描(CT)及核磁共振成像(MRI) 等产生的图像均属于这一类型。MC算法的基本思想是逐个处理数据场中的体素，分类出与等值面相交的体素，采用插值计算出等值面与体素棱边的交点。根据体素中每一顶点与等值面的相对位置，将等值面与立方体边的交点按一定方式连接生成等值面，作为等值面在该立方体内的一个逼近表示。在计算出关于体数据场内等值面的有关参数后山常用的图形软件包或硬件提供的面绘制功能绘制出等值面。

图2.1 离散的三维空间规则数据场中的一个体素
2.1.1 MC算法的主要步骤
1． 确定包含等值面的体素
离散的三维空间规则数据场中的一个体素可以用图2.1表示。8个数据点分别位于该体素的8个角点上。MC算法的基本假设是:沿着体素的边其数据场呈局部连续线性变化，根据这个假设，可认为，如果两个相邻采样点一个为正点，一个为负点，则它们连成的边上一定存在且仅有一个等值点 (设等值面值为c)。如果得到了体素各条边上的等值点，就可以以这些点为顶点，用一系列的三角形拟合出该体素中的等值面。因此确定立方体体素中等值面的分布是该算法的基础。
首先对体素的8个顶点进行分类，以判断其顶点是位于等值面之外，还是位于等值面之内。再根据8个顶点的状态，确定等值面的剖分模式。顶点分类规则为:
1. 如体素顶点的数据值大于或等于等值面的值，则定义该顶点位于等值面之外， 记为正点，即“1“
2. 如体素顶点的数据值小于等值面的值，则定义该顶点位于等值面之内，记为负点， 即“0"
由于每个体素共有8个顶点，且每个顶点有正负两种状态，所以等值面可能以 =256种方式与一个体素相交。通过列举出这256种情况，就能创建一张表格，利用它可以查出任意体素中的等值面的三角面片表示。如果考虑互补对称性，将等值面的值和8个角点的函数值的大小关系颠倒过来，即将体素的顶点标记置反(0变为1, 1变为0)，这样做不会影响该体素的8个角点和等值面之间的拓扑结构，可将256种方式简化成128种。其次，再利用旋转对称性，可将这128种构型进一步简化成15种。图3.2给出了这15种基本构型[131其中黑点标记为“1”的角点。

图2.2 分布状态表

图2.3 体素角点分布不同情况
基于上面的分析，MC算法中用一个字节的空间构造了一个体素状态表，如图2.2所示，该状态表中的每一位可表示出该体元中的一个角点的0或1的状态。根据这一状态表，就可知道当前体素属于图2.3中哪一种情况，以及等值面将与哪一条边相交。
2.求等值面与体元边界的交点
在确定体素内三角剖分模式后，就要计算三角片顶点位置。当三维离散数据场的密度较高时，即当体素很小时，可以假定函数值沿体素边界呈线性变化，这就是MC算法的基本假设。因此，根据这一基本假设，可以直接用线性插值计算等值面与体素边的交点。
对于当前被处理体素的某一条边，如果其两顶点 ， 的标记值不同，那么等值面一定与此边相交，且仅有一个交点。交点为 其中P代表等值点坐标, , 代表两个端点的坐标， ， 代表两个端点的灰度值，v代表域值。求出等值面与体素棱边的交点以后，就可以将这些交点连接成三角形或多边形，形成等值面的一部分。
3.等值面的法向量计算
为了利用图形硬件显示等值面图象，必须给出形成等值面的各三角面片的法向分量，选择适当的局部面光照模型进行光照计算，以生成真实感图形。
对于等值面上的每一点，其沿面的切线方向的梯度分量应该是零，因此，该点的梯度矢量的方向也就代表了等值面在该点的法向量，当梯度值非零。所幸的是等值面往往是由两种具有不同密度的物质的分解面，因此其上的每点的梯度矢量均不为零，即
Mc算法采用中心差分方法求采样点p〔m ,n, k ) 处的梯度矢量，公式如下:

Gx=〔g(i+1,j,k)-g(i-1,j,k)〕/2dx
Gy=〔g(i,j+1,k)-g(i,j-1,k)〕/2dy
Gz=〔g(i,j,k+1)-g(i,j,k-1)〕/2dz
其中D（i，j ,k）是切片k在像素（i，j）的密度， , ， 是立方体边的长度。对g进行归一化，得到(gx/｜g｜，gy/｜g｜，gz/｜g｜）作为（i,j,k)上的单位法向量。然后，对体素八个顶点上法向量进行线性插值就可得到位于体素棱边上的三角片的各个顶点上的法向量。设计算得到的某个三角片的三个顶点上的单位法向量分别为( , 和 )，这个三角片的几何重心为 ，则该三角片的法向量起始于 ，终止于 。代入Gourand光照模型公式，就可计算出小三角片表面的光强(灰度)。将其投影在某个特定的二维平面上进行显示，从而显示出物体富有光感的整个表面形态。其中我们在内存中保留四个切片来计算立方体中所有顶点梯度。
2.1.2 MC算法流程
1、将三维离散规则数据场分层读入内存;
2、扫描两层数据，逐个构造体素，每个体素中的8个角点取自相邻的两层;
3、将体素每个角点的函数值与给定的等值面值c做比较，根据比较结果，构造
该体素的状态表;
4、根据状态表，得出将与等值面有交点的边界体素;
5、通过线性插值方法计算出体素棱边与等值面的交点;
6、利用中心差分方法，求出体素各角点处的法向量，再通过线性插值方法，求出三角面片各顶点处的法向;
7,根据各三角面片上各顶点的坐标及法向量绘制等值面图像。
========================================================
MC代码
MarchingCubes(float lowThreshold,float highThreshold,float XSpace,float YSpace,float ZSpace)
{
//记录生成的顶点数和面数初始时应该为0
m_vNumber=0;
m_fNumber=0;
//当前Cube中生成的顶点和面数
int vertPos,facePos;
//包围盒的尺寸 用于绘制程序计算整个场景的包围盒，用于调整观察位置，以使整个场景尽可能占满整个窗口。
float min[3],max[3];
min[0]=min[1]=min[2]=max[0]=max[1]=max[2]=0;//初始化
//当前扫描层的切片数据和一个临时的切片数据
short *pSliceA,*pSliceB,*pSliceC,*pSliceD,*tempSlice;

pSliceA=pSliceB=pSliceC=tempSlice=NULL;
int imageWidth,imageHeight,imageSize,sliceNumber;
imageWidth=imageHeight=512;//我们是512×512的数据
imageSize=imageWidth*imageHeight;
sliceNumber=m_FileCount-1;
if((highThreshold*lowThreshold)==0)
{
return 0;
}
pSliceD =new short [imageSize];
//因为等值面是每相邻两层切片为单位进行提取的，所以在处理后两层切片时难免生成前两层切片已经生成的顶点，这时候就用下面的数组记录哪些边上的顶点已经生成了，如果遇到已经生成的顶点就不再重复计算而是直接使用记录的索引，否则就生成新的顶点。
long *bottomXEdge=new long[imageSize];
long *bottomYEdge=new long[imageSize];
long *topXEdge=new long[imageSize];
long *topYEdge=new long[imageSize];
long *zEdge=new long[imageSize];

tempSlice=new short [imageSize];

if(bottomXEdge==NULL||bottomYEdge==NULL||
topXEdge==NULL||topYEdge==NULL||
zEdge==NULL||tempSlice==NULL)
{
return 0;//错误
}
//初始化数据
memset(bottomXEdge,-1,sizeof(long)*imageSize);
memset(bottomYEdge,-1,sizeof(long)*imageSize);
memset(topXEdge,-1,sizeof(long)*imageSize);
memset(topYEdge,-1,sizeof(long)*imageSize);
memset(zEdge,-1,sizeof(long)*imageSize);
memset(tempSlice,0,sizeof(short)*imageSize);

//计算某一层顶点和三角时所需要的一些变量
//一些循环变量
int i,j,k,w,r;
//cube类型
unsigned char cubeType(0);
//计算法向量
float dx[8],dy[8],dz[8],squaroot;
//记录某个Cube生成
float vertPoint[12][6];
int cellVerts[12]; //what use
int triIndex[5][3]; //每个cube最多生成5条边

//用于记录已生成顶点索引的临时变量
int offset;
//当前cube8个顶点的灰度值
short cubegrid[8];
long *edgeGroup;
//得到数据

pSliceD=m_volumeData;
pSliceB=tempSlice;
pSliceA=tempSlice;

int tt,tt1;

//扫描4层切片的顺序
/*
-----------------------D |
-----------------------B |
-----------------------C |
-----------------------A |
V
*/
//marching cubes 算法开始实行 ？第一次循环时，只读入一个切片？
for(i=0;i<=(sliceNumber);i++)
{
pSliceC=pSliceA;
pSliceA=pSliceB;
pSliceB=pSliceD;
if(i>=sliceNumber-2)
{
pSliceD=tempSlice;
}
else
{

pSliceD+=imageSize;
}
for(j=0;j<imageHeight-1;++j)
for(k=0;k<imageWidth-1;++k)
/* for(j=10;j<imageHeight-5;j++)//调试用
for(k=10;k<imageWidth-5;k++)*/
{
//得到八个顶点的灰度值step0
cubegrid[0]=pSliceA[j*imageWidth+k];
cubegrid[1]=pSliceA[j*imageWidth+k+1];
cubegrid[2]=pSliceA[(j+1)*imageWidth+k+1];
cubegrid[3]=pSliceA[(j+1)*imageWidth+k];
cubegrid[4]=pSliceB[j*imageWidth+k];
cubegrid[5]=pSliceB[j*imageWidth+k+1];
cubegrid[6]=pSliceB[(j+1)*imageWidth+k+1];
cubegrid[7]=pSliceB[(j+1)*imageWidth+k];
//计算cube的类型
cubeType=0;
for(w=0;w<8;w++)
{

if((cubegrid[w]>lowThreshold)&&(cubegrid[w]<highThreshold))//需要画的点

{
cubeType|=(1<<w);
}
}//end for计算cube的类型
if(cubeType==0||cubeType==255)
{
continue;
}
for(w=0;w<12;w++) //初始化cellVerts表到零
{
cellVerts[w]=-1;
}
//计算6个方向相邻点的象素差值（用于计算法向量）
if(k==0)
{
dx[0]=pSliceA[j*imageWidth+1];
dx[3]=pSliceA[(j+1)*imageWidth+1];
dx[4]=pSliceB[j*imageWidth+1];
dx[7]=pSliceB[(j+1)*imageWidth+1];
}
else
{
dx[0]=pSliceA[j*imageWidth+k+1]
-pSliceA[j*imageWidth+k-1];
dx[3]=pSliceA[(j+1)*imageWidth+k+1]
-pSliceA[(j+1)*imageWidth+k-1];
dx[4]=pSliceB[j*imageWidth+k+1]
-pSliceB[j*imageWidth+k-1];
dx[7]=pSliceB[(j+1)*imageWidth+k+1]
-pSliceB[(j+1)*imageWidth+k-1];
}
if(k==imageWidth-2)
{
dx[1]=-pSliceA[j*imageWidth+imageWidth-2];
dx[2]=-pSliceA[(j+1)*imageWidth+imageWidth-2];
dx[5]=-pSliceB[j*imageWidth+imageWidth-2];
dx[6]=-pSliceB[(j+1)*imageWidth+imageWidth-2];
}
else
{
dx[1]=pSliceA[j*imageWidth+k+2]
-pSliceA[j*imageWidth+k];
dx[2]=pSliceA[(j+1)*imageWidth+k+2]
-pSliceA[(j+1)*imageWidth+k];
dx[5]=pSliceB[j*imageWidth+k+2]
-pSliceB[j*imageWidth+k];
dx[6]=pSliceB[(j+1)*imageWidth+k+2]
-pSliceB[(j+1)*imageWidth+k];
}
if(j==0)
{
dy[0]=pSliceA[imageWidth+k];
dy[1]=pSliceA[imageWidth+k+1];
dy[4]=pSliceB[imageWidth+k];
dy[5]=pSliceB[imageWidth+k+1];
}
else
{
dy[0]=pSliceA[(j+1)*imageWidth+k]
-pSliceA[(j-1)*imageWidth+k];
dy[1]=pSliceA[(j+1)*imageWidth+k+1]
-pSliceA[(j-1)*imageWidth+k+1];
dy[4]=pSliceB[(j+1)*imageWidth+k]
-pSliceB[(j-1)*imageWidth+k];
dy[5]=pSliceB[(j+1)*imageWidth+k+1]
-pSliceB[(j-1)*imageWidth+k+1];
}
if(j==imageHeight-2)
{
dy[2]=-pSliceA[(imageHeight-2)*imageWidth+k+1];
dy[3]=-pSliceA[(imageHeight-2)*imageWidth+k];
dy[6]=-pSliceB[(imageHeight-2)*imageWidth+k+1];
dy[7]=-pSliceB[(imageHeight-2)*imageWidth+k];
}
else
{
dy[2]=pSliceA[(j+2)*imageWidth+k+1]-pSliceA[j*imageWidth+k+1];
dy[3]=pSliceA[(j+2)*imageWidth+k]-pSliceA[j*imageWidth+k];
dy[6]=pSliceB[(j+2)*imageWidth+k+1]-pSliceB[j*imageWidth+k+1];
dy[7]=pSliceB[(j+2)*imageWidth+k]-pSliceB[j*imageWidth+k];
}
dz[0]=pSliceB[j*imageWidth+k]
-pSliceC[j*imageWidth+k];
dz[1]=pSliceB[j*imageWidth+k+1]
-pSliceC[j*imageWidth+k+1];
dz[2]=pSliceB[(j+1)*imageWidth+k+1]
-pSliceC[(j+1)*imageWidth+k+1];
dz[3]=pSliceB[(j+1)*imageWidth+k]
-pSliceC[(j+1)*imageWidth+k];
dz[4]=pSliceD[j*imageWidth+k]
-pSliceA[j*imageWidth+k];
dz[5]=pSliceD[j*imageWidth+k+1]
-pSliceA[j*imageWidth+k+1];
dz[6]=pSliceD[(j+1)*imageWidth+k+1]
-pSliceA[(j+1)*imageWidth+k+1];
dz[7]=pSliceD[(j+1)*imageWidth+k]
-pSliceA[(j+1)*imageWidth+k];

//计算三角形顶点的坐标和梯度
vertPos=0;
facePos=0;
for(w=0;w<12;w++)
{
if(g_EdgeTable[cubeType]&(1<<w)) //what …..
{
//根据g_edgeTable[256]对应值判断cube的那一条边与等值面有交点
switch(w)
{
case 0:
offset=j*imageWidth+k;
edgeGroup=bottomXEdge;
break;
case 1:
offset=j*imageWidth+k+1;
edgeGroup=bottomYEdge;
break;
case 2:
offset=(j+1)*imageWidth+k;
edgeGroup=bottomXEdge;
break;
case 3:
offset=j*imageWidth+k;
edgeGroup=bottomYEdge;
break;
case 4:
offset=j*imageWidth+k;
edgeGroup=topXEdge;
break;
case 5:
offset=j*imageWidth+k+1;
edgeGroup=topYEdge;
break;
case 6:
offset=(j+1)*imageWidth+k;
edgeGroup=topXEdge;
break;
case 7:
offset=j*imageWidth+k;
edgeGroup=topYEdge;
break;
case 8:
offset=j*imageWidth+k;
edgeGroup=zEdge;
break;
case 9:
offset=j*imageWidth+k+1;
edgeGroup=zEdge;
break;
case 10:
offset=(j+1)*imageWidth+k+1;
edgeGroup=zEdge;
break;
case 11:
offset=(j+1)*imageWidth+k;
edgeGroup=zEdge;
break;

}//对应switch的{。。。end for//根据g_EdgeTable对应值判断cube的那一条边与等值面有交点
//该边上的顶点是否已经在上一层中生成
if(edgeGroup[offset]==-1)
{
int index1,index2;
short s1,s2,s;
float x1,y1,z1,nx1,ny1,nz1;
float x2,y2,z2,nx2,ny2,nz2;
//得到该边两端点的索引进而得到两点的灰度值
index1=g_CoordTable[w][3];
index2=g_CoordTable[w][4];
s1=cubegrid[index1];
s2=cubegrid[index2];
if(s1<highThreshold&&s1>lowThreshold)
{
if(s2>=highThreshold)
{
s=highThreshold;
}
else if(s2<=lowThreshold)
{
s=lowThreshold;
}
}
else if(s2<highThreshold&&s2>lowThreshold)
{
if(s1>=highThreshold)
{
s=highThreshold;
}
else if(s1<=lowThreshold)
{
s=lowThreshold;
}
}
//计算两端点实际坐标
x1=(k+g_CoordVertex[index1][0])*XSpace;
y1=(j+g_CoordVertex[index1][1])*YSpace;
z1=(i+g_CoordVertex[index1][2])*ZSpace;
x2=(k+g_CoordVertex[index2][0])*XSpace;
y2=(j+g_CoordVertex[index2][1])*YSpace;
z2=(i+g_CoordVertex[index2][2])*ZSpace;
//计算两端点的法向量

nx1=dx[index1]/XSpace;
ny1=dy[index1]/YSpace;
nz1=dz[index1]/ZSpace;
nx2=dx[index2]/XSpace;
ny2=dy[index2]/YSpace;
nz2=dz[index2]/ZSpace;
float factor=((float)(s-s1))/((float)(s2-s1));

//插值计算交点坐标
vertPoint[vertPos][0]=factor*(x2-x1)+x1;
vertPoint[vertPos][1]=factor*(y2-y1)+y1;
vertPoint[vertPos][2]=factor*(z2-z1)+z1;
//计算法向量
vertPoint[vertPos][3]=factor*(nx1-nx2)-nx1;
vertPoint[vertPos][4]=factor*(ny1-ny2)-ny1;
vertPoint[vertPos][5]=factor*(nz1-nz2)-nz1;
//法向量归一化
squaroot=sqrt(vertPoint[vertPos][3]*vertPoint[vertPos][3]+vertPoint[vertPos][4]*vertPoint[vertPos][4]
+vertPoint[vertPos][5]*vertPoint[vertPos][5]);
if(squaroot<=0)squaroot=1.0;
vertPoint[vertPos][3]/=squaroot;
vertPoint[vertPos][4]/=squaroot;
vertPoint[vertPos][5]/=squaroot;
//更新包围盒数据
if(min[0]>vertPoint[vertPos][0])
{
min[0]=vertPoint[vertPos][0];
}
if(min[1]>vertPoint[vertPos][1])
{
min[1]=vertPoint[vertPos][1];
}
if(min[2]>vertPoint[vertPos][2])
{
min[2]=vertPoint[vertPos][2];
}
if(max[0]<vertPoint[vertPos][0])
{
max[0]=vertPoint[vertPos][0];
}
if(max[1]<vertPoint[vertPos][1])
{
max[1]=vertPoint[vertPos][1];
}
if(max[2]<vertPoint[vertPos][2])
{
max[2]=vertPoint[vertPos][2];
}
//记录新生成的顶点索引
cellVerts[w]=m_vNumber;
edgeGroup[offset]=cellVerts[w];
m_vNumber++;
vertPos++;
} //end if(edgeGroup[offset]==-1) ////
else
{
//若该点已经在上一层生成，则直接得到其索引
cellVerts[w]=edgeGroup[offset];
}
} // end对应if(g_EdgeTable[cubeType]&(1<<w)) //

} //对应for(w=0;w<12;w++)
//保存当前cubes 顶点和法向量
tt1=m_vNumber-vertPos;
for(tt=0;tt<vertPos;tt++)
{
vPointNomal[tt1+tt][0]=vertPoint[tt][0];
vPointNomal[tt1+tt][1]=vertPoint[tt][1];
vPointNomal[tt1+tt][2]=vertPoint[tt][2];
vPointNomal[tt1+tt][3]=vertPoint[tt][3];
vPointNomal[tt1+tt][4]=vertPoint[tt][4];
vPointNomal[tt1+tt][5]=vertPoint[tt][5];
}
// memcpy(vPointNomal+6*(m_vNumber-vertPos) ,vertPoint,sizeof(float)*6*vertPos);
//记录新生成的三角面片数据
w=0;
while (g_TriTable[cubeType][w]!=-1)
{
for(r=0;r<3;r++)
{
triIndex[facePos][r]=cellVerts[g_TriTable[cubeType][w++]];
if(triIndex[facePos][r]<0)
{
AfxMessageBox("有问题",MB_OK,0);
}
}
facePos++;
m_fNumber++;

} //end 对应while (g_TriTable[cubeType][w]!=-1)
//保存面数据
tt1=m_fNumber-facePos;
for(tt=0;tt<facePos;tt++)
{
pFace[tt1+tt][0]=triIndex[tt][0];
pFace[tt1+tt][1]=triIndex[tt][1];
pFace[tt1+tt][2]=triIndex[tt][2];
}
// memcpy(pFace+3*(m_fNumber-facePos)*sizeof(long),triIndex,sizeof(int)*3*facePos);
} memcpy(bottomXEdge,topXEdge,sizeof(short)*imageSize);
memcpy(bottomYEdge,topYEdge,sizeof(short)*imageSize);
memset(topXEdge,-1,sizeof(short)*imageSize);
memset(topYEdge,-1,sizeof(short)*imageSize);
memset(zEdge,-1,sizeof(short)*imageSize);
}
delete []tempSlice;
delete []bottomXEdge;
delete []bottomYEdge;
delete []topXEdge;
delete []topYEdge;
delete []zEdge;

return 1;
}
在OnDraw
glBegin(GL_TRIANGLES);
for(i=0;i<pDoc->m_fNumber;i++)
{
glNormal3fv(&(pDoc->vPointNomal[pDoc->pFace[ i][0]][3]));
glVertex3fv(&(pDoc->vPointNomal[pDoc->pFace[ i][0]][0]));

glNormal3fv(&(pDoc->vPointNomal[pDoc->pFace[ i ][1]][3]));
glVertex3fv(&(pDoc->vPointNomal[pDoc->pFace[ i ][1]][0]));

glNormal3fv(&(pDoc->vPointNomal[pDoc->pFace[ i ][2]][3]));
glVertex3fv(&(pDoc->vPointNomal[pDoc->pFace[ i ][2]][0]));
}

glEnd();
以上代码只用于理解，未测试

2. 《我的世界1.1》中地形生成的算法是什么

我来试着简要说一下地形生成吧。其实我一段时间也没怎么关心过地形生成，但是最近我在翻译一个开发文档的时候，那个开发文档提到了一些关于Minecraft地形生成的细节，所以我就被迫了解了一些关于Minecraft地形生成的知识。目前这里只介绍主世界正常情况下的生成，下界和末界或者超平坦什么的再说

?

3. 请问什么是MC算法如何在C++中应用这种算法

MC算法 -- 移动立方体算法
C++中调用OpenGL库是捷径。

微电子学与计算机 22卷9期 “基于改进MC算法的三维表面重建”：
MC算法是面绘制中构造等值面的方法中最具代表性的方法之一，已经得到了许多完善及改进。介绍了改进MC算法的实现步骤，并介绍了用Delphi7.0为开发平台，利用OpenGL三维图形软件包开发基于改进MC算法的三维表面重建系统的设计过程。

4. MCMC基本原理与应用（一）

贝叶斯统计推断是一个非常 艺术 的东西，他的先验结果（prior knowledge）是一种专家的经验，而 贝叶斯公式 给出的后验分布是贝叶斯推断的基本工具。每个大学生都应该学过数理统计，其中提及的贝叶斯统计推断都是一些基于简单先验和简单后验的结果。

但当我们一旦碰到复杂的先验（比如高维参数、复杂先验分布），我们用贝叶斯公式得出的后验分布将变得非常复杂，在计算上会非常困难。为此，先人提出了 MCMC算法 方便我们可以对任何后验分布进行计算或推断。其思想就是其名字：两个MC。

第一个MC: Monte Carlo(蒙特卡洛)。这个简单来说是让我们使用随机数（随机抽样）来解决计算问题。在MCMC中意味着：后验分布作为一个随机样本生成器，我们利用它来生成样本(simulation)，然后通过这些样本对一些感兴趣的计算问题（特征数，预测）进行估计。

第二个MC：Markov Chain(马尔科夫链)。第二个MC是这个方法的关键，因为我们在第一个MC中看到，我们需要利用 后验分布 生成随机样本，但后验分布太复杂，一些package中根本没有相应的随机数生成函数（如 rnorm(),rbinom() ）怎么办？答案是我们可以利用Markov Chain的 平稳分布 这个概念实现对复杂后验分布的抽样。

为了能顺利阐述MCMC算法，这里就简单地讲一下所涉及到的马尔科夫链概念。因为都是我个人的理解，这里所讲不敢说都是正确的，希望各位童鞋能够独立思考。

随机过程${X_n},X_n$的状态空间为有限集或可列集，比如当$X_n=i$,就称过程在n处于状态i。
定义：如果对任何一列状态$i_0,i_1,...,i_{n-1},i,j$，及对任何n≥0，随机过程${X_n}$满足Markov性质：
$$P{X_{n+1}=j|X_0=i_0,...X_{n-1}=i_{n-1},X_n=i}=P{X_{n+1}=j|X_n=i}$$

转移概率
$P{x_{n+1}=j|X_n=i}$成为Markov链的一步转移概率$P_{ij}^{n,n+1}$,当这个概率与n无关,则记之为$P_{ij}$

转移概率矩阵P
这个矩阵的元素就是一步转移概率$P_{ij}$

结论
一个Markov链可以由它的初始状态以及转移概率矩阵P完全确定。（证明略，自行网络或翻书）

所谓n步转移概率就是从状态i走n步正好到状态j的概率，我们记为$P_{ij}^{(n)}$。
利用概率分割的思想，由基础概率论中 全概率公式 可以得到$$P_{ij} ^{{(n)}=Sigma_{k=0}} {infty}P_{ik}P_{kj}^{(n-1)}$$ 写成矩阵形式就是$$P^{(n)}=P imes P^{(n-1)}$$
进一步推广，我们就推出了着名的Chapman-Kolmogorov方程：$$P_{ij} ^{{(n+m)}=Sigma_{k=0}} {infty}P_{ik} ^{(n)}P_{kj} {(m)}$$写成矩阵形式就是$$P ^{(m+n)}=P {(m)} imes P^{(n)}$$

现在我们已经有了n步转移概率的概念，一个很简单的想法就是如果n趋向无穷会怎么样？这个问题也是后面 极限分布 以及 平稳分布 的来源，更是MCMC算法中第二个MC的关键。
要回答这个问题首先要掌握几个关键概念，我先列出来，如果不熟悉的可以自行网络或翻书：

几个重要结论（证明略，自行网络，或者call me）

是时候回答上面那个问题了！（摘自网络共享PPT）

另外对于一个正常返、非周期的状态(也称为遍历ergodic)，我们有结论:$$对所有ileftrightarrow j,有limP_{ji} ^{{(n)}=limP_{ii}} {(n)}=frac{1}{mu_i}$$

什么是平稳分布？它和求极限概率分布有什么关系呢？

定义：Markov链有转移概率矩阵P，如果有一个概率分布${pi_i }满足pi_j =Sigma_{i=0}^{infty}pi_i P_{ij}$,则称为这个Markov链的平稳分布。这个定义用矩阵形式写出来就是π*P=π.

重要定理

若一个markov链中所有的状态都是遍历的，则对所有i,j有$limP_{ij} ^{{(n)}=pi_j存在，且pi={pi_j,j≥0}就是平稳分布$。反之，拖一个不可约Markov链只存在一个平稳分布，且这个Markov链的所有状态都是遍历的，则这个平稳分布就是该Markov链的极限概率分布。$$limP_{ij}} {(n)}=pi_j$$

5. 等值面的等值面生成算法

Cuberille方法
· Cuberrille等值面方法又称Opaque Cube算法，最初由Herman等人提出，后来又多次改进。算法主要分为两个步骤：
· (1) 确定边界单元
对于规则网格数据，其网格单元可看成是正六面体单元，整个三维数据就是由这种正六面体组成的，这种组成三维图象的基本正六面体单元称为体元。对于给定的阈值Ft，遍历体数据中的各个单元，将组成体元8个顶点上的值与Ft进行比较，找出顶点值跨越Ft的所有体元，即体元中有的顶点值大于阈值，有的顶点值小于阈值，因此体元内包含等值面片，这就是边界单元。
(2) 绘制各边界单元的6个多边形面，即将等值面看成是由各单元的六个外表面拼合而成
每个单元均为一正六面体，包括6个多边形面。对组成所有边界体元的多边形面进行绘制，即可产生最终的图象结果。在绘制多边形过程中应采用合适的光照模型和消隐技术。
如果在具有硬件深度缓存(Z-buffer)功能的计算机上运行立方体方法，可以将这组多边形不分次序地提交给硬件，由硬件完成消除隐藏面的任务。如果以软件方式执行立方体方法，在算法中必须考虑多边形的遮挡问题。一个有效的方法是把遍历体元集合与显示两个步骤合二为一，遍历体元集合时采用从后至前的次序。发现一个边界体元，就立刻显示它的6个面。后显示到屏幕上去的多边形将覆盖先显示的多边形，这样就达到了消除隐藏面的目的，这就是画家算法的思想。
Marching Cubes(MC)方法
Marching Cubes（移动立方体）方法是由W.E.Lorenson和H.E.Cline在1987年提出来的。由于这一方法原理简单，易于实现，目前已经得到了较为广泛的应用，成为三维数据等值面生成的经典算法，Marching Cubes算法又简称为MC算法。MC方法的原理如下：
在Marching Cubes方法中，假定原始数据是离散的三维空间规则数据，一个体元定义为由相邻层上的8个顶点组成的一个长方体。为了在三维数据中构造等值面，应先给定所求等值面的值，该方法的基本原理是逐个处理所有的体元，将体元各顶点处的值与给定的阈值进行比较，首先找出与等值面相交的体元，然后通过插值求等值面与体元棱边的交点，并将各交点连成三角形来构成等值面片，所有体元中的三角形集合就构成了等值面。由于这一方法是逐个处理所有的体元，因此被称为Marching Cubes方法。 在W.E.Lorenson和H.E.Cline于1987年提出Marching Cubes方法之后不久，他们就发现，当离散三维数据的密度很高时，由Marching Cubes方法在体元中产生的小三角面片常常很小，在图像空间中的投影面积与屏幕上一个像素点的大小差不多，甚至还要小，因此，通过插值来计算小三角面片是不必要的。随着新一代CT和MRI等设备的出现，二维切片中图象的分辨率不断提高，断层不断变薄，已经接近并超过计算机屏幕显示的分辨率。在这种情况下，常用于三维表面生成的Marching Cubes方法已不适用。于是，在1988年，仍由W.E.Lorenson和H.E.Cline两人提出了剖分立方体(Dividing Cubes)方法。

6. 关于西方经济学，LONG RUN里的MC 与ATC关系

关于西方经济学，LONG RUN里的MC 与ATC关系是MC交ATC于ATC最小的点上。
1、MC=dTC/dQ也即MC是TC对Q求导后得出的，由于ATC=AFC+AVC，而AFC也即平均固定成本是常数，求导后为0，所以MC=dAVC/dQ是平均可变成本的导数。
2、平均成本AC的最小点，也即极小值点，是它的一阶导数条件为0时取得，也即
dAC/dQ=d(ATC/Q)/dQ=(QATC'-ATC)/Q^2=0
QATC'=ATC
QMC=ATC
MC=ATC/Q=AVC
3、也即AVC取最小值的时候MC交于这个点。
补充资料：
MC边际成本指的是每一单位新增生产的产品（或者购买的产品）带来的总成本的增量。 这个概念表明每一单位的产品的成本与总产品量有关。ATC平均成本是指一定范围和一定时期内成本耗费的平均水平。平均成本总是针对一定的产品或劳务而言的。一定时期产品生产或劳务提供平均成本的变化，往往反映了一定范围内成本管理总体水平的变化。
(6)mc算法实现扩展阅读：
一、什么是西方经济学
西方经济学是指流行于西欧北美资本主义发达国家的经济理论和政策主张，它是15世纪西方经济学产生，18世纪西方经济学建产以来，特别是19世纪70年代以后一直到目前为止认为是能够说明当代资本主义市场经济运行和国家调节的重要理论、概念、政策主张和分析方法进行了综合和系统化形成的。被称为“社会科学之王”。另外，《西方经济学》是我国高等院校财经类和管理类专业必开的一门专业基础课。
二、西方经济学学科分类
1、微观经济学-- 研究家庭、厂商和市场合理配置经济资源的科学 -- 以单个经济单位的经济行为为对象；以资源的合理配置为解决的主要问题；以价格理论为中心理论；以个量分析为方法；其基本假定是市场出清、完全理性、充分信息。
2、宏观经济学- 研究国民经济的整体运行中充分利用经济资源的科学 — 以国民经济整体的运行为对象；以资源的充分利用为解决的主要问题；以收入理论为中心理论；以总量分析为方法，其基本假定为市场失灵、政府有效。

7. 强化学习基础篇（十八）TD与MC方法的对立统一

前面介绍TD的过程中，我们已经提到过一些TD和MC的区别，例如在Bias与Variance角度看：

a. MC具有高方差，为无偏估计

b. TD具有低方差，为有偏估计

a. TD可以在知道最终结果之前就进行学习。

b. TD也可以在没有最终输出的场景下进行学习。

假设只有有限的经验，比如10幕数据或100个时间步。在这种情况下，使用增量学习方法的一般方式是反复地呈现这些经验，直到方法最后收敛到一个答案为止。给定近似价值函数 ，在访问非终止状态的每个时刻 ，使用下面两式计算相应的增量但是价值函数仅根据所有增量的和改变一次。

然后，利用新的值函数再次处理所有可用的经验，产生新的总增量，依此类推，直到价值函数收敛。我们称这种方法为批量更新，因为只有在处理了整批的训练数据后才进行更新。

在批量更新下，如果经验趋于无穷多，只要选择足够小的步长参数 ， 就能确定地收敛到与 无关的唯一结果。常数 MC方法在相同条件下也能确定地收敛，但是会收敛到不同的结果。

当然，这是在经验趋于无穷（也即无数次试验）的情况下达到的理想情况，但是实际中我们不可能达到，那如果我们利用有限的经验来对值函数进行估计将得到什么样的结果？比方说，对于下面这 个 episode：

如果我们重复从这 个episode中进行采样，对于某一次采样得到的样本 应用MC或者 方法，会得到什么样的结论呢？先来看一个例子。

假设在一个强化学习问题中有A和B两个状态，模型未知，不涉及策略和行为，只涉及状态转换和即时奖励，衰减系数为1。现有如下表所示8个完整状态序列的经历，其中除了第1个状态序列发生了状态转移外，其余7个完整的状态序列均只有一个状态构成。现要求根据现有信息计算状态A、 B的价值分别是多少？

考虑分别使用MC算法和TD算法来计算状态A、 B的价值：

对于MC算法：

在8个完整的状态序列中，只有第一个序列中包含状态A，因此A价值仅能通过第一个序列来计算：

所以可以得到 。

状态B的价值，则需要通过状态B在8个序列中的收获值来平均：

可以得到 。

对于TD算法 ，

再来考虑应用TD算法。TD算法试图利用现有的episode经验构建一个MDP（如下图），由于存在一个episode使得状态A有后继状态B，因此状态A的价值是通过状态B的价值来计算的，同时经验表明A到B的转移概率是100%，且A状态的即时奖励是0，并且没有衰减，因此A的状态价值等于B的状态价值。

其计算过程如下：

因此在TD算法下 。

AB Example体现了通过批量 和批量蒙特卡洛方法计算得到的估计值之间的差别。批量蒙特卡洛方法总是找出最小化训练集上均方误差的估计，而批量 总是找出完全符合马尔科夫过程模型的最大似然估计参数。一个参数的最大似然估计是使得生成训练数据的概率最大的参数值。

TD与MC的另一个差异

现在为止所阐述的MC学习算法、TD学习算法和DP算法都可以用来计算状态价值。它们的特点也是十分鲜明的，MC和TD是两种在不依赖模型的情况下的常用方法，这其中又以MC学习需要完整的状态序列来更新状态价值，TD学习则不需要完整的状态序列；DP算法则是基于模型的计算状态价值的方法，它通过计算一个状态 所有可能的转移状态 及其转移概率以及对应的即时奖励来计算这个状态 的价值。

下图，非常直观的体现了三种算法的区别。

综合上述三种学习方法的特点，可以小结如下：

当使用单个采样，同时不经历完整的状态序列更新价值的算法是TD学习； 当使用单个采样，但依赖完整状态序列的算法是MC学习； 当考虑全宽度采样，但对每一个采样经历只考虑后续一个状态时的算法是DP学习； 如果既考虑所有状态转移的可能性，同时又依赖完整状态序列的，那么这种算法是穷举（exhausive search）法。

需要说明的是：DP利用的是整个MDP问题的模型，也就是状态转移概率，虽然它并不实际利用采样经历，但它利用了整个模型的规律，因此也被认为是全宽度（full width） 采样的。

8. 我的世界生态环境算法怎么实现的

用变量，测试生态环境

9. 该算法是怎么实现的（mc[prop]=((mc.y+2*disy)/disy-1)/2*(n2-n1)+n1）

prop:String,n1:Number=.5, n2:Number=1):void
{