How can I efficiently calculate region centroids with non-uniform densities? [duplicate]

Question

Say I have a rectangle of length 100 and width 1. I can very quickly find its centroid with RegionCentroid. But suppose that the farther to the right you are on the rectangle, the denser it is. To find the centroid in this case (imagine finding a point such that you could balance the rectangle on the tip of a pencil at that point), you have to account for the densities. Now I can no longer use builtin RegionCentroid, but the math is quite simple:

Timing[
  length = 100;
  squaretable = Table[Rectangle[{i - 1, 0}], {i, 1, length}];
  bigrectangle = RegionUnion[squaretable];
  squareMF = RegionMember[#] & /@ squaretable;
  dens = Range[length]; 
  densAtPoint[x_?NumericQ, y_?NumericQ] := 
    First[dens[[FirstPosition[squareMF[[#]][{x, y}] & /@ Range[length], True]]]];
  centroid = 
    NIntegrate[{x, y}*densAtPoint[x, y], {x, y} ∈ bigrectangle] /
    NIntegrate[densAtPoint[x, y], {x, y} ∈ bigrectangle]]

 (*{2.65625, {66.5, 0.5}}*)

The problem is that this appears to be highly inefficient. If I increase the length by a factor of ten, the computation time increases by a factor of roughly 100. How can I calculate this centroid more efficiently when I have larger regions (with more subregions with their own densities)?

I'm not looking for an answer that exploits the particular predictability of the density in the MWE that I've given -- please don't tell me to rewrite densAtPoint[x_?NumericQ, y_?NumericQ] := Ceiling[x];!

Answer 1

I don't know if this is helpful, but one idea is to add a dimension representing the density to your regions. Then you can use RegionCentroid of this $n+1$ dimensional object and take the first $n$ coordinates: Here is an example:

cond = Or @@ Table[i-1<x<i && 0<y<1 && 0<z<i, {i, 100}];
Most @ RegionCentroid @ ImplicitRegion[cond, {x,y,z}]

{133/2, 1/2}