三维骨架提取Main.docx
《三维骨架提取Main.docx》由会员分享，可在线阅读，更多相关《三维骨架提取Main.docx（19页珍藏版）》请在冰豆网上搜索。
三维骨架提取Main
使用广义势场法计算三维骨架
主要结构体
typedefstruct//点类型（立体的点坐标）
{
shortx;
shorty;
shortz;
}VoxelPosition;
typedefstruct//向量类型
{
doublexd;
doubleyd;
doublezd;
}Vector;
enumCriticalPointType//关键点的类型
{
CPT_SADDLE=1,//1鞍点
CPT_ATTRACTING_NODE,//2引力点
CPT_REPELLING_NODE,//3斥力点
CPT_UNKNOWN//4其他
};
structVoxelPositionDouble//点
{
doublex;
doubley;
doublez;
};
structCriticalPoint
{
VoxelPositionDoubleposition;//关键点的坐标
CriticalPointTypetype;//关键点的类型
Vectorevect[3];//3个特征向量数组
doubleeval[3];//3个特征值
};
structSkeleton
{
VoxelPositionDouble*Points;//存储骨架上的点
intsizePoints;//初始化时分配的最大储存空间本程序是50000
intnumPoints;//骨架上的点的个数
int**Segments;//骨架段为二维数组
//numSegments行4列的数组LEFTRIGHTFIRSTLAST
intsizeSegments;//初始化分配的最大存储空间，本程序是5000
intnumSegments;//骨架段的个数
};
骨架段的说明：
LEFT
RIGHT
FIRST
LAST
第一段
第一个骨架点index
最后一个骨架点的index
第二段
第三段
……
立体长宽高：
L=123M=116N=114;
计算势场用到的距离的次方数fieldStrenght=8;
计算散度值时和阈值有关的percHDPts=3;
vfin=false;
vfout=false;
interactive=false;
//这些变量代表的实际意义不清楚
#defineSURF100//surfacevoxel
#defineBOUNDARY110//boundaryvoxel-participatesinpotentialfield//calculation
#defineINTERIOR200//interiorvoxel
#definePADDING_MIN210//addedvoxelsinordertothicktheobject
#defineNR_PAD_VALUES40//areinthisrange:
PADDING_MINtoPADDING_MIN+NR_PAD_VALUES
#defineEXTERIOR0//background（exteriortotheobject）voxel（air）
貌似INTERIOR点是立体的体素点，EXTERIOR点是空气，
SURF也是立体的点，是暴露在空气中的点？
？
？
具体算法：
定义了Skeleton变量Skel，
函数AllocateSkeleton对它进行初始化
为Skel分配50000个点5000个段
VoxelPositionDouble*Points;//初始分配50000个空间
intsizePoints;//点最多50000
intnumPoints;//0
Skel:
int**Segments;//初始分配5000指针地址
intsizeSegments;//段最大为5000
intnumSegments;//0
计算骨架：
调用函数boolpfSkel（
char*volFileName,//[in]立体文件的地址
intL,intM,intN,//[in]立体的长宽高
intdistCharges,//[in]点电荷距离物体边界的距离
intfieldStrenght,//[in]势场强度介于4-9，程序为8
floatpHDPts,//[in]高散度点的比例
Skeleton*Skel,//[out]包含骨架点的数组
cbChangeParameterspfnChangeParams,//[in]回调函数，改变参数用
void*other,//[in]回调函数的地址，输出骨架地址
char*vfInFile,//[in]向量场输入文件程序为NULL
char*vfOutFile//[in]向量场输出文件程序为NULL
）
unsignedchar*f;//存储立体的点值
检查distCharges>=0,以及fieldStrenght是否介于1和10之间
读取立体的文件ReadVolume（volFileName,L,M,N,&f）
f分配了L*M*N的空间，存储了立体的值
问：
立体的值代表什么含义
确保立体在三个方向填充了足够的空的平面，所以在距离边界特定的位置放置点电荷
CheckVolumePadding（&f,&L,&M,&N,distCharges）问：
这个函数到底做什么用的
这个函数首先是对f中的值划分为两类
f中的值
0INTERIOR200
非0和i=0,L-1j=0,M-1k=0,N-1EXTERIOR0外部点是指外面的点吗
求出点值为EXTERIOR的点的坐标范围
minxmaxx
minymaxy
minzmaxz
满足以下条件则返回真问：
这个是什么意思，为什么这么处理
minx-distCharges>=0
miny-distCharges>=0
minz-distCharges>=0
maxx+distCharges<=L
maxy+distCharges<=M
maxz+distCharges<=N
下面执行的函数是，填充立体
作用是Makesurethevolumedoesnothaveholesinit问：
这个函数的作用又是什么
MakeSolidVolume（L,M,N,f,EXTERIOR,INTERIOR）;
遍历所有的点，将不是EXTERIOR的点，也就是全部的INTERIOR
设置为OBJECT_VAL=1
for（i=0;i{
FloodFillZPlane（i,L,M,N,vol）;
}
用OUTSIDE_1=2填充每个Z面
//判断该点周围的4个点，如果是0（EXTERIOR）的话，设置为OUTSIDE_1=2
//↑
//←→
//↓
填充函数，遍历XOY平面，设置bool变量anyChange开始为false
从原点开始遍历，原点的值设置为OUTSIDE_1=2次序为左→右，上→下，如果4点中出现了EXTERIOR，将EXTERIOR点设置为OUTSIDE_1=2，设置anyChange开始为true，否则不变，如果anychange的值为true，则从右下点开始，遍历次序为右→左，下→上，设置同上。
接下来填充Y平面，函数算法同Z，对值是EXTERIOR=0和OUTSIDE_1=2的值设置为OUTSIDE_2=3
for（i=0;i{
FloodFillYPlane（i,L,M,N,vol）;
}
接下来填充X平面，函数算法同Z，对值是EXTERIOR=0和OUTSIDE_1=2和OUTSIDE_2=3的值设置为OUTSIDE_3=4
for（i=0;i{
FloodFillXPlane（i,L,M,N,vol）;
}
下边是设置立体的值，分成三类SURFINTERIOREXTERIOR
boolFlagVolume（unsignedchar*vol,intL,intM,intN）
1遍历所有的点（L*M*N），将点的值为0的设置为EXTERIOR=0，非0的设置为INTERIOR=200
2遍历所有的INTERIOR点（假设B点），如果该点的邻接有EXTERIOR点，则将B点设置为SURF，这里的邻接点是B点的上下左右前后六个点
计算电势场
force=newVector[L*M*N];//存储每个点的电势场
//thickentheobjectwithlayersofextravoxelsand
//computethepotentialfield
将点电荷放置在物体的外边
ExpandNLayers（L,M,N,f,distCharges）;
检查distCharges的值不能超过NR_PAD_VALUES=40
IfdistCharges==0不做任何操作return//程序中设定的是0
IfdistCharges>0做以下操作，在立体的SURF外面加上distCharge层，每一层的值依次加1。
ExpandVolume（L,M,N,f,SURF,PADDING_MIN）;
for（i=1;i{
ExpandVolume（L,M,N,f,PADDING_MIN+（i-1）,PADDING_MIN+i）;
}
先遍历SURF点，假如SURF点的周围（6邻接点）是EXTERIOR点（点B）
把点B值设置为PADDING_MIN
同理对于下边的for循环是依次在立体的外边加上一层，同样是对EXTERIOR点设置为新的值。
IfdistCharges<0则剥掉nrLayers层，可能会造成立体的断开
for（i=0;i>nrLayers;i--）
{
PeelVolume（f,L,M,N）;
}
PeelVolume（f,L,M,N）;是将SURF点变成EXTERIOR点，再求出新的SURF点
计算势场的函数
CalculatePotentialField（L,M,N,f,fieldStrenght,force）;//fieldStrenght=8
函数
1先CheckVolumePadding（f,L,M,N）//fastcheck
//Makesurethevolumeispaddedwithatlease1emptyplaneinall3directions
确保立体的四周是由EXTERIOR点填充的
2VoxelPosition*Bound;
Bound=newVoxelPosition[BOUND_SIZE]）BOUND_SIZE=1200000
遍历立体点（1~L-1,1~M-1,1~N-1）,将与EXTERIOR点邻接的INTERIOR点赋为SURF，接着把SURF点赋为BOUND=110，点的坐标存到BOUND中去。
3对边界点进行排序，顺序为ZYX
SortBoundaryArray（numBound,Bound）;
对Z排序:
SortByZ（0,numBound-1,Bound）;//Z次序为主次序，排列之后点是按照Z由小到大的顺序
对Y排序将Z值相同的划分为一组，组内对Y由小到大排序SortByY（st,i-1,Bound）;
对X排序，对于ZY值均相同的划分为一组，组内对X由小到大排序SortByX（st,i-1,Bound）;
4开始计算每个点的势场值（势场是由BOUNDARY点产生的）
只计算不是SURFBOUNDEXTERIOR的点处的值
PF_THRESHOLD=100
假设要计算A点处的势场值，|BOUNDARY.x-A.x|<=PF_THRESHOLD
|BOUNDARY.y-A.y|<=PF_THRESHOLD
|BOUNDARY.z-A.z|<=PF_THRESHOLD
满足这些条件的BOUNDARY点才对A点的势场值起作用，其他的忽略
问：
这种筛选方法有问题吗
计算公式是
for（k=0;k{
根据PF_THRESHOLD得到zStartIndexzEndIndex的值
for（j=0;j{
得到yStartIndex，yEndIndex的值
for（i=0;i{
idx++;
if（f[idx]==0）continue;
if（f[idx]==SURF）continue;
if（f[idx]==BOUNDARY）continue;
得到startIndex，endIndex的值
for（s=startIndex;s<=endIndex;s++）//计算该点处的势场值
{
v1=i-Bound[s].x;
v2=j-Bound[s].y;
v3=k-Bound[s].z;
r=sqrt（v1*v1+v2*v2+v3*v3）;//欧几里得度量
force[idx].xd=force[idx].xd+（v1/r^fieldStrenght）;//fieldStrenght=8
force[idx].yd=force[idx].yd+（v2/r^fieldStrenght）;force[idx].zd=force[idx].zd+（v3/r^fieldStrenght）;
｝
｝
｝
｝
force向量场归一化
下面计算SURF和BOUNDARY点处的电势值，通过点的26邻接点的电势值的平均值来确定，但是26个点中的SURFBOUNDARYEXTERIOR点忽略
得到关键点
GetCriticalPoints（force,L,M,N,f,&CritPts,&numCritPts）;CritPts存储得到的关键点
为关键点数组分配2000的空间
（*CritPts）=newCriticalPoint[MAX_NUM_CRITPTS]）MAX_NUM_CRITPTS=2000
（1~L-1,1~M-1,1~N-1）检查每个点（假设B点）的邻接点，如果邻接点中出现了SURF，BOUNDARY，EXTERIOR的情况，跳过点B，continue
问：
为什么选取这几个邻接点
inds[0]=idx;//B点
inds[1]=idx+1;//以下为B的邻接点
inds[2]=idx+L;
inds[3]=idx+L+1;
inds[4]=idx+slsz;
inds[5]=idx+slsz+1;
inds[6]=idx+slsz+L;
inds[7]=idx+slsz+L+1;
判断点（i,j,k）是否是关键点
FindCriticalPointInIntCell（i,j,k,inds,L,M,N,ForceField,&critPt）
A--对于势场值为0的if（（IS_ZERO（ForceField[inds[0]].xd））&&（IS_ZERO（ForceField[inds[0]].yd））&&（IS_ZERO（ForceField[inds[0]].zd）））
则将点添加到关键点。
B—如果该点和7个邻接点有正负变号if（ChangeInSign（ForceField,inds,8）），则该点有可能成为关键点，将这个cell分成8个subcells（二阶魔方）
for（kk=0;kk<2;kk++）
{
for（jj=0;jj<2;jj++）
{
for（ii=0;ii<2;ii++）
{
if（FindCriticalPointInFloatCell（x+ii/2.00,y+jj/2.00,z+kk/2.00,0.50,sX,sY,sZ,ForceField,critPt））
{
returntrue;
}
}
}
}
判断变化的函数ChangeInSign，
forces[indices[0]].xd的符号与其他7个邻接点forces[indices[1-7]].xd的符号比较，如果不同则count++，如果比较完没有相反的，但是forces[indices[0]].xd=0，count++
forces[indices[0]].yd和forces[indices[0]].zd同样这样处理，最后得到的count>=3返回true
如果有方向的改变，则执行下列操作
将（x,y,z）→（x+1,y+1,z+1）的立体划分成8块，每个立体用最靠近原点的那个点表示
for（kk=0;kk<2;kk++）
{
for（jj=0;jj<2;jj++）
{
for（ii=0;ii<2;ii++）
{
if（
FindCriticalPointInFloatCell（
x+ii/2.00,y+jj/2.00,z+kk/2.00,0.50,
sX,sY,sZ,ForceField,critPt））
{
returntrue;
}
}
}
}
boolFindCriticalPointInFloatCell
（doublex,doubley,doublez,doublecellsize,intsX,intsY,intsZ,Vector*ForceField,VoxelPositionDouble*critPt
）这个函数判断一个cell是否是关键点，这里的（x,y,z）表示每个cell的最靠近原点的点，用插值的方法得到立体的每个顶点的ForceField值这里cellsize的值是0.5
cv[0]→（x,y,z）
cv[1]→（x+cellsize,y,z）
cv[2]→（x,y+cellsize,z）
cv[3]→（x+cellsize,y+cellsize,z）
cv[4]→（x,y,z+cellsize）
cv[5]→（x+cellsize,y,z+cellsize）
cv[6]→（x,y+cellsize,z+cellsize）
cv[7]→（x+cellsize,y+cellsize,z+cellsize）
判断cv[0]的IS_ZERO（cv[0].xd））&&（IS_ZERO（cv[0].yd））&&（IS_ZERO（cv[0].zd）））
满足则返回true
如果不是则判断if（ChangeInSign（cv,8））如果有变号，则这个cell是候选的
判断cell是否足够小if（cellsize<=（1.00/1048576）），则将其cell的中心作为关键点
（*critPt）.x=x+（cellsize/2.00）;
（*critPt）.y=y+（cellsize/2.00）;
（*critPt）.z=z+（cellsize/2.00）;返回true
如果这个是一个候选的cell，但是cellsize不够小，则再将cell划分为更小的8个cell。
。
。
。
每个划分后的cell再执行FindCriticalPointInFloatCell这个函数
得到关键点后，再确定关键点的类型
vdist=1.00/1048576
在每个关键点的上下左右前后插入六个点
//根据插入点的位置来得到插入到该点的向量的大小，
cv[0]→（x+vdist,y,z）
cv[1]→（x-vdist,y,z）
cv[2]→（x,y+dist,z）
cv[3]→（x,y-vdist,z）
cv[4]→（x,y,z+vdist）
cv[5]→（x,y,z-vdist）
构造雅克比行列式
TNT:
:
Array2DJac（3,3,0.00）;二维数组放雅克比行列式
TNT:
:
Array2DEigVals（3,3,0.0）;三维数组放特征值
TNT:
:
Array2DEigVects（3,3,0.0）;存放特征向量
下边计算出雅克比行列式的特征值和特征向量
JAMA:
:
Eigenvalueev（Jac）;ev.getD（EigVals）;ev.getV（EigVects）;
EigVals存放的是特征值，举例
EigVects存放的是对应的特征向量
根据得到的特征值和特征向量判断出关键点的类型
对特征向量进行归一化
存储每个关键点的特征值和特征向量到（*CritPts）[i]中，
得到关键点的坐标类型特征值特征向量后，下边是得到第一骨架
GetLevel1Skeleton（force,L,M,N,CritPts,numCritPts,Skel）;
定义一些相关变量
//该数组用来标记某个关键点是否已经是骨架的一部分了-1表示还不是
critPtInSkeleton=newint[numCritPts]，初始化为-1
从鞍点出发，正特征值对应的特征向量为方向，直到下一点是生成骨架的点或者是关键点
for（i=0;i
升级会员 