[Date Index]
[Thread Index]
[Author Index]
RE: Mathematica Code to Build and Access KD Trees
*To*: mathgroup at smc.vnet.net
*Subject*: [mg45887] RE: [mg45713] Mathematica Code to Build and Access KD Trees
*From*: "Wolf, Hartmut" <Hartmut.Wolf at t-systems.com>
*Date*: Thu, 29 Jan 2004 05:34:32 -0500 (EST)
*Sender*: owner-wri-mathgroup at wolfram.com
Richard,
here is something I did last Friday, when I had a day of rest after an
operation. I did it just for fun, while reading (my now old version of)
"Algorithms" by Robert Sedgewick (from which I took the suggestion for the
pictures). Also found an error in Gonnet (and some subleties which are not
reported here). The Mathematica implementation is mine; including all the
algorithms (but they certainly are reinventions of something out elsewhere).
(*******************************************************)
Off[General::"spell1"]
SeedRandom[123456]
(* bascic definition*)
K = 2; (* the dimension, a global variable *)
Clear[makeKDnode]
makeKDnode[key_] /; Length[key] == K :=
Module[{leftKD, rightKD}, {key, leftKD, rightKD}]
tt = makeKDnode[{1., 2.}]
KD::"insert" = "key `1` already present";
ClearAll[insertKD]
insertKD[t_, indx_, key_] :=
If[Head[t] === Symbol, t = makeKDnode[key],
If[key == t[[1]],
Message[KD::"insert", key],
If[key[[indx]] >
t[[1, indx]],
insertKD[t[[-1]],
Mod[indx + 1, K, 1], key],
insertKD[t[[-2]],
Mod[indx + 1, K, 1], key]
]]]
Clear[tt]
Remove["leftKD*", "rightKD*"]
pts = Table[{Random[], Random[]}, {7}]
insertKD[tt, 1, #] & /@ pts
tt
(* drawing the regions *)
draw[t_, indx_, {{xlow_, xhi_}, {ylow_, yhi_}}] :=
Module[{indx2 = Mod[indx + 1, 2, 1]}, {Point[t[[1]]],
Line[
ReplacePart[{{xlow, ylow}, {xhi, yhi}},
t[[1, indx]], {{1, indx}, {2, indx}}]],
If[Head[t[[-#]]] =!= Symbol,
draw[t[[-#]], indx2,
ReplacePart[{{xlow, xhi}, {ylow, yhi}},
t[[1, indx]], {indx, #}]], {}] & /@ {2, 1}
}]
Show[Graphics[{PointSize[.02], draw[tt, 1, {{0, 1}, {0, 1}}]}],
PlotRange -> {{0, 1}, {0, 1}}, AspectRatio -> Automatic, Frame -> True]
(* balanced kd-tree *)
balanceKD[pts_, indx_] :=
If[pts == {}, {},
Block[{mid = Ceiling[Length[pts]/2], order = Ordering[pts[[All,
indx]]]},
Join[pts[[order[[{mid}]]]],
balanceKD[pts[[Take[order, {1, mid - 1}]]], Mod[indx + 1, K, 1]],
balanceKD[pts[[Take[order, {mid + 1, -1}]]], Mod[indx + 1, K, 1]]]]]
pts = Table[{Random[], Random[]}, {170}];
Clear[tt]
insertKD[tt, 1, #] & /@ pts;
pts2 = balanceKD[pts, 1];
Clear[tt2]
insertKD[tt2, 1, #] & /@ pts2;
gg = Show[Graphics[{PointSize[.02], draw[tt, 1, {{0, 1}, {0, 1}}]}],
PlotRange -> {{0, 1}, {0, 1}}, AspectRatio -> Automatic, Frame -> True,
DisplayFunction -> Identity]
gg2 = Show[Graphics[{PointSize[.02], draw[tt2, 1, {{0, 1}, {0, 1}}]}],
PlotRange -> {{0, 1}, {0, 1}}, AspectRatio -> Automatic, Frame -> True,
DisplayFunction -> Identity]
Show[GraphicsArray[{{gg, gg2}}]]
(* range search *)
rangeSearchKD[t_, indx_, lowkey_, uppkey_] :=
If[Head[t] =!= Symbol,
Join[If[And @@ Thread[lowkey <= t[[1]] <= uppkey], {t[[1]]}, {}],
If[lowkey[[indx]] <= t[[1, indx]],
rangeSearchKD[t[[-2]], Mod[indx + 1, K, 1], lowkey, uppkey], {}],
If[t[[1, indx]] < uppkey[[indx]],
rangeSearchKD[t[[-1]], Mod[indx + 1, K, 1], lowkey, uppkey], {}]],
{}]
region = Sequence[{.2, .4}, {.35, .6}];
rpts = rangeSearchKD[tt, 1, region]
rangeSearchKD[tt2, 1, region]
Sort[%] == Sort[%%]
Show[Graphics[{GrayLevel[.6], Rectangle[region]}], gg,
Graphics[{Hue[0], PointSize[.02], Point /@ rpts}],
DisplayFunction -> $DisplayFunction, Options[gg]]
Show[Graphics[{GrayLevel[.6], Rectangle[region]}], gg2,
Graphics[{Hue[0], PointSize[.02], Point /@ rpts}],
DisplayFunction -> $DisplayFunction, Options[gg2]]
(* Nearest point search *)
(* search cell *)
findCellKD[t_, indx_, lowkey_, uppkey_, pnt_] :=
If[Head[t] === Symbol, Return[{lowkey, uppkey}],
If[pnt[[indx]] <= t[[1, indx]],
findCellKD[t[[-2]], Mod[indx + 1, K, 1], lowkey,
ReplacePart[uppkey, t[[1, indx]], indx], pnt],
findCellKD[t[[-1]], Mod[indx + 1, K, 1],
ReplacePart[lowkey, t[[1, indx]], indx], uppkey, pnt]]]
cell = findCellKD[tt, 1, {0, 0}, {1, 1}, {.4, .5}]
Show[Graphics[{Hue[0, .5, 1], Rectangle @@ cell}], gg,
DisplayFunction -> $DisplayFunction, Options[gg]]
Show[Graphics[{Hue[0, .5, 1],
Rectangle @@ findCellKD[tt, 1, {0, 0}, {1, 1}, {0.62, 0.51}]}], gg,
DisplayFunction -> $DisplayFunction, Options[gg]]
Show[Graphics[{Hue[0, .5, 1],
Rectangle @@ findCellKD[tt, 1, {0, 0}, {1, 1}, {0.27, 0.51}]}], gg,
DisplayFunction -> $DisplayFunction, Options[gg]]
Show[Graphics[{Hue[0, .5, 1],
Rectangle @@ findCellKD[tt, 1, {0, 0}, {1, 1}, {-1., -1.}]}], gg,
DisplayFunction -> $DisplayFunction, Options[gg]]
Show[Graphics[{Hue[0, .5, 1],
Rectangle @@ findCellKD[tt, 1, {0, 0}, {1, 1}, {-1., 2.}]}], gg,
DisplayFunction -> $DisplayFunction, Options[gg]]
(* nearest neighbor *)
nearestPointSearchKD[tree_, pnt_] :=
Module[{nspKD, lowkey = {-\[Infinity], -\[Infinity]},
uppkey = {\[Infinity], \[Infinity]}, d = \[Infinity], d2, nearestPnt},
nspKD[t_, indx_] := If[Head[t] =!= Symbol,
Block[{pp = t[[1]]},
If[(d2 = Sqrt[(pp - pnt).(pp - pnt)]) < d,
d = d2; nearestPnt = pp; lowkey = pnt - d; uppkey = pnt + d;
AppendTo[dlist, d]];
If[lowkey[[indx]] <= pp[[indx]],
nspKD[t[[-2]], Mod[indx + 1, K, 1]]];
If[pp[[indx]] < uppkey[[indx]],
nspKD[t[[-1]], Mod[indx + 1, K, 1]]]
]];
dlist = {};
nspKD[tree, 1];
nearestPnt]
pt = {.5, .5}
nearestPointSearchKD[tt, pt]
dlist
With[{dds = Sqrt[(# - pt).(# - pt)] & /@ pts},
i = Position[dds, Min[dds]][[1, 1]]; ptx = pts[[i]] ]
Show[gg, Graphics[{Hue[.3], PointSize[.02], Point[pt], Hue[0], Point[ptx]}],
DisplayFunction -> $DisplayFunction];
pt = {.5, .48}
nearestPointSearchKD[tt, pt]
dlist
With[{dds = Sqrt[(# - pt).(# - pt)] & /@ pts},
i = Position[dds, Min[dds]][[1, 1]]; ptx = pts[[i]] ]
Show[gg, Graphics[{Hue[.3], PointSize[.02], Point[pt], Hue[0], Point[ptx]}],
DisplayFunction -> $DisplayFunction];
pt = {.1, .5}
nearestPointSearchKD[tt, pt]
dlist
With[{dds = Sqrt[(# - pt).(# - pt)] & /@ pts},
i = Position[dds, Min[dds]][[1, 1]]; ptx = pts[[i]] ]
Show[gg, Graphics[{Hue[.3], PointSize[.02], Point[pt], Hue[0], Point[ptx]}],
DisplayFunction -> $DisplayFunction];
pt = {.15, .5}
nearestPointSearchKD[tt, pt]
dlist
With[{dds = Sqrt[(# - pt).(# - pt)] & /@ pts},
i = Position[dds, Min[dds]][[1, 1]]; ptx = pts[[i]] ]
Show[gg, Graphics[{Hue[.3], PointSize[.02], Point[pt], Hue[0], Point[ptx]}],
DisplayFunction -> $DisplayFunction];
(*******************************************************)
etc. you see how you may play with. Just paste everything between the marker
lines into a notebook an press the button.
I was lazy to write a makeBalancedTree, just reordered the points such that
mapping my insertKD would result in a balanced tree. (And such my nearest
Point Search keeps valid, it must be changed a (very) little in case of a
balanced insert (can you see it?)). BTW in most cases a random shuffle of
the point would do. Other functions have to be implemented, but it all
depends on your real needs, and I will not go further here. Perhaps a
n-nearest Points Search might be nice, and fixing up nearestPointSearch,
such that if the point in question is part of the structure it wouldn't find
itself, but the next nearest. But I think this is relatively easy. As you
see, there is no payload at the "points" == "keys"; and I did not delve into
partial key search. Finally I spend no thoughts on performance and storage
use (beyond my routine). At the end I have no idea what a "vintage point
tree" is, and I did not read the original papers.
As stated, this is all play work, for production got to C.
--
Hartmut
>-----Original Message-----
>From: Richard Palmer [mailto:dickp at bellatlantic.net]
To: mathgroup at smc.vnet.net
>Sent: Tuesday, January 27, 2004 7:43 PM
>To: Wolf, Hartmut
>Subject: [mg45887] Re: [mg45713] Mathematica Code to Build and Access KD Trees
>
>
>I'm doing a prototype of something. If I have to turn it into
>a production
>version. I'm willing to do the work. I have found the references you
>mentioned earlier. I couldn't find any math code that
>addressed any version
>of the problem. Thanks for the heads up!
>
>Richard
>
>On 1/27/04 12:10 PM, "Wolf, Hartmut"
><Hartmut.Wolf at t-systems.com> wrote:
>
>>
>>> -----Original Message-----
>>> From: Richard Palmer [mailto:dickp at bellatlantic.net]
To: mathgroup at smc.vnet.net
>>> Sent: Wednesday, January 21, 2004 10:55 AM
>>> To: mathgroup at smc.vnet.net
>>> Subject: [mg45887] [mg45713] Mathematica Code to Build and Access KD Trees
>>>
>>>
>>> Anybody know where some published code exists.
>>> --
>>>
>>>
>>
>> Some code, hm, you must be joking!
>>
>> There are many different flavors of k-d trees, much depending on the
>> application.
>>
>> In my private library I find something in
>>
>> - G.H.Gonnet: Handbook of algorithms and Data Structures,
>1st ed. 1984
>>
>> - Robert Sedgewick: Algorithms
>> (there are quite a lot versions out now, I have the oldest
>one from 1983,
>> with coding in Pascal)
>>
>> - Steven S.Skiena: The Algorithm Design Manual, 1998
>>
>> Gonnet and Skiena give you pointers to the literature.
>>
>> There also is something to be found at the WRI site:
>> http://library.wolfram.com/infocenter/MathSource/462/
>>
>> But this is only a study of an original algorithm as given
>by J.L.Bently
>> (ACM Transactions, Vol.3 No.3, 1977) using Mathematica as a rapid
>> prototyping language. For production purpose it's useless.
>>
>>
>> Dealing with functional data structures (and keep
>performance) is not easy.
>> Perhaps you might like to look at Jens-Peer Kuska's
>contribution to this
>> newsgroup
>> http://forums.wolfram.com/mathgroup/archive/2000/Apr/msg00304.html
>>
>> and also at my response to it
>> http://forums.wolfram.com/mathgroup/archive/2000/Apr/msg00424.html
>>
>> The example there is for OctTrees, but similar principles
>also apply to
>> kd-trees in Mathematica. If you want to study problems of
>that kind, look
>> into
>>
>> - Chris Okasaki: Purely Functional Data Structures, 1998
>> (this is for ML and Haskell)
>>
>> However, if you intend to do more than a toy application you
>preferably
>> should code in C or buy or steal a library.
>>
>>
>> If you want to get reasonable help with Mathematica coding,
>however, you
>> must show more from your application problem, and show your
>tries, i.e. give
>> a challenge to the community. Your quest doesn't.
>>
>>
>> --
>> Hartmut Wolf
>>
>>
>>
>
>--
>
>
Prev by Date:
**Re: RE: GridPlot**
Next by Date:
** Re: MathML Namespace + Mathematica**
Previous by thread:
**RE: Mathematica Code to Build and Access KD Trees**
Next by thread:
**copying**
| |