MathGroup Archive 2000

[Date Index] [Thread Index] [Author Index]

Search the Archive

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:
  • 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