how to find the geometric centroid of a trapezoid? I'm doing a project involving a trapezoidal prism

I need to find the center of mass. To do this, I need to first find the geometric centroid of the trapezoid.
I have found the formulas online. According to Wolfram Mathworld,
$$\begin{align}
\bar{x} &= \frac{b}{2} + \frac{(2a+b)(c^2-d^2)}{6(b^2-a^2)}\\
\bar{y} &= \frac{b+2a}{3(a+b)}h
\end{align}$$
https://mathworld.wolfram.com/Trapezoid.html
The project I'm doing is too serious for me to just take these formulas as given off the internet. I want to derive them from scratch.
However, I cannot find a derivation.
Can someone derive the geometric centroid of a trapezoid? If it's a long proof, then just sketch the proof for me and I can derive it myself.

If you choose a coordinate system where $b$ is on the positive $x$ axis, between origin $(0, 0)$ and $(b, 0)$, the four vertices of the trapezoid are
$$(0, 0), \quad (b, 0), \quad (f+a, h), \quad (f, h)$$
Since trapezoids are simple polygons, we can use the shoelace formula for the area of a 2D polygon,
$$A = \displaystyle \frac{1}{2} \sum_{i=0}^{n-1} x_i y_{i+1} - x_{i+1} y_i \tag{1}\label{None1}$$
and the centroid of a simple polygon:
$$\left\lbrace ~ \begin{aligned}
\overline{x} &= \displaystyle \frac{1}{6 A} \sum_{i=0}^{n-1} (x_i + x_{i+1})(x_i y_{i+1} - x_{i+1} y_i) \\
\overline{y} &= \displaystyle \frac{1}{6 A} \sum_{i=0}^{n-1} (y_i + y_{i+1})(x_i y_{i+1} - x_{i+1} y_i) \\
\end{aligned} \right . \tag{2}\label{None2}$$
where $x_0 = 0$, $x_1 = b$, $x_2 = f + a$, $x_3 = f$, and $x_4 = x_0 = 0$; and $y_0 = 0$, $y_1 = 0$, $y_2 = h$, $y_3 = 0$, and $y_4 = y_0 = 0$.
These yield
$$\begin{aligned}
\overline{x} &= \displaystyle \frac{f (b + 2 a) + (a + b)^2 - a b}{3 (a + b)} \\
\overline{y} &= \displaystyle h \frac{b + 2 a}{3 (a + b)} \\
\end{aligned} \tag{3}\label{None3}$$
You can solve $f$ and $h$ by splitting the trapezoid into two right triangles and a rectangle in between, and solving the system of three equations and three unknowns (the third being $g$, such that $f + a + g = b$).  It yields
$$\begin{aligned}
f &= \displaystyle \frac{b - a}{2} + \frac{c^2 - d^2}{2 (b - a)} \\
g &= \displaystyle \frac{b - a}{2} - \frac{c^2 - d^2}{2 (b - a)} \\
h &= \displaystyle \frac{\sqrt{ 2 (d^2 + c^2)(b - a)^2 - (d^2 - c^2)^2 - (b - a)^4}}{2 ( b - a ) } \\
\end{aligned} \tag{4}\label{None4}$$
Substituting $f$ in the centroid $x$ coordinate indeed yields, after simplification, $\overline{x} = (b/2) + (2 a + b)(c^2 - d^2) / (6 (b^2 - a^2))$.
As to references, the above links to Wikipedia have references you can use.

Or, you can start with the integral definition.
Given the characteristic function $g(x, y)$ of the shape, $g(x,y)$ being $1$ inside the shape and $0$ outside, the centroid is
$$\begin{aligned}
\overline{x} &= \displaystyle \frac{ \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} x g(x, y) ~ d x ~ d y }{ \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} g(x, y) ~ d x ~ d y } \\
\overline{y} &= \displaystyle \frac{ \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} y g(x, y) ~ d x ~ d y }{ \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} g(x, y) ~ d x ~ d y } \\
\end{aligned} \tag{5}\label{None5}$$
The characteristic function is the weight function, so one might consider this the definition of the 2D centroid.  (For references, look for "centroid integral".)
For the $y$ axis, the integral simplifies to
$$\begin{aligned}
\overline{y} &= \frac{ \int_{0}^{h} y \left(a + \frac{b - a}{h} y\right) ~ dy }{ \int_{0}^{h} a + \frac{b - a}{h} y ~dy } \\
~ &= \frac{ ~ \frac{ 2 b + a }{ 6 } h^2 ~ }{ ~ \frac{ b + a }{2} h ~ } \\
~ &= \frac{ 2 b + a }{ 3 ( b + a ) } h \\
\end{aligned} \tag{6a} \label{None6a}$$
For the $x$ axis, we need to split the integral into three parts. From above, we already know the divisor integral evaluates to $h (b + a) / 2$:
$$\begin{aligned}
\overline{x} &= \frac{2}{h (b + a)} \biggr( \int_{0}^{f} x \left( \frac{h}{f} x \right) ~ d x ~ + \int_{f}^{f+a} x h ~ d x ~ + \int_{f+a}^{b} x \left(\frac{b - x}{b - f - a} h \right) ~ d x \biggr) \\
~ &= \frac{2}{h (b + a)} \biggr( \frac{f^2 h}{3} ~ + ~ \frac{2 a f h + a^2 h}{2} ~ + ~ \frac{ h }{6} \left( f ( b - 4 a - 2 f ) + b^2 + a b - 2 a^2 \right) \biggr) \\
~&= \frac{(2 a + b) f + (a + b)^2 - a b}{3 (a + b)} \\
\end{aligned}$$
which is exactly the same as we got in the polygon form.
