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
{