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>
- Re: Vector Union