[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Center of mass???
- Subject: Re: Center of mass???
- From: "Dick Jackson" <dickjackson(at)dfanning.com>
- Date: Tue, 09 Nov 1999 19:00:00 GMT
- Newsgroups: comp.lang.idl-pvwave
- Organization: @Home Network Canada
- References: <3826E4D0.6E84F53B@mrc.ks.se> <MPG.email@example.com> <38279336.C13C76F8@ucsd.edu> <MPG.firstname.lastname@example.org>
- Xref: news.doit.wisc.edu comp.lang.idl-pvwave:17191
Paul Hick (email@example.com) writes:
> Reasoning by analogy to the 2D case, this should work, I think:
> xcm = Total( Total(Total(array,3),2) * Indgen(s)) / totalMass
> ycm = Total( Total(Total(array,3),1) * Indgen(s)) / totalMass
> zcm = Total( Total(Total(array,2),1) * Indgen(s)) / totalMass
Right, but as a wise man once told me (Dr. Coyote, or something like that),
the fastest dimension to run across with things like Total is usually the
second dimension. Empirical testing confirms this, so I propose the
following, now extended to do 2D or 3D. I also pulled the "/ totalMass"
inside a bit to keep the numbers closer to 1, lessen the possibility of
overflow and perhaps maintain more precision.
FUNCTION CenterOfMass, array
s = Size(array, /Dimensions)
totalMass = Total(array)
CASE Size(array, /N_Dimensions) OF
xcm = Total(Total(array,2) / totalMass * Indgen(s))
ycm = Total(Total(array,1) / totalMass * Indgen(s))
Return, [xcm, ycm]
totalAcross2 = Total(array, 2) ; 2 is fastest dim to total across
xcm = Total(Total(totalAcross2, 2) / totalMass * Indgen(s))
ycm = Total(Total(Total(array,1), 2) / totalMass * Indgen(s))
zcm = Total(Total(totalAcross2, 1) / totalMass * Indgen(s))
Return, [xcm, ycm, zcm]
; Time testing was done as follows:
array = findgen(200, 200, 200)
t0 = systime(1)
My timings went from 1.2 seconds (with the previous approach) down to 0.8.
I'd love to see the CASE statement disappear. Who will dare to generalize
this to N dimensions, while ensuring that we total over dimension 2 wherever
Dick Jackson Fanning Software Consulting, Canadian Office
firstname.lastname@example.org Calgary, Alberta Voice/Fax: (403) 242-7398
Coyote's Guide to IDL Programming: http://www.dfanning.com/