三维骨架提取Main教材.docx

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

2、ouble z;struct CriticalPointVoxelPositionDoubleCriticalPointTypeVectordoubleposition;type;evect3;eval3;/关键点的坐标/关键点的类型/3 个特征向量数组/3 个特征值;struct SkeletonVoxelPositionDouble *Points; int sizePoints;/存储骨架上的点/初始化时分配的最大储存空间本程序是50000int numPoints;int* Segments;/骨架上的点的个数/骨架段为二维数组/numSegments行4 列的数组LEFT RIGHT

3、 FIRST LASTint sizeSegments; int numSegments;/初始化分配的最大存储空间，本程序是/骨架段的个数5000;骨架段的说明：LEFTRIGHTFIRSTLAST第一段第一个骨架点index最后一个骨架点的 index第二段第三段立体长宽高：L=123M=116N计算势场用到的距离的次方数fieldStrenght = 8;计算散度值时和阈值有关的percHDPts= 3;= 114;vfin = false;vfout = false;interactive = false;/这些变量代表的实际意义不清楚#define SURF100/ surface

4、voxel#define BOUNDARY110/ boundary voxel - participates in potential field/calculation#define INTERIOR200/ interior voxel#define PADDING_MIN210/ added voxels in order to thick the object#define NR_PAD_VALUES40/are inthisrange: PADDING_MINtoPADDING_MIN + NR_PAD_V ALUES#define EXTERIOR0/ background (e

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

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

7、调函数的地址，输出骨架地址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的空间，存储了立体的值问：立体的值代表什么含义确保立体在三个方向填充了足够的空的平面，所以在距离边界特定的位置放置点电荷CheckVolumePadding(&f, &L, &M,

8、 &N, distCharges)问：这个函数到底做什么用的这个函数首先是对f 中的值划分为两类f 中的值0INTERIOR200求出非 0 和 i=0,L-1 j=0,M-1 k=0,N-1点值为 EXTERIOR 的点的坐标范围EXTERIOR0外部点是指外面的点吗minxmaxxminymaxyminzmaxz满足以下条件则返回真问：这个是什么意思，为什么这么处理minx - distCharges= 0miny - distCharges = 0minz - distCharges = 0maxx + distCharges = Lmaxy + distCharges = Mmaxz

9、+ distCharges= N下面执行的函数是，填充立体作用是 Make sure the volume does not have holes in it 问：这个函数的作用又是什么 MakeSolidV olume(L, M, N, f, EXTERIOR, INTERIOR);遍历所有的点，将不是 EXTERIOR 的点，也就是全部的 INTERIOR 设置为 OBJECT_V AL=1for(i=0; i N; i+)FloodFillZPlane(i, L, M, N, vol);用 OUTSIDE_1=2 填充每个 Z 面/ 判断该点周围的4 个点，如果是0（ EXTERIOR)

10、的话，设置为OUTSIDE_1=2/填充函数，遍历 XOY 平面，设置从原点开始遍历，原点的值设置为bool 变量 anyChange 开始为 falseOUTSIDE_1=2次序为左右，上下，如果4 点中出现了 EXTERIOR ，将 EXTERIOR 点设置为 OUTSIDE_1=2 ，设置 anyChange 开始为 true，否则不变，如果 anychange 的值为 true，则从右下点开始，遍历次序为右左，下上，设置同上。接下来填充Y 平面，函数算法同Z ，对值是EXTERIOR=0和 OUTSIDE_1=2的值设置为OUTSIDE_2=3for(i=0; i M; i+)Floo

11、dFillYPlane(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);下边是设置立体的值，分成三类SURFINTERIOREXTERIORbool FlagV olume(unsigned char *vol, int L, int M, int N)1 遍历所有的点（L*M*N），将点的值为0 的设置为EXTERIOR=0，非0 的设置为INTER

12、IOR=2002 遍历所有的 INTERIOR 点（假设 B 点），如果该点的邻接有置为 SURF，这里的邻接点是 B 点的上下左右前后六个点EXTERIOR点，则将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);检查 distCharges 的值不能超过NR_PAD_V ALUES=40IfI

13、fdistCharges=0 distCharges0不做任何操作return做以下操作，在立体的/程序中设定的是0SURF 外面加上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 邻接点）是EXTERIOR把点 B 值设置为 PADDING_MIN同理对于下边的for 循环是依次在

14、立体的外边加上一层，同样是对的值。点（点 B）EXTERIOR点设置为新If distCharges nrLayers; i-)nrLayers 层，可能会造成立体的断开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 check/ Make sure

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

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

17、ESHOLD| BOUNDARY.y-A.y|= PF_THRESHOLD| BOUNDARY.z-A.z |= PF_THRESHOLD满足这些条件的 BOUNDARY 点才对 A 点的势场值起作用，其他的忽略问：这种筛选方法有问题吗计算公式是FPiCi PR8for (k = 0; k N; k+)/遍历所有的点根据 PF_THRESHOLD得到zStartIndexzEndIndex的值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=8 forceidx.yd =

19、 forceidx.yd + (v2 / r fieldStrenght);forceidx.zd = forceidx.zd + (v3/ r fieldStrenght);force 向量场归一化( xd, yd, zd)xd,yd,zdyd2xd2yd2zd2yd2zd 2xd2zd2xd2下面计算 SURF 和 BOUNDARY点处的电势值，通过点的26 邻接点的电势值的平均值来确定，但是 26 个点中的 SURFBOUNDARYEXTERIOR点忽略得到关键点GetCriticalPoints(force, L, M, N, f, &CritPts, &numCritPts);Cri

20、tPts 存储得到的关键点为关键点数组分配 2000 的空间(*CritPts)=new CriticalPointMAX_NUM_CRITPTS)MAX_NUM_CRITPTS=2000SURF，(1L-1,1M-1,1N-1) 检查每个点（假设B 点）的邻接点，如果邻接点中出现了1BOUNDARY ， EXTERIOR 的情况，跳过点B， continue问：为什么选取这几个邻接点inds0 = idx;/B 点inds1 = idx + 1;/ 以下为 B 的邻接点inds2 = idx + L;inds3 = idx + L + 1;inds4 = idx + slsz;inds5 =

21、 idx + slsz + 1;inds6 = idx + slsz + L;inds7 = idx + slsz + L + 1;2判断点 (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(ChangeInS

22、ign(ForceField, inds, 8) ，则该点有可能成为关键点，将这个 cell 分成 8 个 subcells（二阶魔方）for(kk=0; kk 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

23、, y + jj/2.00, z + kk/2.00, 0.50,sX, sY, sZ, ForceField, critPt)return true;bool FindCriticalPointInFloatCell(double x, double y, double z, double cellsize,int sX, int sY, int sZ, Vector* ForceField, VoxelPositionDouble *critPt) 这个函数判断一个cell 是否是关键点，这里的(x,y,z) 表示每个 cell 的最靠近原点的点，用插值的方法得到立体的每个顶点的Force

24、Field 值这里 cellsize 的值是 0.5cv0 （ x,y,z）cv1 (x+cellsize,y,z)cv2(x,y+cellsize,z)cv3(x+cellsize,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

25、(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 再执行 FindCriticalPointInFloat

26、Cell这个函数cell。得到关键点后，再确定关键点的类型vdist = 1.00 / 1048576在每个关键点的上下左右前后插入六个点/根据插入点的位置来得到插入到该点的向量的大小，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 TNT:Array2D TNT:Array2DJac(3, 3, 0.00); EigVals(3, 3, 0.0); EigVects(3, 3, 0.0);二维数组

27、放雅克比行列式三维数组放特征值存放特征向量下边计算出雅克比行列式的特征值和特征向量JAMA:Eigenvalueev(Jac); ev.getD(EigVals);ev.getV(EigVects);EigVals 存放的是特征值，举例0.1533440.0000000.0000000.0000001.6078370.0000000.0000000.0000000.668626EigVects 存放的是对应的特征向量0.6983430.1196000.0779060.3096230.9623100.6449990.6453300.7242250.836460根据得到的特征值和特征向量判断出关键

28、点的类型EigVals000所有特征值 0EigVals110该点为吸引点CPT_ATTRACTING_NODEEigVals220所有特征值 0该点为斥点CPT_REPELLING_NODEEigVals220特征值正负都有该点为鞍点CPT_SADDLE对特征向量进行归一化存储每个关键点的特征值和特征向量到(*CritPts)i中，得到关键点的坐标类型特征值特征向量后，下边是得到第一骨架GetLevel1Skeleton(force, L, M, N, CritPts, numCritPts, Skel);定义一些相关变量/该数组用来标记某个关键点是否已经是骨架的一部分了critPtInSk

29、eleton = new intnumCritPts，初始化为 -1-1 表示还不是从鞍点出发，正特征值对应的特征向量为方向，直到下一点是生成骨架的点或者是关键点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, numCr

30、itPts,critPtInSkeleton,Skel);函数具体实现步骤：遍历所有的鞍点， 每条骨架段都是从鞍点开始的，第一次循环沿着方向向量作为方向得到下一点的坐标，接下来的循环将使用当前点（x,y,z）的势场值作为方向（注：（ x,y,z）可能不是整形的坐标，采用插值来求得其势场向量值），循环结束的条件很多1.step 过多，认为是Most likely we are chasing our own tail.2.遇到了骨架 A 中的点 a，则将该点加入到当前段，并且作为段结束的端点，将A 段在点a 处分开， A 段变成了 2 段3.遇到了关键点，则直接加入当前段，并且作为段结束的端点。

31、4.计算出来的 nextpos 和 startpos 值是相等的，不做处理，结束5.段到了立体的边界，将当前段删除以下是相关的函数s = AddSkeletonPoint(Skel, &Startpos);函数参数为骨架，待添加关键点的坐标返回值为新添加的点的索引（0 开始）将点 Startpos 添加到 Skel 的 Points 数组中添加新的骨架段crtSegment = AddNewSkelSegment(Skel, s);函数参数s 为该点在Skel-Point 的索引， 返回值是crtSegment 为段 Skel-Segments 的索引（ 0 开始）自动添加新的段(Skel-S

32、egmentsSkel-numSegments = new int4初始化为点 s 的索引LEFTRIGHTFIRSTLASTcrtSegmentssss函数 sp=CloseToSkel(&Startpos,oCPSkelPos,Skel, nrAlreadyInSkel,SP_SMALL_DISTANCE,crtSegLength);SP_SMALL_DISTANCE=0.5判断该点是否和一个骨架中的点离得很近（也就是遇到了骨架点，可能骨架的当前段在该点结束），如果是的话 ,返回骨架点的位置(inSkel ） i ，特别的，假如骨架点不是端点，但是和端点离得近，这里不返回骨架点的位置，返回

33、-1wait ait to get to oneofthe segment end points不是的话，返回-1详细：遍历所有的骨架点，计算距离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）if(c

34、rtSegmentLength = 10) / 当前段的长度小，在这里停下return i;return i;函数 seg = GetSkelSegment(Skel, i, &endpoint)给出一个骨架点，返回它的段号，并且返回它是否是端点的信息详细：遍历每一段，如果i= =LEFT 或 i= =RIGHT, 返回段号， (*endpoint) = true;如果 FIRSTiLAST,返回段号，(*endpoint) = false;LEFTRIGHTFIRSTLAST第一段函数 cp = CloseToCP( &Startpos, originCP, CritPts, numCrit

35、Pts, SP_SMALL_DISTANCE); SP_SMALL_DISTANCE=0.5/是否和关键点距离近是，则则返回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

36、_ini); / 经过插值后的势场值 OutForce 归一化 Nextpos= Startpos+OutForce* STEP_SIZE,得到第一层骨架后，先把骨架的段保存到二维数组numSegmentsLevel1for(i=0; i numSegmentsLevel1; i+)for(j=0; j Segmentsij;中去然后在计算低散度值点函数 GetHighDivergencePoints(force, L, M, N, f, percHDPts, &HDPts, &numHDPts);得到低散度点的坐标怀疑函数名字写错了函数的详细步骤：定义数组 (*HDPts) = new V

37、oxelPositionDoubleMAX_NUM_HDPTS初始化 5000 个点遍历立体点 (1L-1,1M-1,1N-1), 遇到EXTERIOR, 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 和最大值maxDiv计算阈值