IJ中的FindMaxima算法
某个图像处理的任务需求被转化为一个寻找图像局部最大值的问题，例如下面的图像是一段来自医学投影的组织影像，其中用肉眼能够辨识出一些灰度值比较大的团块，这些团块都是代表着一定医学含义的组织。需求是找出类似这种团块的个数。
 |
 |
| 原图 | 需要找的团块 |
后来在ImageJ中发现了一个这个比较奇妙的算法，叫做FindMaxima，这个算法支持提供一个ErorrTolerance参数之后，根据这个参数来调整寻找到极大位置的个数。例如选择ErorrTolerance为下面三个不同的数值后，算出来的团块数是不同的。
 |
 |
 |
| IJ中的算法 | 算法参数面板 | 算法结果 |
作者的描述
根据ImageJ中提供的信息可以找到这个算法的原作者，叫做Michael Schmid，在论坛上可以找到他对自己写的这个算法的简单描述：
|
Hi, sorry, there is no publication on it (at least none that I am aware of, probably I have reinvented the wheel), but it is rather simple code: (1) Find the local maxima (2) Sort them in descending sequence (3) For each local maximum, do a flood fill algorithm with the gray level tolerance (without modifying the original, it is done on a temporary scratch image). Maxima where flood filling reaches a previously filled area (i.e., area of other maximum within the tolerance) are discarded. (4) In case of 'single points' output, if there are several points having the highest value inside the flood-filled area, use the one that is closest to their geometric center.
Segmentation is a bit more difficult (it stems from older ImageJ versions), and I am currently about to see whether it can't be done faster and with better accuracy... ------Michael |
总结这个描述，可以按照如下思路简单的实现这个算法：
- 寻找到所有比8邻域像素值都大的像素点，存入数组L
- 对L按像素值从高到低排序
- 对每一个L中的像素P
- 以P为种子点执行泛洪法，包含原则为与P的像素值相差在tolerance之内，若遇到比P像素值更高的点或者遇到了被标记为Maxima的点，则此P点不为Maxima，否则P被标记为Maxima。
- 将L中标记为Maxima的体素输出
根据描述实现一个轻量版
根据上述思路用C#实现一个版本，其代码如下：
public struct Int16Double { public int X; public int Y; public Int16Double(int x, int y) { X = x; Y = y; } } public struct Int16DoubleWithValue:IComparable<Int16DoubleWithValue> { public int X; public int Y; public float V; public Int16DoubleWithValue(int x, int y, float value) { X = x; Y = y; V = value; } public int CompareTo(Int16DoubleWithValue other) { return -this.V.CompareTo(other.V); } } public class BitMap2d { public float[] data; public int width; public int height; public BitMap2d(int width, int height, float v) { this.width = width; this.height = height; data = new float[width * height]; for (int i = 0; i < width * height; i++) data[i] = v; } public void SetPixel(int x, int y, byte v) { data[x + y * width] = v; } public float GetPixel(int x, int y) { return data[x + y * width]; } public void ReadRaw(string path) { FileStream fs = new FileStream(path, FileMode.Open, FileAccess.Read); BinaryReader sr = new BinaryReader(fs); for (int i = 0; i < width * height; i++) { byte[] floatBytes = sr.ReadBytes(4); // swap the bytes byte temp = floatBytes[0]; floatBytes[0] = floatBytes[3]; floatBytes[3] = temp; temp = floatBytes[1]; floatBytes[1] = floatBytes[2]; floatBytes[2] = temp; // get the float from the byte array float value = BitConverter.ToSingle(floatBytes, 0); data[i] = value; } sr.Close(); fs.Close(); return; } public Bitmap MakeBmp() { float min=float.MaxValue; float max=float.MinValue; for (int i = 0; i < width; i++) { for (int j = 0; j < height; j++) { float r=this.GetPixel(i,j); if(r>max) max=r; if(r<min) min=r; } } float delta=max-min; Bitmap bmp = new Bitmap(this.width, this.height, System.Drawing.Imaging.PixelFormat.Format32bppArgb); for (int i = 0; i < width; i++) { for (int j = 0; j < height; j++) { float r=this.GetPixel(i,j); int b=(int)(255*(r-min)/delta); Color c=Color.FromArgb((byte)b,(byte)b,(byte)b); bmp.SetPixel(i, j, c); } } return bmp; } } public class FlagMap2d { public int action_set_count; public int action_get_count; public int width; public int height; byte[] flags; public FlagMap2d(int width, int height,byte v) { this.width = width; this.height = height; action_get_count = 0; action_set_count = 0; flags = new byte[width * height]; for(int i=0;i<width*height;i++) flags[i] = v; } public void SetFlagOn(int x, int y, byte v) { flags[x + y * width] = v; action_set_count++; } public byte GetFlagOn(int x, int y) { action_get_count++; return flags[x + y * width]; } }
public class MaximunFinder { BitMap2d bmp; float torlerance; const byte UNPROCESSED = 0; const byte VISITED = 1; const byte PROCESSED = 2; static Int16Double[] Delta = new Int16Double[8] { new Int16Double(-1,-1), new Int16Double(-1,0), new Int16Double(-1,1), new Int16Double(0,-1), new Int16Double(0,1), new Int16Double(1,-1), new Int16Double(1,0), new Int16Double(1,1), }; public MaximunFinder(BitMap2d bmp,float torlerance) { this.bmp = bmp; this.torlerance = torlerance; } public List<Int16DoubleWithValue> FindMaxima() { List<Int16DoubleWithValue> list = FindLocalMaxima(); list.Sort(); FlagMap2d flag=new FlagMap2d(bmp.width,bmp.height,0); List<Int16DoubleWithValue> r=new List<Int16DoubleWithValue>(); List<Int16Double> temp=new List<Int16Double>(); for (int i = 0; i < list.Count; i++) { if (flag.GetFlagOn(list[i].X, list[i].Y) == UNPROCESSED) { bool ret = FloodFill(list[i].X, list[i].Y,temp,flag); if (ret) { r.Add(list[i]); MarkAll(temp, PROCESSED, flag); } else { MarkAll(temp, UNPROCESSED, flag); flag.SetFlagOn(list[i].X, list[i].Y, PROCESSED); } temp.Clear(); } } return r; } private List<Int16DoubleWithValue> FindLocalMaxima() { List<Int16DoubleWithValue> list = new List<Int16DoubleWithValue>(); for (int i = 1; i < bmp.width - 1; i++) { for (int j = 1; j < bmp.height - 1; j++) { if (IsMaxima(i, j)) { list.Add(new Int16DoubleWithValue(i, j,bmp.GetPixel(i,j))); } } } return list; } private bool IsMaxima(int i, int j) { float v = bmp.GetPixel(i, j); bool b1 = v > bmp.GetPixel(i - 1, j - 1); bool b2 = v > bmp.GetPixel(i, j - 1); bool b3 = v > bmp.GetPixel(i +1, j - 1); bool b4 = v > bmp.GetPixel(i - 1, j); bool b5 = v > bmp.GetPixel(i + 1, j); bool b6 = v > bmp.GetPixel(i - 1, j + 1); bool b7 = v > bmp.GetPixel(i, j + 1); bool b8 = v > bmp.GetPixel(i + 1, j + 1); return b1 && b2 && b3 && b4 && b5 && b6 && b7 && b8; } private bool FloodFill(int x, int y,List<Int16Double> ret,FlagMap2d flag) { ret.Clear(); Queue<Int16Double> queue = new Queue<Int16Double>(); ret.Add(new Int16Double(x, y)); float pvalue = bmp.GetPixel(x, y); flag.SetFlagOn(x, y, VISITED); queue.Enqueue(new Int16Double(x, y)); while (queue.Count != 0) { Int16Double p = queue.Dequeue(); for (int i = 0; i < 8; i++) { int tx = p.X + Delta[i].X; int ty = p.Y + Delta[i].Y; if(InRange(tx, ty)) { byte f= flag.GetFlagOn(tx,ty); if(f==PROCESSED) return false; else { bool minum = false; if (IncludePredicate(tx, ty, pvalue,ref minum) && f == UNPROCESSED) { if (minum) return false; Int16Double t = new Int16Double(tx, ty); queue.Enqueue(t); flag.SetFlagOn(tx, ty, VISITED); ret.Add(t); } } } } } return true; } private bool InRange(int tx, int ty) { return tx >= 0 && tx < bmp.width && ty >= 0 && ty < bmp.height; } private bool IncludePredicate(int x, int y, float pv, ref bool min) { float v = bmp.GetPixel(x, y); if (pv < v) min = true; return pv - v <= torlerance; } private void MarkAll(List<Int16Double> ret, byte v,FlagMap2d flag) { for (int i = 0; i < ret.Count; i++) { flag.SetFlagOn(ret[i].X, ret[i].Y, v); } } }
算法输出预览
原图和输出图像如下表所示：
 |
 |
|
| 原图 | 本文实现的输出 | IJ的输出 |
代码地址：https://github.com/chnhideyoshi/SeededGrow2d/tree/master/FindMaxima
爬网的太疯狂了，转载本文要注明出处啊：http://www.cnblogs.com/chnhideyoshi/