Re: Re: Vector Union

*To*: mathgroup at smc.vnet.net*Subject*: [mg23263] Re: [mg23141] Re: Vector Union*From*: Hartmut Wolf <hwolf at debis.com>*Date*: Sat, 29 Apr 2000 22:05:05 -0400 (EDT)*Organization*: debis Systemhaus*References*: <8djl4o$7jh@smc.vnet.net> <200004200720.DAA14499@smc.vnet.net>*Sender*: owner-wri-mathgroup at wolfram.com

Jens-Peer Kuska schrieb: > > Hi Russel, > > I overlooked you question, so my answer will be a bit late. > I can't give you a n-dimensional version but the extension is > easy if you know the dimension. The fastest version to find the > unique $n$-dimensional points in a set is to use a tree of $2^n$ > > In three dimensions one can use > > $SamePrecision = 0.1*^-6; > > pointCount = 0; > > getIndex[p_List, p1_List] := > Module[{index}, > index = 1; > If[p[[1]] < p1[[1]], index += 1]; > If[p[[2]] < p1[[2]], index += 2]; > If[p[[3]] < p1[[3]], index += 4]; > index > ] > > (* pointCount is fluid from MakeOctTree[] *) > > OctTreeInsert[p_List, {}] := {p, pointCount++, Table[{}, {8}], {}} > > OctTreeInsert[p_List, {p1_, i_Integer, l_, poly_List}] := > > Module[{index, delta, l1}, > delta = p - p1; > If[Abs[Dot[delta, delta]] < $SamePrecision, > Return[{p1, i, l, poly}] > ]; > index = getIndex[p, p1]; > l1 = OctTreeInsert[p, l[[index]]]; > {p1, i, ReplacePart[l, l1, index], poly} > ] > > and for the points pnts > > pointCount = 0; > > uniquePoints=Cases[Fold[OctTreeInsert[#2, #1] &, {}, > pnts], {_?NumericQ, _?NumericQ, _?NumericQ}, Infinity] > > The speed gain depends on the number of points. MathGL3d uses an octtree > to make a indexed point set form the vertex data. Typical for 10^4 > points the > usual search is to slow. A second advantage may be, that the octtree is > incremental. You can add a k+1 point to your set and find quickly if it > is > already in the set. The Union/Sort[] versions need always the full set > of > points and can't be incremented without a new Sort[] call that consume > (k+1)*Log[k+1] > operations where the octtree need only Log[k+1] operations. > > For small point sets the overhead of the tree creation is to huge and > the other versions maye work well. Also the memory for the tree is a > constrain, because for a surface a maximum of 4 links may be used and > four links are useless memory. In higher dimensions this becomes more > dramatical. > > How ever for more than 10^4 points (with MathGL3d I have tryed 10^7) > the octtree works excelent. > > Regards > Jens > > Russell Towle wrote: > > > > The problem of finding a function to return the union of a list of > > n-vectors which is at once *fast* and yet *unfailing* has beguiled me for > > several years. I contrived a slow but unfailing 'vunion' function about > > five years ago; made a not-so-slow and unfailing version a year ago; and > > then found it to occasionally fail when Chop was applied to the vectors > > without a second argument. > > Dear Jens, not meant as a point of critique -- a quick answer often is of higher value -- but as general remark for use of data structures with Mathematica: I observed that with your QuadTree program, when a new node is added, the whole path from the added leaf up to the root has to be rebuilt. This might seem to be uneconomic. Normally, when programming in C e.g., one uses pointers, only one of which is reset at the insertion point. In Mathematica we do not have pointers, however there is a surrogate for that: i.e. (unevaluated) Symbols. So I reworked your program using this kluge: In[9]:= $SamePrecisionSquared = (0.5 10^-6)^2; In[10]:= ClearAll[OctTreeInsert, MakeNode, getIndex2] In[11]:= getIndex2[p_, p1_] := {1, 2, 4}.(Thread[p < p1] /. {True -> 1, False -> 0}) + 1 In[12]:= MakeNode[p_] := ( pointCount++; {p, Table[Unique[OctTree], {8}]}) In[13]:= OctTreeInsert2[p_, node_Symbol] := node = MakeNode[p] In[14]:= OctTreeInsert2[p_, thisnode : {p1_, daughters_}] := Module[{delta = p - p1}, If[Abs[delta.delta] < $SamePrecisionSquared, Return[thisnode]]; OctTreeInsert2[p, daughters[[getIndex2[p, p1]]]] ] Trial run (I used the vects as I defined in my answer to Russell): In[15]:= pointCount = 0; Clear[tree]; Map[OctTreeInsert2[#, tree] &, vects]; // Timing Out[16]= {13.48 Second, Null} In[17]:= uniquePoints2 = Cases[tree, {_?NumericQ, _?NumericQ, _?NumericQ}, Infinity]; // Timing Out[17]= {0.03 Second, Null} [Your version (having removed the scrap) with the same data was {17.075 Second, Null}.] In[18]:= {pointCount, Length[uniquePoints2]} Out[18]= {158, 158} In[19]:= Shallow[tree, 6] Out[19]//Shallow= {{-7.105233285472004`*^-7, 1.2297443354282358`*^-6, 2.1442483040006357`*^-6}, {{{4.708766454944861`, 4.786779659984509`, 5.928100903105975`}, {{<<2>>}, {<<2>>}, {<<2>>}, OctTree$11990, {<<2>>}, OctTree$11992, {<<2>>}, {<<2>>}}}, {{-0.31739244789510035`, 0.7331739702018931`, 0.11293200180872479`}, {OctTree$7201, {<<2>>}, {<<2>>}, OctTree$7204, OctTree$7205, {<<2>>}, OctTree$7207, {<<2>>}}}, {{-5.405432541097434`*^-7, 4.771981519312292`*^-7, 2.248561517851597`*^-6}, {OctTree$3393, OctTree$3394, {<<2>>}, OctTree$3396, OctTree$3397, OctTree$3398, OctTree$3399, OctTree$3400}}, {{-7.29552608995409`*^-7, 7.327158225945956`*^-7, 2.225459337623892`*^-6}, {OctTree$3455, OctTree$3456, OctTree$3457, {<<2>>}, OctTree$3459, OctTree$3460, OctTree$3461, {<<2>>}}}, {{-5.368389831293919`*^-7, 1.7807576570054204`*^-6, 1.9614778091772764`*^-6}, {OctTree$3316, OctTree$3317, OctTree$3318, OctTree$3319, {<<2>>}, OctTree$3321, OctTree$3322, OctTree$3323}}, {{-6.511382857212354`, 0.47816750596348706`, -2.3640961082723035`}, {{<<2>>}, {<<2>>}, {<<2>>},{<<2>>}, {<<2>>}, OctTree$9630, OctTree$9631, OctTree$9632}}, {{-6.708995711564143`*^-7, 5.623158954916112`*^-7, 2.0236280131769387`*^-6}, {OctTree$3424, OctTree$3425, OctTree$3426, OctTree$3427, {<<2>>}, {<<2>>}, {<<2>>}, OctTree$3431}}, {{-7.132906373155661`*^-7, 5.444596867953894`*^-7, 2.0197355418986647`*^-6}, {OctTree$3435, OctTree$3436, OctTree$3437, {<<2>>}, OctTree$3439, {<<2>>}, OctTree$3441, {<<2>>}}}}} The undefined symbols are Nil-Pointers. Kind regards, Hartmut

**References**:**Re: Vector Union***From:*Jens-Peer Kuska <kuska@informatik.uni-leipzig.de>