种子点生长算法下——三维种子点生长

上一篇文章“二维种子点生长研究”介绍了二维种子点生长，对于平面图像，二维种子点生长其实很容易理解，其典型应用是“油漆桶”。例如在一个具有多片联通区域的二维图上，我们能用油漆桶寻找所需要的单片联通区域。

　　至于三维种子点生长，其目的与二维种子点生长一样，都是要找出图像中的联通区域。三维图像由于其多了一维，所以三维的联通区域就不太好用二维图片来查看了。关于三位图片的介绍可以参考第一篇“图像数据的组织方式”。不过借助一些三维可视化软件，如ParaView，我们也可以观察三维图像的内容。

　　下面的图片是显示在ParaView下打开一张三维图片Lobster.raw的结果,我们通过调节不同体素值的透明度，就能凸显出三维图像中所需要关心的内容。

图像预览
像素透明度曲线

　　从上图看，由于三维图像本身是个长方体，各个像素有自己的颜色，所以只能看到外表面的样子。顺便提一下，在三维图像中，我们一般不使用“像素”这个词，而是使用“体素”来形容组成图像的这些单元。我们可以使用ParaView体绘制属性中调节透明度的选项，把一部分值较小的体素设为透明，就能看到这个三维图像中，值较大的体素就像一个个小砖头一样堆成了一个龙虾的形状。也就是说，如果我们在这些体素中找一个种子点，利用三维图像种子点生长算法，就能找到所有组成这个龙虾的体素。

　　上一篇文章已经介绍了二维种子点生长算法的相关知识，介绍了三种种子点生长算法—泛洪法、扫描线法和区段法，并在最后分析了他们的性能。本文的核心部分是三维种子点点生长的讨论，也就是把上一篇文章的算法扩展到三维，那么本文也按照相同的结构来说明这三种算法是如何扩展到三维的。

泛洪法
扫描线法
区段法
算法分析与实验

一、泛洪法

　　泛洪法的基本思想上篇文章详细介绍过，在三维图像中，需要扩展的主要是邻域的范围。邻域从二维扩展到三维后，就不再是四邻域和八邻域而是六邻域和二十六邻域。示意图如下：

　　其中绿色的点表示为当前点（蓝色点）的6邻域，绿色+红色的点为这个像素的26邻域。

　　对于三维点P，P的6邻域可以表示为：

　　　　S=Adj6(P)=｛P0(P.X-1,P.Y,P.Z)，P1(P.X+1,P.Y,P.Z)，P2(P.X,P.Y-1,P.Z)，P3(P.X,P.Y+1P.Z)，P4(P.X,P.Y,P.Z-1)，P5(P.X,P.Y,P.Z+1)｝。

　　P的26邻域可以表示为：

　　　　S=Adj26(P)=｛P0(P.X-1,P.Y,P.Z)，P1(P.X+1,P.Y,P.Z)，P2(P.X,P.Y-1,P.Z)，P3(P.X,P.Y+1,P.Z)，

　　　　　　　　　　　　P4(P.X-1,P.Y-1,P.Z)，P5(P.X+1,P.Y+1,P.Z)，P6(P.X+1,P.Y-1,P.Z)，P7(P.X-1,P.Y+1,P.Z)，

　　　　　　　　　　　　P8(P.X-1,P.Y-1,P.Z-1)，P9(P.X+1,P.Y+1,P.Z-1)，P10(P.X+1,P.Y-1,P.Z-1)，P11(P.X-1,P.Y+1,P.Z-1)

　　　　　　　　　　　　P12(P.X,P.Y-1,P.Z-1)，P13(P.X,P.Y+1,P.Z-1)，P14(P.X+1,P.Y,P.Z-1)，P15(P.X-1,P.Y,P.Z-1)，

　　　　　　　　　　　　P16(P.X,P.Y,P.Z-1)，P17(P.X,P.Y,P.Z+1)，P18(P.X+1,P.Y-1,P.Z+1)，P19(P.X-1,P.Y+1,P.Z+1)，

　　　　　　　　　　　　P20(P.X-1,P.Y-1,P.Z+1)，P21(P.X+1,P.Y+1,P.Z+1)，P22(P.X,P.Y-1,P.Z+1)，P23(P.X,P.Y+1,P.Z+1)

　　　　　　　　　　　　P24(P.X-1,P.Y,P.Z+1)，P25(P.X+1,P.Y,P.Z+1)｝。

　　那么实际上，三维泛洪法采用和二维泛洪法一样的逻辑。

FloodFill(seed,bitmap,includePredicate,process)set all postions of flagMap to falseset container Q to empty.push seed into Qset seed position in flagMap to trueprocess(seed)while Q is not emptypop a point P from Qforeach point T in Adj(P)if  includePredicate(T) is trueset position of T in flagMap to truepush T into Qprocess(T)end

　　注意在三维的情况下，相应的点类、图像类、位图标记表类等基础数据结构都需要分别增加一维。相应类的c#定义如下：

    public struct Int16Triple{public int X;public int Y;public int Z;public Int16Triple(int x, int y,int z){X = x;Y = y;Z = z;}}public class BitMap3d{public const byte WHITE = 255;public const byte BLACK = 0;public byte[] data;public int width;public int height;public int depth;public BitMap3d(int width, int height,int depth, byte v){this.width = width;this.height = height;this.depth = depth;data = new byte[width * height * depth];for (int i = 0; i < width * height*depth; i++)data[i] = v;}public BitMap3d(byte[] data, int width, int height,int depth){this.data = data;this.width = width;this.height = height;this.depth = depth;}public void SetPixel(int x, int y,int z,byte v){data[x + y * width+z*width*height] = v;}public byte GetPixel(int x, int y,int z){return data[x + y * width+z*width*height];}public void ReadRaw(string path){FileStream fs = new FileStream(path, FileMode.Open, FileAccess.Read);fs.Read(data, 0, width*height*depth);fs.Close();}public void SaveRaw(string path){FileStream fs = new FileStream(path, FileMode.OpenOrCreate, FileAccess.Write);fs.Write(data, 0, data.Length);fs.Close();}}public class FlagMap3d{public int width;public int height;public int depth;BitArray flags;public FlagMap3d(int width, int height,int depth){this.width = width;this.height = height;this.depth = depth;flags = new BitArray(width * height*depth, false);}public void SetFlagOn(int x, int y, int z,bool v){flags[x + y * width+z*width*height] = v;}public bool GetFlagOn(int x, int y,int z){return flags[x + y * width+z*width*height];}
}

　　相应实现的c#代码如下：

class FloodFill3d{protected BitMap3d bmp;protected FlagMap3d flagsMap;protected Container<Int16Triple> container;protected int count = 0;public byte min;public byte max;protected bool IncludePredicate(Int16Triple p){byte v = bmp.GetPixel(p.X, p.Y, p.Z);return v > min && v < max;}protected void Process(Int16Triple p){count++;return;}protected virtual void ExcuteFloodFill(BitMap3d data, Int16Triple seed){this.bmp = data;data.ResetVisitCount();flagsMap = new FlagMap3d(data.width, data.height,data.depth);Int16Triple[] adjPoints6 = new Int16Triple[6];flagsMap.SetFlagOn(seed.X, seed.Y, seed.Z,true);container.Push(seed);Process(seed);while (!container.Empty()){Int16Triple p = container.Pop();InitAdj6(ref adjPoints6, ref p);for (int adjIndex = 0; adjIndex < 6; adjIndex++){Int16Triple t = adjPoints6[adjIndex];if (t.X < data.width && t.X >= 0 && t.Y < data.height && t.Y >= 0&&t.Z<data.depth&&t.Z>=0){if (!flagsMap.GetFlagOn(t.X, t.Y,t.Z) && IncludePredicate(t)){flagsMap.SetFlagOn(t.X, t.Y,t.Z, true);container.Push(t);Process(t);}}}}return;}protected void InitAdj6(ref Int16Triple[] adjPoints6, ref Int16Triple p){adjPoints6[0].X = p.X - 1;adjPoints6[0].Y = p.Y;adjPoints6[0].Z = p.Z;adjPoints6[1].X = p.X + 1;adjPoints6[1].Y = p.Y;adjPoints6[1].Z = p.Z;adjPoints6[2].X = p.X;adjPoints6[2].Y = p.Y - 1;adjPoints6[2].Z = p.Z;adjPoints6[3].X = p.X;adjPoints6[3].Y = p.Y + 1;adjPoints6[3].Z = p.Z;adjPoints6[4].X = p.X;adjPoints6[4].Y = p.Y;adjPoints6[4].Z = p.Z-1;adjPoints6[5].X = p.X;adjPoints6[5].Y = p.Y;adjPoints6[5].Z = p.Z+1;}}

二、扫描线法

　　二维的扫描线法，是从一个点弹出堆栈后，对其左右扫描出两个端点，然后再对端点范围的两侧即Y-1和Y+1方向进行CheckRange。算法拓展到三维，只需增加检查Z-1和Z+1方向即可。也就是“两侧”变成了“四侧”，即:

（xleft，p.Y-1，p.Z）~（xright，p.Y-1，p.Z）
（xleft，p.Y+1，p.Z）~（xright，p.Y+1，p.Z）
（xleft，p.Y，p.Z-1）~（xright，p.Y，p.Z-1）
（xleft，p.Y，p.Z+1）~（xright，p.Y，p.Z+1）

　　扫描这四条线等价于6向泛洪法，假如扫描下面八条线，就相当于26向泛洪法

（xleft-1，p.Y-1，p.Z）~（xright+1，p.Y-1，p.Z）
（xleft-1，p.Y+1，p.Z）~（xright+1，p.Y+1，p.Z）
（xleft-1，p.Y，p.Z-1）~（xright+1，p.Y，p.Z-1）
（xleft-1，p.Y，p.Z+1）~（xright+1，p.Y，p.Z+1）
（xleft-1，p.Y-1，p.Z-1）~（xright+1，p.Y-1，p.Z-1）
（xleft-1，p.Y+1，p.Z+1）~（xright+1，p.Y+1，p.Z+1）
（xleft-1，p.Y+1，p.Z-1）~（xright+1，p.Y+1，p.Z-1）
（xleft-1，p.Y-1，p.Z+1）~（xright+1，p.Y-1，p.Z+1）

　　扫描线法拓展到三维，基本操作FindXleft，FindXright，CheckRange分别修改至相应的三维即可。其作用没有任何变化。

　　下面是三维扫描线线的c#代码：

class ScanlineFill3d{protected int count = 0;protected Container<Int16Triple> container;//这个容器可以是Queue和Stack中任意一种，这里抽象成一个Containerprotected BitMap3d bmp;public FlagMap3d flagsMap;protected virtual void ExcuteScanlineFill(BitMap3d data, Int16Triple seed){this.bmp = data;data.ResetVisitCount();flagsMap = new FlagMap3d(data.width, data.height, data.depth);flagsMap.SetFlagOn(seed.X, seed.Y, seed.Z, true);container.Push(seed);Process(seed);while (!container.Empty()){Int16Triple p = container.Pop();int xleft = FindXLeft(p);int xright = FindXRight(p);if (p.Y - 1 >= 0)CheckRange(xleft, xright, p.Y - 1, p.Z);if (p.Y + 1 < data.height)CheckRange(xleft, xright, p.Y + 1, p.Z);if (p.Z - 1 >= 0)CheckRange(xleft, xright, p.Y, p.Z - 1);if (p.Z + 1 < data.depth)CheckRange(xleft, xright, p.Y, p.Z + 1);}}//该函数为扫描线法主体protected void CheckRange(int xleft, int xright, int y, int z){for (int i = xleft; i <= xright; ){if ((!flagsMap.GetFlagOn(i, y, z)) && IncludePredicate(i, y, z)){int rb = i + 1;while (rb <= xright && (!flagsMap.GetFlagOn(rb, y, z)) && IncludePredicate(rb, y, z)){rb++;}rb--;Int16Triple t = new Int16Triple(rb, y, z);flagsMap.SetFlagOn(rb, y, z, true);container.Push(t);Process(t);i = rb + 1;}else{i++;}}}//CheckRange操作protected int FindXLeft(Int16Triple p){int xleft = p.X - 1;while (true){if (xleft == -1 || flagsMap.GetFlagOn(xleft, p.Y, p.Z)){break;}else{byte value = bmp.GetPixel(xleft, p.Y, p.Z);if (IncludePredicate(xleft, p.Y, p.Z)){Int16Triple t = new Int16Triple(xleft, p.Y, p.Z);flagsMap.SetFlagOn(xleft, p.Y, p.Z, true);Process(t);xleft--;}else{break;}}}return xleft + 1;}//FindXLeft操作protected int FindXRight(Int16Triple p){int xright = p.X + 1;while (true){if (xright == bmp.width || flagsMap.GetFlagOn(xright, p.Y, p.Z)){break;}else{byte value = bmp.GetPixel(xright, p.Y, p.Z);if (IncludePredicate(xright, p.Y, p.Z)){Int16Triple t = new Int16Triple(xright, p.Y, p.Z);flagsMap.SetFlagOn(xright, p.Y, p.Z, true);Process(t);xright++;}else{break;}}}return xright - 1;}//FindXRight操作public byte min;public byte max;protected bool IncludePredicate(int x,int y,int z){byte v = bmp.GetPixel(x, y,z);return v > min && v < max;}protected void Process(Int16Triple p){count++;}}

三、区段法

　　区段法扩展到三维也遵循和扫描线法一样的原则，相应的基本操作FindXleft，FindXright，CheckRange需要分别修改到三维的情况。同时在主函数逻辑中，监测相应的相邻区段时要注意ParentDirection所指示的父区段的特殊性。

　　下面是三维区段法的C#代码：

enum ParentDirections{Y0 = 1, Y2 = 2,Z0=3,Z2=4, Non = 5}enum ExtendTypes{LeftRequired = 1, RightRequired = 2, AllRez = 3, UnRez = 4}struct Span{public int XLeft;public int XRight;public int Y;public int Z;public ExtendTypes Extended;public ParentDirections ParentDirection;}class SpanFill3d{protected int count = 0;protected BitMap3d bmp;public FlagMap3d flagsMap;protected Container<Span> container;//以Span为单位的Queue或Stack容器protected virtual void ExcuteSpanFill(BitMap3d data, Int16Triple seed){this.bmp = data;data.ResetVisitCount();flagsMap = new FlagMap3d(data.width, data.height,data.depth);Process(seed);flagsMap.SetFlagOn(seed.X, seed.Y,seed.Z, true);Span seedspan = new Span();seedspan.XLeft = seed.X;seedspan.XRight = seed.X;seedspan.Y = seed.Y;seedspan.Z = seed.Z;seedspan.ParentDirection = ParentDirections.Non;seedspan.Extended = ExtendTypes.UnRez;container.Push(seedspan);while (!container.Empty()){Span span = container.Pop();#region AllRezif (span.Extended == ExtendTypes.AllRez){if (span.ParentDirection == ParentDirections.Y2){if (span.Y - 1 >= 0)//[spx-spy,y-1,z]CheckRange(span.XLeft, span.XRight, span.Y - 1, span.Z, ParentDirections.Y2);if (span.Z - 1 >= 0)//[spx-spy,y,z-1]CheckRange(span.XLeft, span.XRight, span.Y, span.Z - 1, ParentDirections.Z2);if (span.Z + 1 < data.depth)//[spx-spy,y,z+1]CheckRange(span.XLeft, span.XRight, span.Y, span.Z + 1, ParentDirections.Z0);continue;}if (span.ParentDirection == ParentDirections.Y0){if (span.Y + 1 < bmp.height)//[spx-spy,y+1,z]CheckRange(span.XLeft, span.XRight, span.Y + 1, span.Z, ParentDirections.Y0);if (span.Z - 1 > 0)//[spx-spy,y,z-1]CheckRange(span.XLeft, span.XRight, span.Y, span.Z - 1, ParentDirections.Z2);if (span.Z + 1 < data.depth)//[spx-spy,y,z+1]CheckRange(span.XLeft, span.XRight, span.Y, span.Z + 1, ParentDirections.Z0);continue;}if (span.ParentDirection == ParentDirections.Z2){if (span.Y - 1 >= 0)//[spx-spy,y-1,z]CheckRange(span.XLeft, span.XRight, span.Y - 1, span.Z, ParentDirections.Y2);if (span.Y + 1 < bmp.height)//[spx-spy,y+1,z]CheckRange(span.XLeft, span.XRight, span.Y + 1, span.Z, ParentDirections.Y0);if (span.Z - 1 >= 0)//[spx-spy,y,z-1]CheckRange(span.XLeft, span.XRight, span.Y, span.Z - 1, ParentDirections.Z2);continue;}if (span.ParentDirection == ParentDirections.Z0){if (span.Y - 1 >= 0)CheckRange(span.XLeft, span.XRight, span.Y - 1, span.Z, ParentDirections.Y2);if (span.Y + 1 < bmp.height)CheckRange(span.XLeft, span.XRight, span.Y + 1, span.Z, ParentDirections.Y0);if (span.Z + 1 < data.depth)CheckRange(span.XLeft, span.XRight, span.Y, span.Z + 1, ParentDirections.Z0);continue;}throw new Exception();}#endregion#region UnRezif (span.Extended == ExtendTypes.UnRez){int xl = FindXLeft(span.XLeft, span.Y,span.Z);int xr = FindXRight(span.XRight, span.Y,span.Z);if (span.ParentDirection == ParentDirections.Y2){if (span.Y - 1 >= 0)CheckRange(xl, xr, span.Y - 1, span.Z, ParentDirections.Y2);if (span.Y + 1 < bmp.height){if (xl != span.XLeft)CheckRange(xl, span.XLeft, span.Y + 1, span.Z, ParentDirections.Y0);if (span.XRight != xr)CheckRange(span.XRight, xr, span.Y + 1, span.Z, ParentDirections.Y0);}if (span.Z - 1 >= 0)CheckRange(xl, xr, span.Y, span.Z - 1, ParentDirections.Z2);if (span.Z + 1 < bmp.depth)CheckRange(xl, xr, span.Y, span.Z + 1, ParentDirections.Z0);continue;}if (span.ParentDirection == ParentDirections.Y0){if (span.Y + 1 < bmp.height)CheckRange(xl, xr, span.Y + 1, span.Z, ParentDirections.Y0);if (span.Y - 1 >= 0){if (xl != span.XLeft)CheckRange(xl, span.XLeft, span.Y - 1, span.Z, ParentDirections.Y2);if (span.XRight != xr)CheckRange(span.XRight, xr, span.Y - 1, span.Z, ParentDirections.Y2);}if (span.Z - 1 >= 0)CheckRange(xl, xr, span.Y, span.Z - 1, ParentDirections.Z2);if (span.Z + 1 < bmp.depth)CheckRange(xl, xr, span.Y, span.Z + 1, ParentDirections.Z0);continue;}if (span.ParentDirection == ParentDirections.Z2){if (span.Y - 1 >= 0)CheckRange(xl, xr, span.Y - 1, span.Z, ParentDirections.Y2);if (span.Y + 1 < bmp.height)CheckRange(xl, xr, span.Y + 1, span.Z, ParentDirections.Y0);if (span.Z - 1 >= 0)CheckRange(xl, xr, span.Y, span.Z - 1, ParentDirections.Z2);if (span.Z + 1 < bmp.depth){if (xl != span.XLeft)CheckRange(xl, span.XLeft, span.Y, span.Z+1, ParentDirections.Z0);if (span.XRight != xr)CheckRange(span.XRight, xr, span.Y, span.Z+1, ParentDirections.Z0);}continue;}if (span.ParentDirection == ParentDirections.Z0){if (span.Y - 1 >= 0)CheckRange(xl, xr, span.Y - 1, span.Z, ParentDirections.Y2);if (span.Y + 1 < bmp.height)CheckRange(xl, xr, span.Y + 1, span.Z, ParentDirections.Y0);if (span.Z - 1 >= 0){if (xl != span.XLeft)CheckRange(xl, span.XLeft, span.Y, span.Z -1, ParentDirections.Z2);if (span.XRight != xr)CheckRange(span.XRight, xr, span.Y, span.Z - 1, ParentDirections.Z2);}if (span.Z + 1 < bmp.depth)CheckRange(xl, xr, span.Y, span.Z + 1, ParentDirections.Z0);continue;}if (span.ParentDirection == ParentDirections.Non){if (span.Y + 1 < bmp.height)CheckRange(xl, xr, span.Y + 1, span.Z, ParentDirections.Y0);if (span.Y - 1 >= 0)CheckRange(xl, xr, span.Y - 1, span.Z, ParentDirections.Y2);if (span.Z - 1 >= 0)CheckRange(xl, xr, span.Y, span.Z - 1, ParentDirections.Z2);if (span.Z + 1 < bmp.depth)CheckRange(xl, xr, span.Y, span.Z + 1, ParentDirections.Z0);continue;}throw new Exception();}#endregion#region LeftRequiredif (span.Extended == ExtendTypes.LeftRequired){int xl = FindXLeft(span.XLeft, span.Y,span.Z);if (span.ParentDirection == ParentDirections.Y2){if (span.Y - 1 >= 0)CheckRange(xl, span.XRight, span.Y - 1, span.Z, ParentDirections.Y2);if (span.Y + 1 < bmp.height && xl != span.XLeft)CheckRange(xl, span.XLeft, span.Y + 1, span.Z, ParentDirections.Y0);if (span.Z - 1 >= 0)CheckRange(xl, span.XRight, span.Y , span.Z-1, ParentDirections.Z2);if (span.Z + 1 < bmp.depth)CheckRange(xl, span.XRight, span.Y , span.Z+1, ParentDirections.Z0);continue;}if (span.ParentDirection == ParentDirections.Y0){if (span.Y - 1 >= 0 && xl != span.XLeft)CheckRange(xl, span.XLeft, span.Y - 1, span.Z, ParentDirections.Y2);if (span.Y + 1 < bmp.height)CheckRange(xl, span.XRight, span.Y + 1, span.Z, ParentDirections.Y0);if (span.Z - 1 >= 0)CheckRange(xl, span.XRight, span.Y, span.Z - 1, ParentDirections.Z2);if (span.Z + 1 < bmp.depth)CheckRange(xl, span.XRight, span.Y, span.Z + 1, ParentDirections.Z0);continue;}if (span.ParentDirection == ParentDirections.Z2){if (span.Y - 1 >= 0 )CheckRange(xl, span.XRight, span.Y - 1, span.Z, ParentDirections.Y2); if (span.Y + 1 < bmp.height)CheckRange(xl, span.XRight, span.Y + 1, span.Z, ParentDirections.Y0);if (span.Z - 1 >= 0)CheckRange(xl, span.XRight, span.Y, span.Z - 1, ParentDirections.Z2);if (span.Z + 1 < bmp.depth && xl != span.XLeft)CheckRange(xl, span.XLeft, span.Y, span.Z + 1, ParentDirections.Z0);continue;}if (span.ParentDirection == ParentDirections.Z0){if (span.Y - 1 >= 0)CheckRange(xl, span.XRight, span.Y - 1, span.Z, ParentDirections.Y2);if (span.Y + 1 < bmp.height)CheckRange(xl, span.XRight, span.Y + 1, span.Z, ParentDirections.Y0);if (span.Z - 1 >= 0 && xl != span.XLeft)CheckRange(xl, span.XLeft, span.Y, span.Z - 1, ParentDirections.Z2);if (span.Z + 1 < bmp.depth )CheckRange(xl, span.XRight, span.Y, span.Z + 1, ParentDirections.Z0);continue;}throw new Exception();}#endregion#region RightRequiredif (span.Extended == ExtendTypes.RightRequired){int xr = FindXRight(span.XRight, span.Y,span.Z);if (span.ParentDirection == ParentDirections.Y2){if (span.Y - 1 >= 0)CheckRange(span.XLeft, xr, span.Y - 1,span.Z, ParentDirections.Y2);if (span.Y + 1 < bmp.height && span.XRight != xr)CheckRange(span.XRight, xr, span.Y + 1,span.Z, ParentDirections.Y0);if (span.Z - 1 >= 0)CheckRange(span.XLeft, xr, span.Y, span.Z-1, ParentDirections.Z2);if (span.Z + 1 < bmp.depth)CheckRange(span.XLeft, xr, span.Y, span.Z+1, ParentDirections.Z0); continue;}if (span.ParentDirection == ParentDirections.Y0){if (span.Y + 1 < bmp.height)CheckRange(span.XLeft, xr, span.Y + 1, span.Z, ParentDirections.Y0);if (span.Y - 1 >= 0 && span.XRight != xr)CheckRange(span.XRight, xr, span.Y - 1, span.Z, ParentDirections.Y2);if (span.Z - 1 >= 0)CheckRange(span.XLeft, xr, span.Y, span.Z - 1, ParentDirections.Z2);if (span.Z + 1 < bmp.depth)CheckRange(span.XLeft, xr, span.Y, span.Z + 1, ParentDirections.Z0); continue;}if (span.ParentDirection == ParentDirections.Z2){if (span.Y - 1 >= 0)CheckRange(span.XLeft, xr, span.Y - 1, span.Z, ParentDirections.Y2);if (span.Y + 1 < bmp.height)CheckRange(span.XLeft, xr, span.Y + 1, span.Z, ParentDirections.Y0);if (span.Z - 1 >= 0)CheckRange(span.XLeft, xr, span.Y, span.Z - 1, ParentDirections.Z2);if (span.Z + 1 < bmp.depth && span.XRight != xr)CheckRange(span.XRight, xr, span.Y, span.Z + 1, ParentDirections.Z0);continue;}if (span.ParentDirection == ParentDirections.Z0){if (span.Y - 1 >= 0)CheckRange(span.XLeft, xr, span.Y - 1, span.Z, ParentDirections.Y2);if (span.Y + 1 < bmp.height)CheckRange(span.XLeft, xr, span.Y + 1, span.Z, ParentDirections.Y0);if (span.Z - 1 >= 0 && span.XRight != xr)CheckRange(span.XRight, xr, span.Y, span.Z - 1, ParentDirections.Z2);if (span.Z + 1 < bmp.depth)CheckRange(span.XLeft, xr, span.Y, span.Z + 1, ParentDirections.Z0); continue;}throw new Exception();}#endregion}}protected void CheckRange(int xleft, int xright, int y, int z,ParentDirections ptype){for (int i = xleft; i <= xright; ){if ((!flagsMap.GetFlagOn(i, y,z)) && IncludePredicate(i, y,z)){int lb = i;int rb = i + 1;while (rb <= xright && (!flagsMap.GetFlagOn(rb, y,z)) && IncludePredicate(rb, y,z)){rb++;}rb--;Span span = new Span();span.XLeft = lb;span.XRight = rb;span.Y = y;span.Z = z;if (lb == xleft && rb == xright){span.Extended = ExtendTypes.UnRez;}else if (rb == xright){span.Extended = ExtendTypes.RightRequired;}else if (lb == xleft){span.Extended = ExtendTypes.LeftRequired;}else{span.Extended = ExtendTypes.AllRez;}span.ParentDirection = ptype;for (int j = lb; j <= rb; j++){flagsMap.SetFlagOn(j, y,z, true);Process(new Int16Triple(j, y,z));}container.Push(span);i = rb + 1;}else{i++;}}}//区段法的CheckRange 注意与扫描线的CheckRange的不同protected int FindXRight(int x, int y,int z){int xright = x + 1;while (true){if (xright == bmp.width || flagsMap.GetFlagOn(xright, y,z)){break;}else{if (IncludePredicate(xright, y,z)){Int16Triple t = new Int16Triple(xright, y,z);flagsMap.SetFlagOn(xright, y,z ,true);Process(t);xright++;}else{break;}}}return xright - 1;}protected int FindXLeft(int x, int y,int z){int xleft = x - 1;while (true){if (xleft == -1 || flagsMap.GetFlagOn(xleft, y, z)){break;}else{if (IncludePredicate(xleft, y,z)){Int16Triple t = new Int16Triple(xleft, y,z);flagsMap.SetFlagOn(xleft, y,z, true);Process(t);xleft--;}else{break;}}}return xleft + 1;}public byte min;public byte max;protected bool IncludePredicate(int x, int y, int z){byte v = bmp.GetPixel(x, y, z);return v > min && v < max;}protected void Process(Int16Triple p){count++;}}

四、算法测试与实验

　　算法测试采用来自http://www.volvis.org的5组体数据，注意这里的阈值范围是值一个体素值的范围（min,max），在此范围内的体素才纳入区域。也就是说本文的includePredicate不再是跟上一篇一样是判断像素灰度为255的纳入区域，而是判断是否在这个（min,max）的范围内。

图像预览	数据名称	参数说明
	Lobster	大小：301×324×56 文件：Lobster.raw 种子点：(124, 168, 27) 阈值范围：37~255 描述：CT scan of a lobster contained in a block of resin.
	engine	大小：256×256×128 文件：Engine.raw 种子点：(149, 44, 43) 阈值范围：64~255 描述：CT scan of two cylinders of an engine block.
	backpack	大小：301×324×56 文件：Backpack.raw 种子点：(53, 390, 160) 阈值范围：46~255 描述：CT scan of a backpack filled with items.
	phantom	大小：512×512×442 文件：Phantom.raw 种子点：（256，256，227）阈值范围：42~255 描述：CT scan of a Colon phantom with several different objects and five pedunculated large polyps in the central object.
	cube	大小：200×200×200 文件：Cube.raw 种子点：（50，50，50）阈值范围：0~255 描述：全白色，用于测试全部充满的极端情况

　　测试结果如下：

测试数据	泛洪法（栈式）	泛洪法(队列式)	扫描线法（栈式）	扫描线法（队列式）	区段法（栈式）	区段法（队列式）	区域点数
lobster	86ms	76ms	56ms	64ms	40ms	47ms	277367
engine	329ms	345ms	216ms	254ms	137ms	183ms	1154807
backpack	596ms	638ms	426ms	470ms	319ms	351ms	19115450
phantom	13468ms	14118ms	6520ms	8279ms	3701ms	4066ms	39759257
cube	2120ms	2259ms	1292ms	1404ms	691ms	692ms	800000

　　对于其中的backpack数据的单元访问比例统计如下：

算法	GetPixel/总点数	GetFlagOn/总点数	SetFlagOn/总点数	Push和Pop/总点数	总点数	时间花费
泛洪法（栈式）	1.945581	5.997103	1.0	1.0	1915450	596ms
泛洪法（队列式）	1.945581	5.997103	1.0	1.0	1915450	638ms
扫描线法（栈式）	4.165088	5.337406	1.0	0.2668	1915450	426ms
扫描线法（队列式）	2.715424	5.774523	1.0	0.7589	1915450	470ms
区段法（栈式）	1.911485	3.526728	1.0	0.2147	1915450	319ms
区段法（队列式）	1.931963	3.754333	1.0	0.3628	1915450	351ms

　　通过测试可以看出，相比于二维种子点生长，由于多了一维，所以利用X轴数据连续特性所能达到增加效率的效果明显减弱。既可以从时间结果上看，也可以从单元访问比例上看，都能得到这一结论。不过从整体上看，效率上仍然能够得出如下结论：栈式算法略块于队列式算法；区段法快于扫描线法，扫描线法快于泛洪法。

种子点生长算法下——三维种子点生长相关推荐

生长算法实现点集的三角剖分（Python（Tkinter模块））
生长算法实现点集的三角剖分( Python(Tkinter模块)) 关于三角剖分假设V是二维实数域上的有限点集,边e是由点集中的点作为端点构成的封闭线段, E为e的集合.那么该点集V的一个三角剖分T ...
种子点生长算法（上）——二维种子点生长
本文只用作笔记,方便日后查阅. 种子点生长算法上--二维种子点生长下文提到的种子点生长算法,包括泛洪法,扫描线法,区段法三种.文本先从最简单的泛洪法入手介绍种子点生长算法的相关概念.之后进一步讨论了 ...
python 区域生长算法_多种子的区域生长算法
摘要:多种子的区域生长算法,基于C++编写. 关键字:图像处理, 种子生长, 区域生长 1. 题外话最近需要找一种简单对图像进行分割的算法,作为后续算法的前处理阶段.最开始想到的是聚类,但是聚类会有 ...
一颗种子，一颗小树苗在快速生长长大的过程中，遇到风雨在所难免
合抱之木,生于毫末,这个道理已经流传了千年,而且还将继续流传千万年下去. 所以,鉴往知来. 您要想有合抱之木的事业,必须先要埋下合抱之木的种子,然后才是然后的事情.现在的您,如果已经埋下了种子,那么, ...
多种子的区域生长算法
摘要:多种子的区域生长算法,基于C++编写. 关键字:图像处理, 种子生长, 区域生长本文是我旧博客中的博文,在CSDN图片显示不正常,请移步旧博客查看:https://imlogm.github. ...
matlab下三维dla模型模拟,Matlab下三维DLA模型模拟
Matlab下三维DLA模型模拟 2007-01-11 19:18 function dla3dv5(Nsum,Wstep) %定义dla函数,Nsum为所生成絮体包含的颗粒数,Wstep为计算过程中 ...
叶脉图案以及藤蔓生长算法在houdini里面的实现 Leaf Venation
继续跟着CGworkshops里面的Shawn Lipowski大神学习在houdini里面写Vex,昨天简单看了一下他做植物生长的那一课,了解了一下思路于是开始自己着手写,本人不太习惯一步一步跟着学 ...
C++实现Delaunay三角网生长算法
目录一.概述 1.1 三角网的介绍 1.2 Delaunay三角形二.三角网生长算法 2.1 建立第一个三角形 2.2 扩展TIN 三.各部分代码实现 3.1 数据结构 3.2 算法过程 3.3 ...
c++ 凸包分治算法_三维凸包
缘起众所周知,二维凸包可以使用 Graham 扫描内解决. 所以本文来学习一下三维空间中凸包的一种直观算法--增量算法(increment algorithm) 分析有一条叫 Willy 的苹果 ...
【APF三维路径规划】基于matlab人工势场算法无人机三维路径规划【含Matlab源码 168期】
一.获取代码方式获取代码方式1: 通过订阅紫极神光博客付费专栏,凭支付凭证,私信博主,可获得此代码. 获取代码方式2: 完整代码已上传我的资源:[三维路径规划]基于matlab人工势场算法无人机三维 ...

种子点生长算法下——三维种子点生长

种子点生长算法下——三维种子点生长相关推荐

最新文章

热门文章