8
$\begingroup$

Suppose I have a table with some data that are 'lumped' around some central location, such as the dummy example

i0 = 22.3;
j0 = 34.1;
table = Table[
  Exp[-((i - i0)^2 + (j - j0)^2)/10^2]
  , {i, 1, 50}, {j, 1, 50}
  ]

which kind of looks like

ListPointPlot3D[table, PlotRange -> Full]

Mathematica graphics

I would like to recover the centroid of the "mass" determined by the signal in my table, i.e. in the case above, the tuplet {i0,j0}, in the cleanest way possible. I've managed to do it in the obvious but expensive way with a bunch of explicit sums and products with explicit Ranges, but it feels like there should be a built-in function that will do this ─ and I've as yet been unable to find it.

Can this be done? If so, how?

$\endgroup$

3 Answers 3

12
$\begingroup$
wd = WeightedData[Tuples @ Range @ Dimensions @ table, Join @@ table]
Mean @ wd

{22.3232, 33.9072}

Also

Total[MapIndexed[#2 # &, table / Total[table, 2], {2}], 2] (* and  *)
Dot[Join @@ table , Tuples[Range @ Dimensions @ table]] / Total[table, 2]

{22.3232, 33.9072}

$\endgroup$
9
$\begingroup$
i0 = 22.3; j0 = 34.1;
table = Table[
   Exp[-((i - i0)^2 + (j - j0)^2)/10^2], {i, 1, 50}, {j, 1, 50}];

x = Array[#1 &, Dimensions[table]];
y = Array[#2 &, Dimensions[table]];

pt = {Total@Flatten[x table], Total@Flatten[y table]}/
  Total[table, 2]
(* {22.3232, 33.9072} *)

ListDensityPlot[table, Epilog -> {Red, Point@Reverse[pt]}, 
 PlotRange -> All]

enter image description here

$\endgroup$
1
  • 1
    $\begingroup$ Thanks both. I accepted the 'proper' built-in as that's what the OP asks for, but this is faster for my instance and it is clean enough that it makes for readable code. $\endgroup$ Jul 11, 2018 at 12:19
8
$\begingroup$

I can provide some improvement if it is about speed.

Let's generate a larger data set (I use Compile merely to speed it up a bit).

i0 = 22.3; j0 = 34.1;
m = 2500;
n = 1500;
x = Subdivide[1., 50, m - 1];
y = Subdivide[1., 50, n - 1];
table2 = Partition[#, n] &@
    Compile[{{X, _Real, 1}, {i0, _Real}, {j0, _Real}},
      Exp[-((X[[1]] - i0)^2 + (X[[2]] - j0)^2)/10^2],
      RuntimeAttributes -> {Listable}
      ][Tuples[{x, y}], i0, j0];

Szabolcs' approach

AbsoluteTiming[
 pt = {
    Total@Flatten[Transpose[ConstantArray[x, n]] table],
    Total@Flatten[ConstantArray[y, m] table]
    }/Total[table, 2]
 ]

{0.733425, {22.3288, 33.8733}}

The problem is that before summation, some large arrays have to be constructed and multiplied. But the summations and matrix-matrix products can also be expressed by cheaper matrix-vector products:

AbsoluteTiming[
 pt2 = With[{buffer = ConstantArray[1., m].table},
   {x.(table.ConstantArray[1., n]), 
     buffer.y
   }/(ConstantArray[1., n].buffer)
   ]
 ]

{0.044209, {22.3288, 33.8733}}

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.