       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[] < p1[], index += 1];
>        If[p[] < p1[], index += 2];
>        If[p[] < p1[], 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:= \$SamePrecisionSquared = (0.5 10^-6)^2;

In:= ClearAll[OctTreeInsert, MakeNode, getIndex2]

In:= getIndex2[p_, p1_] :=
{1, 2, 4}.(Thread[p < p1] /. {True -> 1, False -> 0}) + 1

In:= MakeNode[p_] := ( pointCount++; {p, Table[Unique[OctTree],
{8}]})

In:= OctTreeInsert2[p_, node_Symbol] := node = MakeNode[p]

In:= 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:= pointCount = 0; Clear[tree];
Map[OctTreeInsert2[#, tree] &, vects]; // Timing

Out= {13.48 Second, Null}

In:= uniquePoints2 =
Cases[tree, {_?NumericQ, _?NumericQ, _?NumericQ}, Infinity]; //
Timing

Out= {0.03 Second, Null}

[Your version (having removed the scrap) with the same data was {17.075
Second, Null}.]

In:= {pointCount, Length[uniquePoints2]}
Out= {158, 158}

In:= Shallow[tree, 6]

Out//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:
• Prev by Date: Interpolation bugs (feature?) - work arounds & efficiency for large n?
• Next by Date: User's Manual?
• Previous by thread: Re: Vector Union
• Next by thread: breaking up lists into intervals