I want to calculate the surface of a body made of at least 3 overlapping ellipsoids. Below there is a picture of the cross section of the body.
I already know how to calculate the surface of single ellipsoids.
Does someone know a formula to subtract the overlapping parts or some other kind of formula to calculate the surface?
This problem will not in general have a nice closed-form analytic solution, because it depends on many input parameters (semiaxes and angles of all the ellipses, how many ellipses there are and so on). Using analytic final results would maybe make sense for special cases with less parameters (such as two ellipses at right angles), but for more ellipses the analytic solution won't help you much.
I'll show you how to get all the way to the end and you can also use this very elegant representation and then use numerical integration as the last step.
Consider an axis-alligned ellipse,
$$Ax^2+By^2=1$$ Going into polar coordinates $x=r\cos\phi$, $y=r\sin\phi$, you get $$r^2=\frac{1}{A\cos^2\phi+B\sin^2\phi}$$ You will see later why I didn't even bother with the square root. And of course, this is extremely easy to rotate: $$r^2(\phi)=\frac{1}{A\cos^2(\phi-\alpha)+B\sin^2(\phi-\alpha)}$$
This is now a polar representation of the ellipse centered at the coordinate origin, for any chosen combination of $A,B,\alpha$ (you can express $A$ and $B$ with the half-axes, $A=1/a^2$ and the same for $B$).
Now remember how to get the area inside a curve in the polar coordinates:
$$S=\frac12\int_0^{2\pi} r^2(\phi)\,{\rm d}\phi$$
See how conveniently we only need $r^2$, not $r$. For a single ellipse, you know what to do. For any number of ellipses, it's only one simple step further: the polar representation for the overlapping shape is simply taking the maximum radius at each $\phi$! So all you have to do is
$$r^2(\phi)=\max(r^2(\phi,A_1,B_1,\alpha_1),r^2(\phi,A_2,B_2,\alpha_2),\ldots)$$
In general, this can be done by computer and put into the integral, which can be then computed with any existing numerical integration method. In this case, you can have the $\max$ function used explicitly and you never have to care where the ellipses intersect.
If you want an analytic solution, you have to find for each ellipse, on which range in $\phi$ it has the maximum distance from the origin (you need to find the intersections). Then, you just integrate each ellipse on the interval in $\phi$ where it is the important one for the shape, and add the contributions together.
For instance, two ellipses at right angles, $\alpha_1=0$, $\alpha_2=\frac{\pi}{2}$, intersect at $\phi=\pm \frac{\pi}{4}, \pm \frac{3\pi}{4}$. So in that case, the solution would be
$$S=4\cdot\frac12\int_{-\pi/4}^{\pi/4}\frac{1}{A\cos^2\phi+B\sin^2\phi}d\phi=4ab \arctan\frac{a}{b}$$ where I additionally took into account the fact that all four quarters have the same area.
Even better, the indefinite integral of $r^2$ of a single ellipse is analytical (of course it is, the area of the ellipse is not the problem - the circumference would be but we don't need it):
$$s(\phi)=\int \frac{1}{A\cos^2\phi+B\sin^2 \phi}\,{\rm d}\phi=\frac{1}{\sqrt{AB}}\arctan \left(\sqrt{\frac{B}{A}}\tan \phi\right)=ab\arctan\left(\frac{a}{b}\tan \phi\right)$$
I didn't bother with the one-half factor at the front because every ellipse contributes two parts to the area (symmetrically on both sides of the origin) so it cancels out.
The expression is therefore $$S=\sum_{i=0}^N \left[s(\phi_{i,i+1}-\alpha_i,a_i,b_i)-s_i(\phi_{i,i-1}-\alpha_i,a_i,b_i)\right]$$ where $\phi_{i,i+1}$ is the angle where $i$-th and $i+1$-th ellipses intersect. It's important that the arguments $\phi_{i,i+1}-\alpha_i$ stay in the range $-\pi/2,\pi/2$ otherwise the tangent function will cross a pole and give invalid results.
As you can see, the worst part is finding the intersections.
For the intersections, you need to be careful: some ellipses may not even contribute to the shape (if they are so small they are inside the shape). So this part needs some tests (in programming) or visual inspection (on paper). As you see from your sketch, not all the intersections are the visible outer intersections, so it is really important that you find only the intersections that contribute to the shape.
The intersections can be found analytically, because you just need to solve the equation
$$r_i^2=r_j^2$$ $$A_i\cos^2(\phi-\alpha_i)+B_i\sin^2(\phi-\alpha_i)=A_j\cos^2(\phi-\alpha_j)+B_j\sin^2(\phi-\alpha_j)$$
You can apply the usual trig formulas to split into sines and cosines, and divide by $\cos^2 \phi$. You get a quadratic equation for $\tan\phi$. If it has no solutions, the shapes don't intersect (should not happen if you correctly figured out which intersections you need).
In this calculation I used the unrotated version for $s(\phi)$ and subtracted the rotation angles later. If you "unrotate" the ellipses so that the first one is unrotated when you find the intersection, it's even nices because you get $\tan\phi$ from the quadratic equation for the intersection and that is then exactly what you need in the $s(\phi)$ - no need to go to angles and back to tangents.
I do hope any of this helps. You can decide: do you want to bother with finding the correct intersections and be careful with the angles fitting in the interval between $\pm \pi/2$, or you just let the $\max$ function do everything and you just put the numbers in and get the number out.