三维骨架提取Main

1、使用广义势场法计算三维骨架主要结构体typedef struct /点类型（立体的点坐标）short x;short y;short z; VoxelPosition;typedef struct /向量类型 double xd; double yd;double zd; Vector;enum CriticalPointType /关键点的类型 CPT_SADDLE = 1, / 1 鞍点 CPT_ATTRACTING_NODE, /2 引力点CPT_REPELLING_NODE, /3 斥力点CPT_UNKNOWN /4 其他;struct VoxelPositionDouble /点do

2、uble x;double y;double z;struct CriticalPoint VoxelPositionDoubleposition; /关键点的坐标CriticalPointTypetype; /关键点的类型 Vector evect3; /3个特征向量数组 double eval3; /3个特征值 ;struct Skeleton VoxelPositionDouble *Points; /存储骨架上的点 int sizePoints; /初始化时分配的最大储存空间本程序是50000 int numPoints; /骨架上的点的个数 int * Segments; /骨架段为

3、二维数组 /numSegments行 4列的数组 LEFT RIGHT FIRST LAST int sizeSegments; /初始化分配的最大存储空间，本程序是5000 int numSegments; /骨架段的个数;骨架段的说明：LEFTRIGHTFIRSTLAST第一段第一个骨架点index最后一个骨架点的index第二段第三段立体长宽高： L = 123 M = 116 N = 114; 计算势场用到的距离的次方数 fieldStrenght = 8; 计算散度值时和阈值有关的 percHDPts= 3; vfin = false; vfout = false; interact

4、ive = false;/这些变量代表的实际意义不清楚#define SURF 100/ surface voxel#define BOUNDARY 110/ boundary voxel - participates in potential field /calculation #define INTERIOR 200/ interior voxel#define PADDING_MIN210/ added voxels in order to thick the object#define NR_PAD_VALUES 40/ are in this range: PADDING_MIN

5、to PADDING_MIN + NR_PAD_VALUES#define EXTERIOR 0/ background (exterior to the object) voxel (air)貌似 INTERIOR点是立体的体素点，EXTERIOR点是空气，SURF也是立体的点，是暴露在空气中的点？具体算法：定义了Skeleton变量Skel，函数AllocateSkeleton对它进行初始化为Skel分配50000个点 5000个段 VoxelPositionDouble *Points; /初始分配50000个空间 int sizePoints; /点最多50000 int numPoi

6、nts; /0Skel: int * Segments; /初始分配5000指针地址 int sizeSegments; /段最大为5000 int numSegments; /0计算骨架：调用函数bool pfSkel(char* volFileName, / in 立体文件的地址int L, int M, int N, / in立体的长宽高int distCharges, / in 点电荷距离物体边界的距离int fieldStrenght, / in 势场强度 介于4-9，程序为8float pHDPts, / in高散度点的比例 Skeleton *Skel, / out 包含骨架点的

7、数组 cbChangeParameters pfnChangeParams, / in 回调函数，改变参数用 void *other, / in 回调函数的地址，输出骨架地址 char *vfInFile, / in 向量场输入文件 程序为NULLchar *vfOutFile / in 向量场输出文件 程序为NULL)unsigned char *f; /存储立体的点值检查distCharges=0,以及fieldStrenght是否介于1和10之间读取立体的文件 ReadVolume(volFileName, L, M, N, &f)f分配了L*M*N的空间，存储了立体的值问：立体的值代表

8、什么含义确保立体在三个方向填充了足够的空的平面，所以在距离边界特定的位置放置点电荷CheckVolumePadding(&f, &L, &M, &N, distCharges) 问：这个函数到底做什么用的这个函数首先是对f中的值划分为两类 f中的值0 INTERIOR 200非0和i=0,L-1 j=0,M-1 k=0,N-1 EXTERIOR 0 外部点是指外面的点吗求出 点值为EXTERIOR的点的坐标范围 minx maxx miny maxy minz maxz 满足以下条件则返回真 问：这个是什么意思，为什么这么处理minx - distCharges= 0 miny - distC

9、harges = 0 minz - distCharges = 0 maxx + distCharges = L maxy + distCharges = M maxz + distCharges= N 下面执行的函数是，填充立体作用是Make sure the volume does not have holes in it 问：这个函数的作用又是什么 MakeSolidVolume(L, M, N, f, EXTERIOR, INTERIOR);遍历所有的点，将不是EXTERIOR的点，也就是全部的INTERIOR设置为OBJECT_VAL=1for(i=0; i N; i+) Flood

10、FillZPlane(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，对值是EX

11、TERIOR=0和OUTSIDE_1=2的值设置为OUTSIDE_2=3 for(i=0; i M; i+)FloodFillYPlane(i, L, M, N, vol);接下来填充X平面，函数算法同Z，对值是EXTERIOR=0和OUTSIDE_1=2 和OUTSIDE_2=3的值设置为OUTSIDE_3=4for(i=0; i L; i+) FloodFillXPlane(i, L, M, N, vol);下边是设置立体的值，分成三类SURF INTERIOR EXTERIORbool FlagVolume(unsigned char *vol, int L, int M, int N)

12、1遍历所有的点（L*M*N），将点的值为0的设置为EXTERIOR=0，非0的设置为INTERIOR=2002 遍历所有的INTERIOR点（假设B点），如果该点的邻接有EXTERIOR点，则将B点设置为SURF，这里的邻接点是 B点的上下左右前后六个点计算电势场force = new VectorL*M*N;/ 存储每个点的电势场 / thicken the object with layers of extra voxels and / compute the potential field将点电荷放置在物体的外边ExpandNLayers(L, M, N, f, distCharges)

13、;检查distCharges的值不能超过NR_PAD_VALUES=40If distCharges=0 不做任何操作 return /程序中设定的是0If distCharges0 做以下操作，在立体的SURF外面加上distCharge层，每一层的值依次加1。ExpandVolume(L, M, N, f, SURF, PADDING_MIN); for(i=1; i nrLayers; i+) ExpandVolume(L, M, N, f, PADDING_MIN + (i - 1), PADDING_MIN + i);先遍历SURF点，假如SURF点的周围（6邻接点）是EXTERIO

14、R点（点B）把点B值设置为PADDING_MIN同理对于下边的for循环是依次在立体的外边加上一层，同样是对EXTERIOR点设置为新的值。If distCharges 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) /fast ch

15、eck / Make sure the volume is padded with at lease 1 empty plane in all 3 directions 确保立体的四周是由EXTERIOR点填充的2 VoxelPosition* Bound; Bound = new VoxelPositionBOUND_SIZE) BOUND_SIZE=1200000遍历立体点(1L-1,1M-1,1N-1),将与EXTERIOR点邻接的INTERIOR点赋为SURF，接着把SURF点赋为BOUND=110，点的坐标存到BOUND中去。3 对边界点进行排序，顺序为ZYXSortBoundary

16、Array(numBound, Bound);对Z排序:SortByZ(0, numBound-1, Bound); /Z次序为主次序，排列之后点是按照Z由小到大的顺序对Y排序将Z值相同的划分为一组，组内对Y由小到大排序SortByY(st, i-1, Bound);对X排序，对于Z Y值均相同的划分为一组，组内对X由小到大排序SortByX(st, i-1, Bound);4 开始计算每个点的势场值（势场是由BOUNDARY点产生的）只计算不是 SURF BOUNDEXTERIOR的点处的值 PF_THRESHOLD=100假设要计算A点处的势场值， | BOUNDARY.x-A.x|= P

17、F_THRESHOLD | BOUNDARY.y-A.y|= PF_THRESHOLD | BOUNDARY.z-A.z|= PF_THRESHOLD 满足这些条件的BOUNDARY点才对A点的势场值起作用，其他的忽略问：这种筛选方法有问题吗计算公式是for (k = 0; k N; k+) /遍历所有的点根据PF_THRESHOLD 得到zStartIndex zEndIndex的值 for (j = 0; j M; j+) 得到yStartIndex ，yEndIndex 的值for (i = 0; i L; i+) idx+; if(fidx = 0) continue; if(fidx

18、 = SURF) continue;if(fidx = BOUNDARY) continue;得到startIndex，endIndex的值 for (s = startIndex; s = endIndex; s+) /计算该点处的势场值 v1 = i - Bounds.x;v2 = j - Bounds.y;v3 = k - Bounds.z;r = sqrt(v1*v1 + v2*v2 + v3*v3);/欧几里得度量 forceidx.xd = forceidx.xd + (v1 / r fieldStrenght); / fieldStrenght=8forceidx.yd = fo

19、rceidx.yd + (v2 / r fieldStrenght); forceidx.zd = forceidx.zd + (v3/ r fieldStrenght);force向量场归一化 下面计算SURF和BOUNDARY点处的电势值，通过点的26邻接点的电势值的平均值来确定，但是26个点中的SURF BOUNDARY EXTERIOR点忽略得到关键点 GetCriticalPoints(force, L, M, N, f, &CritPts, &numCritPts); CritPts存储得到的关键点为关键点数组分配2000的空间(*CritPts)=new CriticalPoin

20、tMAX_NUM_CRITPTS) MAX_NUM_CRITPTS=2000(1L-1,1M-1,1N-1)检查每个点（假设B点）的邻接点，如果邻接点中出现了SURF，BOUNDARY，EXTERIOR的情况，跳过点B，continue问：为什么选取这几个邻接点inds0 = idx; /B点inds1 = idx + 1;/以下为B的邻接点inds2 = idx + L;inds3 = idx + L + 1;inds4 = idx + slsz;inds5 = idx + slsz + 1;inds6 = idx + slsz + L;inds7 = idx + slsz + L + 1;

21、判断点(i,j,k)是否是关键点FindCriticalPointInIntCell(i, j, k, inds, L, M, N, ForceField, &critPt)A-对于势场值为0的if(IS_ZERO(ForceFieldinds0.xd) &(IS_ZERO(ForceFieldinds0.yd) &(IS_ZERO(ForceFieldinds0.zd)则将点添加到关键点。B如果该点和7个邻接点有正负变号if(ChangeInSign(ForceField, inds, 8)，则该点有可能成为关键点，将这个cell分成8个subcells（二阶魔方）for(kk=0; kk

22、2; kk+)for(jj=0; jj 2; jj+) for(ii=0; ii =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)return true;bool Find

23、CriticalPointInFloatCell(double x, double y, double z, double cellsize,int sX, int sY, int sZ, Vector* ForceField, VoxelPositionDouble *critPt) 这个函数判断一个cell是否是关键点，这里的(x,y,z)表示每个cell的最靠近原点的点，用插值的方法得到立体的每个顶点的ForceField值 这里cellsize的值是0.5cv0 （x, y, z）cv1 (x+cellsize, y, z)cv2 (x,y+cellsize,z)cv3 (x+cell

24、size,y+cellsize,z)cv4 (x,y,z+cellsize)cv5 (x+cellsize,y,z+cellsize)cv6 (x,y+cellsize,z+cellsize)cv7 (x+cellsize,y+cellsize,z+cellsize)判断cv0的IS_ZERO(cv0.xd) &(IS_ZERO(cv0.yd) &(IS_ZERO(cv0.zd) 满足 则返回true如果不是 则判断if(ChangeInSign(cv, 8) 如果有变号，则这个cell是候选的 判断cell是否足够小if(cellsize = (1.00 / 1048576)，则将其cell

25、的中心作为关键点 (*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在每个关键点的上下左右前后插入六个点/根据插入点的位置来得到插入到该点的向量的大小

26、， cv0(x+vdist, y, z)cv1(x-vdist,y,z)cv2(x,y+dist,z)cv3(x,y-vdist,z)cv4(x,y,z+vdist)cv5(x,y,z-vdist)构造雅克比行列式 TNT:Array2D Jac(3, 3, 0.00); 二维数组放雅克比行列式 TNT:Array2DEigVals(3, 3, 0.0); 三维数组放特征值 TNT:Array2DEigVects(3, 3, 0.0); 存放特征向量下边计算出 雅克比行列式的特征值和特征向量JAMA:Eigenvalueev(Jac); ev.getD(EigVals); ev.getV(Ei

27、gVects);EigVals存放的是特征值，举例EigVects 存放的是对应的特征向量根据得到的特征值和特征向量判断出关键点的类型对特征向量进行归一化存储每个关键点的特征值和特征向量到(*CritPts)i中，得到关键点的坐标 类型 特征值 特征向量后，下边是得到第一骨架GetLevel1Skeleton(force, L, M, N, CritPts, numCritPts, Skel);定义一些相关变量/该数组用来标记某个关键点是否已经是骨架的一部分了 -1表示还不是critPtInSkeleton = new intnumCritPts ，初始化为-1从鞍点出发，正特征值对应的特征向

28、量为方向，直到下一点是生成骨架的点或者是关键点for(i = 0; i numCritPts; i+) 该点不是鞍点，则continue; for(j=0; j 0) /特征值为正的 沿着特征向量的方向FollowStreamlines(); 沿着特征向量的反方向 FollowStreamlines(); 进入函数FollowStreamlines(i, true,&(CritPtsi.position),&whereTo,L, M, N, ForceField, CritPts, numCritPts,critPtInSkeleton,Skel); 函数具体实现步骤：遍历所有的鞍点，每条骨架

29、段都是从鞍点开始的，第一次循环沿着方向向量作为方向得到下一点的坐标，接下来的循环将使用当前点（x,y,z）的势场值作为方向（注：（x,y,z）可能不是整形的坐标，采用插值来求得其势场向量值），循环结束的条件很多1. step过多，认为是Most likely we are chasing our own tail.2. 遇到了骨架A中的点a，则将该点加入到当前段，并且作为段结束的端点，将A段在点a处分开，A段变成了2段3. 遇到了关键点，则直接加入当前段，并且作为段结束的端点。4. 计算出来的nextpos和startpos值是相等的，不做处理，结束5. 段到了立体的边界，将当前段删除以下是相

30、关的函数s = AddSkeletonPoint(Skel, &Startpos); 函数参数为 骨架，待添加关键点的坐标 返回值为新添加的点的索引（0开始）将点Startpos添加到Skel的Points数组中添加新的骨架段crtSegment = AddNewSkelSegment(Skel, s); 函数参数s为该点在Skel-Point的索引，返回值是crtSegment为段Skel-Segments的索引（0开始）自动添加新的段(Skel-SegmentsSkel-numSegments = new int4初始化为 点s的索引LEFTRIGHTFIRSTLAST第一段第二段第三段c

31、rtSegmentssss函数sp=CloseToSkel(&Startpos, oCPSkelPos,Skel, nrAlreadyInSkel, SP_SMALL_DISTANCE, crtSegLength); SP_SMALL_DISTANCE=0.5判断该点是否和一个骨架中的点离得很近（也就是遇到了骨架点，可能骨架的当前段在该点结束），如果是的话,返回骨架点的位置(inSkel） i ，特别的，假如骨架点不是端点，但是和端点离得近，这里不返回骨架点的位置，返回-1 wait ait to get to one of the segment end points不是的话，返回-1详细：

32、遍历所有的骨架点，计算距离for(i=0; i nrAlreadyInSkel; i+) if( |x1-x|+|y1-y|+|z1-z| SP_SMALL_DISTANCE ) seg = GetSkelSegment(Skel, i, &endpoint);/得到这个点i在骨架中的段号 if(endpoint) /该点是端点，返回点在骨架中的位置 return i;else /该点不是结束点，看该点是否和某个端点很近 if（该点与所在段的2个端点的距离都很大（大于5） return i;if(crtSegmentLength = 10) /当前段的长度小，在这里停下return i; 函数

33、seg = GetSkelSegment(Skel, i, &endpoint) 给出一个骨架点，返回它的段号，并且返回它是否是端点的信息 详细：遍历每一段，如果i= =LEFT或i= =RIGHT,返回段号，(*endpoint) = true;如果FIRSTiLAST, 返回段号，(*endpoint) = false;LEFTRIGHTFIRSTLAST第一段函数cp = CloseToCP( &Startpos, originCP, CritPts, numCritPts, SP_SMALL_DISTANCE);SP_SMALL_DISTANCE=0.5/是否和关键点距离近是，则则返回

34、CritPts数组中的index否，返回-1遍历CritPts的点，计算距离|x1-x|+|y1-y|+|z1-z| SP_SMALL_DISTANCE函数 rk2( Startpos.x, Startpos.y, Startpos.z, sX, sY, sZ, STEP_SIZE,ForceField, &Nextpos); /Nextpos的坐标是由势场决定的OutForce=interpolation(x,y,z,sizx,sizy,sizz,Force_ini); /经过插值后的势场值OutForce归一化 Nextpos= Startpos+OutForce* STEP_SIZE,得

35、到第一层骨架后，先把骨架的段保存到二维数组numSegmentsLevel1中去for(i=0; i numSegmentsLevel1; i+) for(j=0; j Segmentsij;然后在计算低散度值点函数GetHighDivergencePoints(force, L, M, N, f, percHDPts, &HDPts, &numHDPts); 得到低散度点的坐标 怀疑函数名字写错了函数的详细步骤：定义数组(*HDPts) = new VoxelPositionDoubleMAX_NUM_HDPTS 初始化5000个点遍历立体点(1L-1,1M-1,1N-1),遇到 EXTER

36、IOR, BOUNDARY or SURF，则跳过计算散度 某点上下左右前后六点插值v0 = interpolation(x + vdist, y, z, L, M, N, ForceField);v1 = interpolation(x - vdist, y, z, L, M, N, ForceField);v2 = interpolation(x, y + vdist, z, L, M, N, ForceField);v3 = interpolation(x, y - vdist, z, L, M, N, ForceField);v4 = interpolation(x, y, z + vdist, L, M, N, ForceField);v5 = interpolation(x, y, z - vdist, L, M, N, ForceField);div = (v0.xd - v1.xd) + (v2.yd - v3.yd) + (v4.zd - v5.zd) / (2 * vdist);寻找到散度值最小值minDiv和最大值maxDi