In the book by W.G. Chinn, N.E. Steenrod "First Concepts of Topology" the another remarkable theorem was proved: any two flat bounded regions can be cut by a single straight line so that each of these regions is divided into two regions of equal area (the second  pancake problem). This is an existence theorem which does not provide any way to find this cut. In this post I made an attempt to find such cut for any 2 convex regions on the plane bounded by a piecewise smooth self-non-intersecting curves.
The Each_Into_2_Equal_Areas procedure returns a list of coordinates of 4 endpoints of the cutting segments. This procedure significantly uses my old procedures  Area  and  Picture , which can be found in detail at the link  https://mapleprimes.com/posts/145922-Perimeter-Area-And-Visualization-Of-A-Plane-Figure-  . The formal arguments of the Each_Into_2_Equal_Areas procedure are the lists  L1  and  L2 specifying the boundaries of the regions to be cut. When specifying  L1  and  L2 , the boundary can be passed clockwise or counterclockwise, but it is necessary that the parameter t (when specifying each link) should go in ascending order. This can always be achieved by replacing  t  with  -t  if necessary. The Pic procedure draws a picture of the source regions and cutting segments. For ease of use, the code for the  Area  and  Picture   procedure is also provided. It is also worth noting that the procedure also works for "not too" non-convex regions (see examples below).

restart;
Area := proc(L) 
local i, var, e, e1, e2, P; 
for i to nops(L) do 
if type(L[i], listlist(algebraic)) then 
P[i] := (1/2)*add(L[i, j, 1]*L[i, j+1, 2]-L[i, j, 2]*L[i, j+1, 1], j = 1 .. nops(L[i])-1) else 
var := lhs(L[i, 2]); 
if type(L[i, 1], algebraic) then e := L[i, 1]; 
if nops(L[i]) = 3 then P[i] := (1/2)*(int(e^2, L[i, 2])) else 
if var = y then P[i] := (1/2)*simplify(int(e-var*(diff(e, var)), L[i, 2])) else 
P[i] := (1/2)*simplify(int(var*(diff(e, var))-e, L[i, 2])) end if end if else e1 := L[i, 1, 1]; e2 := L[i, 1, 2]; 
P[i] := (1/2)*simplify(int(e1*(diff(e2, var))-e2*(diff(e1, var)), L[i, 2])) end if end if end do; 
abs(add(P[i], i = 1 .. nops(L))); 
end proc:

Picture := proc(L, C, N::posint := 100, Boundary::list := [linestyle = 1]) 
local i, var, var1, var2, e, e1, e2, P, Q, h; 
global Border;
uses plottools; 
for i to nops(L) do 
if type(L[i], listlist(algebraic)) then P[i] := op(L[i]) else 
var := lhs(L[i, 2]); var1 := lhs(rhs(L[i, 2])); var2 := rhs(rhs(L[i, 2])); h := (var2-var1)/N; 
if type(L[i, 1], algebraic) then e := L[i, 1]; 
if nops(L[i]) = 3 then P[i] := seq(subs(var = var1+h*i, [e*cos(var), e*sin(var)]), i = 0 .. N) else 
P[i] := seq([var1+h*i, subs(var = var1+h*i, e)], i = 0 .. N) fi else e1 := L[i, 1, 1]; e2 := L[i, 1, 2]; P[i] := seq(subs(var = var1+h*i, [e1, e2]), i = 0 .. N) fi; fi; od; 
Q := [seq(P[i], i = 1 .. nops(L))]; Border := curve([op(Q), Q[1]], op(Boundary)); [polygon(Q, C), Border] 
end proc:

Each_Into_2_Equal_Areas:=proc(L1::list, L2::list)
local D, n, m, L10, L20, S1,S2, f, L11, L21, i, j, k, s, A, B, C , sol;

f:=(X,Y)->expand((y-X[2])*(Y[1]-X[1])-(x-X[1])*(Y[2]-X[2]));
L10:=map(p->`if`(type(p,listlist),[[p[1,1]+t*(p[2]-p[1])[1],p[1,2]+t*(p[2]-p[1])[2]],t=0..1],p), L1);
L20:=map(p->`if`(type(p,listlist),[[p[1,1]+t*(p[2]-p[1])[1],p[1,2]+t*(p[2]-p[1])[2]],t=0..1],p), L2);
S1:=Area(L1); S2:=Area(L2);  
n:=nops(L1); m:=nops(L2);

for i from 1 to n do
for j from i to n do

for k from 1 to m do
for s from k to m do

if not ((nops({i,j})=1 and type(L1[i],listlist)) or (nops({k,s})=1 and type(L2[k],listlist))) then

A:=eval(L10[i,1],t=t1): 
B:=eval(L10[j,1],t=t2):
C:=eval(L20[k,1],t=t3): 
D:=eval(L20[s,1],t=t4):

L11:=`if`(j=i,[subsop([2,2]=t1..t2,L10[i]),[B,A]],`if`(j=i+1,[subsop([2,2]=t1..op([2,2,2],L10[i]),L10[i]),subsop([2,2]=op([2,2,1],L10[j])..t2,L10[j]),[B,A]], [subsop([2,2]=t1..op([2,2,2],L10[i]),L10[i]),op(L10[i+1..j-1]),subsop([2,2]=op([2,2,1],L10[j])..t2,L10[j]),[B,A]])):

L21:=`if`(s=k,[subsop([2,2]=t3..t4,L20[k]),[D,C]],`if`(s=k+1,[subsop([2,2]=t3..op([2,2,2],L20[k]),L20[k]),subsop([2,2]=op([2,2,1],L20[s])..t4,L20[s]),[D,C]], [subsop([2,2]=t3..op([2,2,2],L20[k]),L20[k]),op(L20[k+1..s-1]),subsop([2,2]=op([2,2,1],L20[s])..t4,L20[s]),[D,C]])):

sol:=fsolve(simplify({Area(L11)-S1/2,Area(L21)-S2/2,eval(f(A,B),[x=C[1],y=C[2]]),eval(f(A,B),[x=D[1],y=D[2]])}),{t1=op([2,2,1],L10[i])..op([2,2,2],L10[i]),t2=op([2,2,1],L10[j])..op([2,2,2],L10[j]),t3=op([2,2,1],L20[k])..op([2,2,2],L20[k]),t4=op([2,2,1],L20[s])..op([2,2,2],L20[s])});

if type(sol,set(`=`)) then  return eval([A,B,C,D],sol) fi;

fi;
od: od: od: od:
end proc:

Pic := proc(L1,L2,col1,col2,Size:=[800,400])
local P1, P2, P3, T, P;
uses plots, plottools;
P1, P2 := Picture(L1, color=col1), Picture(L2, color=col2):
P3 := line(Sol[1..2][],color=red,thickness=3), line(Sol[3..4][],color=red,thickness=3), line(Sol[1],Sol[4],linestyle=2,thickness=3,color=red):
T:=textplot([[Sol[1][],"A"],[Sol[2][],"B"],[Sol[3][],"C"],[Sol[4][],"D"]], font=[times,18], align=[left,above]);
P:=pointplot(Sol, symbol=solidcircle, color=red, symbolsize=10);
display(P1,P2,P3,T,P, scaling=constrained, size=Size, axes=none);
end proc: 

Examples of use.

local D:
L1:=[[[8,0],[6,7]],[[6,7],[2,5]],[[2,5],[0,2]],[[0,2],[0,0]],[[0,0],[8,0]]]:
L2:=[[[5*cos(t)+16,5*sin(t)],t=Pi/2..Pi],[[5*cos(t)+16,5*sin(t)/2],t=Pi..2*Pi],[[21,0],[16,5]]]:
Sol:=Each_Into_2_Equal_Areas(L1,L2): Points:=[A,B,C,D]:
seq(Points[i]=Sol[i], i=1..4);
Pic(L1,L2,"Yellow","LightBlue",[900,400]);

The specified regions may overlap:

L1:=[[[8,0],[6,7]],[[6,7],[2,5]],[[2,5],[0,2]],[[0,2],[0,0]],[[0,0],[8,0]]]:
L2:=[[[5*cos(t)+9,5*sin(t)],t=Pi/2..Pi],[[5*cos(t)+9,5*sin(t)/2],t=Pi..2*Pi],[[14,0],[9,5]]]:
Sol:=Each_Into_2_Equal_Areas(L1,L2):  Points:=[A,B,C,D]:
seq(Points[i]=Sol[i], i=1..4);
Pic(L1,L2,"Yellow","LightBlue");

If there is a solution for which the cutting segments intersect the boundary of each of the regions at 2 points, then the procedure also works for such non-convex regions:

L1:=[[[cos(t),sin(t)],t=Pi/3..2*Pi-Pi/3],[[cos(-t)+1,sin(-t)],t=-Pi-Pi/3..-Pi+Pi/3]]:
L2:=[[[cos(t)+2,sin(t)],t=-Pi/6..Pi+Pi/6],[[cos(-t)+2,sin(-t)-1],t=-5*Pi/6..-Pi/6]]:
Sol:=Each_Into_2_Equal_Areas(L1,L2): Points:=[A,B,C,D]:
seq(Points[i]=Sol[i], i=1..4);
Pic(L1,L2,"Yellow","LightBlue");

A number of other interesting examples can be found in the attached file.

Each_Into_2_Equal_Areas1.mw