读取shp文件建立二叉树索引的基本思路分析：
1 /// <summary>
2 /// Generates a spatial index for a specified shape file.
3 /// 根据特定的shp文件构建空间索引
4 /// </summary>
5 /// <param name="filename"></param>
6 private QuadTree CreateSpatialIndex(string filename)
7 {
8 List<QuadTree.BoxObjects> objList = new List<QuadTree.BoxObjects>();
9 //Convert all the geometries to boundingboxes
10 uint i = 0;
11 foreach (BoundingBox box in GetAllFeatureBoundingBoxes())
12 {
13 if (!double.IsNaN(box.Left) && !double.IsNaN(box.Right) && !double.IsNaN(box.Bottom) &&
14 !double.IsNaN(box.Top))
15 {
16 QuadTree.BoxObjects g = new QuadTree.BoxObjects();
17 g.box = box;
18 g.ID = i;
19 objList.Add(g);
20 i++;
21 }
22 }
23
24 Heuristic heur;
25 heur.maxdepth = (int) Math.Ceiling(Math.Log(GetFeatureCount(), 2));
26 heur.minerror = 10;
27 heur.tartricnt = 5;
28 heur.mintricnt = 2;
29 return new QuadTree(objList, 0, heur);
30 }
1 // Copyright 2005, 2006 - Morten Nielsen (www.iter.dk)
2 //
3 // This file is part of SharpMap.
4 // SharpMap is free software; you can redistribute it and/or modify
5 // it under the terms of the GNU Lesser General Public License as published by
6 // the Free Software Foundation; either version 2 of the License, or
7 // (at your option) any later version.
8 //
9 // SharpMap is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 // GNU Lesser General Public License for more details.
13
14 // You should have received a copy of the GNU Lesser General Public License
15 // along with SharpMap; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17
18 using System;
19 using System.Collections.Generic;
20 using System.Collections.ObjectModel;
21 using System.IO;
22 using SharpMap.Geometries;
23
24 namespace SharpMap.Utilities.SpatialIndexing
25 {
26 /// <summary>
27 /// Heuristics used for tree generation
28 /// </summary>
29 public struct Heuristic
30 {
31 /// <summary>
32 /// Maximum tree depth
33 /// </summary>
34 public int maxdepth;
35
36 /// <summary>
37 /// Minimum Error metric ?the volume of a box + a unit cube.
38 /// The unit cube in the metric prevents big boxes that happen to be flat having a zero result and muddling things up.
39 /// </summary>
40 public int minerror;
41
42 /// <summary>
43 /// Minimum object count at node
44 /// </summary>
45 public int mintricnt;
46
47 /// <summary>
48 /// Target object count at node
49 /// </summary>
50 public int tartricnt;
51 }
52
53 /// <summary>
54 /// Constructs a Quad-tree node from a object list and creates its children recursively
55 /// 二叉树？四叉树？
56 /// </summary>
57 public class QuadTree : IDisposable
58 {
59 private BoundingBox _box;
60 private QuadTree _child0;
61 private QuadTree _child1;
62
63 /// <summary>
64 /// Nodes depth in a tree
65 /// </summary>
66 protected uint _Depth;
67
68 /// <summary>
69 /// Node ID
70 /// </summary>
71 public uint? _ID;
72
73 private List<BoxObjects> _objList;
74
75 /// <summary>
76 /// Creates a node and either splits the objects recursively into sub-nodes, or stores them at the node depending on the heuristics.
77 /// Tree is built top->down
78 /// </summary>
79 /// <param name="objList">Geometries to index</param>
80 /// <param name="depth">Current depth of tree树的深度</param>
81 /// <param name="heurdata">Heuristics data</param>
82 public QuadTree(List<BoxObjects> objList, uint depth, Heuristic heurdata)
83 {
84 _Depth = depth;
85
86 _box = objList[0].box;
87 for (int i = 0; i < objList.Count; i++)
88 _box = _box.Join(objList[i].box);
89
90 // test our build heuristic探试性的- if passes, make children
91 if (depth < heurdata.maxdepth && objList.Count > heurdata.mintricnt &&
92 (objList.Count > heurdata.tartricnt || ErrorMetric(_box) > heurdata.minerror))
93 {
94 List<BoxObjects>[] objBuckets = new List<BoxObjects>[2]; // buckets of geometries大量的几何对象，一桶几何对象
95 objBuckets[0] = new List<BoxObjects>();
96 objBuckets[1] = new List<BoxObjects>();
97
98 uint longaxis = _box.LongestAxis; // longest axis
99 double geoavg = 0; // geometric average - midpoint of ALL the objects
100
101 // go through all bbox and calculate the average of the midpoints
102 double frac = 1.0f / objList.Count;
103 for (int i = 0; i < objList.Count; i++)
104 geoavg += objList[i].box.GetCentroid()[longaxis] * frac;
105
106 // bucket bbox based on their midpoint's side of the geo average in the longest axis
107 for (int i = 0; i < objList.Count; i++)
108 objBuckets[geoavg > objList[i].box.GetCentroid()[longaxis] ? 1 : 0].Add(objList[i]);
109
110 //If objects couldn't be splitted, just store them at the leaf
111 //TODO: Try splitting on another axis
112 if (objBuckets[0].Count == 0 || objBuckets[1].Count == 0)
113 {
114 _child0 = null;
115 _child1 = null;
116 // copy object list
117 _objList = objList;
118 }
119 else
120 {
121 // create new children using the buckets
122 _child0 = new QuadTree(objBuckets[0], depth + 1, heurdata);
123 _child1 = new QuadTree(objBuckets[1], depth + 1, heurdata);
124 }
125 }
126 else
127 {
128 // otherwise the build heuristic failed, this is
129 // set the first child to null (identifies a leaf)
130 _child0 = null;
131 _child1 = null;
132 // copy object list
133 _objList = objList;
134 }
135 }
136
137 /// <summary>
138 /// This instantiator is used internally for loading a tree from a file
139 /// </summary>
140 private QuadTree()
141 {
142 }
143
144 #region Read/Write index to/from a file
145
146 private const double INDEXFILEVERSION = 1.0;
147
148 /// <summary>
149 /// Loads a quadtree from a file
150 /// </summary>
151 /// <param name="filename"></param>
152 /// <returns></returns>
153 public static QuadTree FromFile(string filename)
154 {
155 FileStream fs = new FileStream(filename, FileMode.Open, FileAccess.Read);
156 BinaryReader br = new BinaryReader(fs);
157 if (br.ReadDouble() != INDEXFILEVERSION) //Check fileindex version
158 {
159 fs.Close();
160 fs.Dispose();
161 throw new ObsoleteFileFormatException(
162 "Invalid index file version. Please rebuild the spatial index by either deleting the index");
163 }
164 QuadTree node = ReadNode(0, ref br);
165 br.Close();
166 fs.Close();
167 return node;
168 }
169
170 /// <summary>
171 /// Reads a node from a stream recursively
172 /// </summary>
173 /// <param name="depth">Current depth</param>
174 /// <param name="br">Binary reader reference</param>
175 /// <returns></returns>
176 private static QuadTree ReadNode(uint depth, ref BinaryReader br)
177 {
178 QuadTree node = new QuadTree();
179 node._Depth = depth;
180 node.Box = new BoundingBox(br.ReadDouble(), br.ReadDouble(), br.ReadDouble(), br.ReadDouble());
181 bool IsLeaf = br.ReadBoolean();
182 if (IsLeaf)
183 {
184 int FeatureCount = br.ReadInt32();
185 node._objList = new List<BoxObjects>();
186 for (int i = 0; i < FeatureCount; i++)
187 {
188 BoxObjects box = new BoxObjects();
189 box.box = new BoundingBox(br.ReadDouble(), br.ReadDouble(), br.ReadDouble(), br.ReadDouble());
190 box.ID = (uint)br.ReadInt32();
191 node._objList.Add(box);
192 }
193 }
194 else
195 {
196 node.Child0 = ReadNode(node._Depth + 1, ref br);
197 node.Child1 = ReadNode(node._Depth + 1, ref br);
198 }
199 return node;
200 }
201
202 /// <summary>
203 /// Saves the Quadtree to a file
204 /// </summary>
205 /// <param name="filename"></param>
206 public void SaveIndex(string filename)
207 {
208 FileStream fs = new FileStream(filename, FileMode.Create);
209 BinaryWriter bw = new BinaryWriter(fs);
210 bw.Write(INDEXFILEVERSION); //Save index version
211 SaveNode(this, ref bw);
212 bw.Close();
213 fs.Close();
214 }
215
216 /// <summary>
217 /// Saves a node to a stream recursively
218 /// </summary>
219 /// <param name="node">Node to save</param>
220 /// <param name="sw">Reference to BinaryWriter</param>
221 private void SaveNode(QuadTree node, ref BinaryWriter sw)
222 {
223 //Write node boundingbox
224 sw.Write(node.Box.Min.X);
225 sw.Write(node.Box.Min.Y);
226 sw.Write(node.Box.Max.X);
227 sw.Write(node.Box.Max.Y);
228 sw.Write(node.IsLeaf);
229 if (node.IsLeaf)
230 {
231 sw.Write(node._objList.Count); //Write number of features at node
232 for (int i = 0; i < node._objList.Count; i++) //Write each featurebox
233 {
234 sw.Write(node._objList[i].box.Min.X);
235 sw.Write(node._objList[i].box.Min.Y);
236 sw.Write(node._objList[i].box.Max.X);
237 sw.Write(node._objList[i].box.Max.Y);
238 sw.Write(node._objList[i].ID);
239 }
240 }
241 else if (!node.IsLeaf) //Save next node
242 {
243 SaveNode(node.Child0, ref sw);
244 SaveNode(node.Child1, ref sw);
245 }
246 }
247
248 internal class ObsoleteFileFormatException : Exception
249 {
250 /// <summary>
251 /// Exception thrown when layer rendering has failed
252 /// </summary>
253 /// <param name="message"></param>
254 public ObsoleteFileFormatException(string message)
255 : base(message)
256 {
257 }
258 }
259
260 #endregion
261
262 /// <summary>
263 /// Determines whether the node is a leaf (if data is stored at the node, we assume the node is a leaf)
264 /// </summary>
265 public bool IsLeaf
266 {
267 get { return (_objList != null); }
268 }
269
270 ///// <summary>
271 ///// Gets/sets the list of objects in the node
272 ///// </summary>
273 //public List<SharpMap.Geometries.IGeometry> ObjList
274 //{
275 // get { return _objList; }
276 // set { _objList = value; }
277 //}
278
279 /// <summary>
280 /// Gets/sets the Axis Aligned Bounding Box
281 /// </summary>
282 public BoundingBox Box
283 {
284 get { return _box; }
285 set { _box = value; }
286 }
287
288 /// <summary>
289 /// Gets/sets the left child node
290 /// </summary>
291 public QuadTree Child0
292 {
293 get { return _child0; }
294 set { _child0 = value; }
295 }
296
297 /// <summary>
298 /// Gets/sets the right child node
299 /// </summary>
300 public QuadTree Child1
301 {
302 get { return _child1; }
303 set { _child1 = value; }
304 }
305
306 /// <summary>
307 /// Gets the depth of the current node in the tree
308 /// </summary>
309 public uint Depth
310 {
311 get { return _Depth; }
312 }
313
314 #region IDisposable Members
315
316 /// <summary>
317 /// Disposes the node
318 /// </summary>
319 public void Dispose()
320 {
321 //this._box = null;
322 _child0 = null;
323 _child1 = null;
324 _objList = null;
325 }
326
327 #endregion
328
329 /// <summary>
330 /// Calculate the floating point error metric
331 /// </summary>
332 /// <returns></returns>
333 public double ErrorMetric(BoundingBox box)
334 {
335 Point temp = new Point(1, 1) + (box.Max - box.Min);
336 return temp.X * temp.Y;
337 }
338
339 /// <summary>
340 /// Searches the tree and looks for intersections with the boundingbox 'bbox'
341 /// </summary>
342 /// <param name="box">Boundingbox to intersect with</param>
343 public Collection<uint> Search(BoundingBox box)
344 {
345 Collection<uint> objectlist = new Collection<uint>();
346 IntersectTreeRecursive(box, this, ref objectlist);
347 return objectlist;
348 }
349
350 /// <summary>
351 /// Recursive function that traverses the tree and looks for intersections with a boundingbox
352 /// </summary>
353 /// <param name="box">Boundingbox to intersect with</param>
354 /// <param name="node">Node to search from</param>
355 /// <param name="list">List of found intersections</param>
356 private void IntersectTreeRecursive(BoundingBox box, QuadTree node, ref Collection<uint> list)
357 {
358 if (node.IsLeaf) //Leaf has been reached
359 {
360 foreach (BoxObjects boxObject in node._objList)
361 {
362 if (box.Intersects(boxObject.box))
363 list.Add(boxObject.ID);
364
365 }
366 /*
367 for (int i = 0; i < node._objList.Count; i++)
368 {
369 list.Add(node._objList[i].ID);
370 }
371 */
372 }
373 else
374 {
375 if (node.Box.Intersects(box))
376 {
377 IntersectTreeRecursive(box, node.Child0, ref list);
378 IntersectTreeRecursive(box, node.Child1, ref list);
379 }
380 }
381 }
382
383 #region Nested type: BoxObjects
384
385 /// <summary>
386 /// BoundingBox and Feature ID structure used for storing in the quadtree
387 /// 包围盒和要素ID数据结构，用于在四叉树中存储
388 /// </summary>
389 public struct BoxObjects
390 {
391 /// <summary>
392 /// Boundingbox
393 /// 包围盒
394 /// </summary>
395 public BoundingBox box;
396
397 /// <summary>
398 /// Feature ID
399 /// 要素ID
400 /// </summary>
401 public uint ID;
402 }
403
404 #endregion
405 }
406 }
今天看了下DotSpatial中读取Shp的代码，真正实现的四叉树空间索引。