Is there a built-in function to get the centroid of a table?

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] 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?

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}

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] • 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. – Emilio Pisanty Jul 11 '18 at 12:19

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[] - i0)^2 + (X[] - 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}}