Double integral computation with triangular function

Question

Suppose we have $6$ constants, $(x_l, x_m, x_u$, $y_l, y_m, y_u)\in \mathbb R^6$ such as:

$$x_l < x_m < x_u$$

$$y_l < y_m < y_u$$

Let $D = \left\{(x, y) : x_l\le x\le x_u, \quad y_l \le y \le y_u, \quad y < x\right\}$

I want to compute:

$$\displaystyle \iint_D (x-y)f(x, x_l, x_m,x_u)f(y,y_l,y_m,y_u)dxdy$$

Where $(a,b,c) \in \mathbb R^3, a<b<c$ and

$\displaystyle f(x, a, b, c) = \begin{cases}\dfrac{2(x-a)}{(c-a)(b-a)} &\text{ if } a \le x < b \\ \dfrac{2(x-c)}{(c-a)(b-c)} &\text{ if }b < x \le c\\ \dfrac2{(c-a)} &\text{ if } x = b\\ 0& \text{ otherwise}\end{cases} $

The following picture is an example of $f$ with $x$ in abscissa where $a = 50, b = 55, c = 60, \dfrac{2}{c - a} = 0.2$. Note that $c-b$ can be different from $b-a$.

Actually, the area under $f$ is always equal to $1$ no matter its arguments (it's a distribution law, but that's not necessary to be known for solving the integral).

I have started to try computing this double integral by separating in several cases:

• When $y_l \ge x_u$, this integral is $0$

• When $x_m \le y_l \le x_u$ and $y_m \ge x_u$

• When $x_l \le y_l \le x_m$ and $y_m \ge x_u$

• When...

And so on until all values of the $6$ constants are matched. But that is particularly inefficient and long as it requires the computation of the double integrals a lot of times (see EDIT 1). Is there a quickest way?

EDIT 1: Actually, lot of times here means 20 times as there are 20 possible combinations for the 6 constants.

EDIT 2: if you do tell me there is no quickest way, I will compute the 20 integrals thanks to this related question I posted before this one.

EDIT 3: I tried to write a code in SAGEMath to have a clue about what to find. The final result is quite expanded:

var('x,y,xl,xm,xu,yl,ym,yu')
assume(yl<ym)
assume(ym<yu)
assume(xl<xm)
assume(xm<xu)
assume(xl<x)
assume(x<xu)
assume(yl<y)
assume(y<yu)
C1x = 2/((xu-xl)*(xm-xl))
C2x = 2/((xu-xl)*(xm-xu))
C1y = 2/((yu-yl)*(ym-yl))
C2y = 2/((yu-yl)*(ym-yu))
fx(x,a,b,c) = unit_step(x-a) * unit_step(b-x) * (C1x * x - a*C1x) + unit_step(x - b)*unit_step(c - x)* (C2x* x - C2x*c)
fy(y,a,b,c) = unit_step(y-a) * unit_step(b-y) * (C1y * y - a*C1y) + unit_step(y - b)*unit_step(c - y)* (C2y* y - C2y*c)
expr1 = integrate((x-y)*fx(x,xl,xm,xu)*fy(y,yl,ym,yu), x, max_symbolic(xl, y),xu)
expr2 = integrate(expr1, y, yl, yu)

SAGEMaths returns the following result which is, I guess, correct, if we assume my code is. However that's not a proof but that would be ok if we are certain it's the correct result. What's more annoying is, this is a really complicated form. Just converting it in correct LaTeX would take a huge amount of time for instance.

1/120*((xl - xm)*xu^5*yl*sgn(-xu + yl) - ((xm - xu)*yl*sgn(-xl + yl) - (xl - xu)*yl*sgn(-xm + yl) + (xl - xm)*yl*sgn(-xu + yl))*ym^5 + 5*((xl - xm)*xu*yl*sgn(-xu + yl) - 8*xl^2*xm + 8*xl*xm^2 - 8*(xl - xm)*xu^2 + (xl*xm - xl*xu)*yl*sgn(-xl + yl) - (xl*xm - xm*xu)*yl*sgn(-xm + yl) + 8*(xl^2 - xm^2)*xu)*ym^4 - 10*((xl - xm)*xu^2*yl*sgn(-xu + yl) - 2*xl^3*xm + 2*xl*xm^3 - 2*(xl - xm)*xu^3 + (xl^2*xm - xl^2*xu)*yl*sgn(-xl + yl) - (xl*xm^2 - xm^2*xu)*yl*sgn(-xm + yl) + 2*(xl^3 - xm^3)*xu - 6*(xl^2*xm - xl*xm^2 + (xl - xm)*xu^2 - (xl^2 - xm^2)*xu)*yl)*ym^3 + 10*((xl - xm)*xu^3*yl*sgn(-xu + yl) + (xl^3*xm - xl^3*xu)*yl*sgn(-xl + yl) - (xl*xm^3 - xm^3*xu)*yl*sgn(-xm + yl) - 4*(xl^3*xm - xl*xm^3 + (xl - xm)*xu^3 - (xl^3 - xm^3)*xu)*yl)*ym^2 + (xl^5*xm - xl^5*xu)*yl*sgn(-xl + yl) - (xl*xm^5 - xm^5*xu)*yl*sgn(-xm + yl) - (20*(xl^2*xm - xl*xm^2 + (xl - xm)*xu^2 - (xl^2 - xm^2)*xu)*yl^3 - 20*(xl^3*xm - xl*xm^3 + (xl - xm)*xu^3 - (xl^3 - xm^3)*xu)*yl^2 + (2*xl^5*xm - 2*xl^5*xu - 2*(xm - xu)*yl^5 + 10*(xl*xm - xl*xu)*yl^4 - 20*(xl^2*xm - xl^2*xu)*yl^3 + 20*(xl^3*xm - xl^3*xu)*yl^2 - 5*(xl^4*xm - xl^4*xu)*yl)*sgn(-xl + yl) - (2*xl*xm^5 - 2*xm^5*xu - 2*(xl - xu)*yl^5 + 10*(xl*xm - xm*xu)*yl^4 - 20*(xl*xm^2 - xm^2*xu)*yl^3 + 20*(xl*xm^3 - xm^3*xu)*yl^2 - 5*(xl*xm^4 - xm^4*xu)*yl)*sgn(-xm + yl) + (2*(xl - xm)*xu^5 - 5*(xl - xm)*xu^4*yl + 20*(xl - xm)*xu^3*yl^2 - 20*(xl - xm)*xu^2*yl^3 + 10*(xl - xm)*xu*yl^4 - 2*(xl - xm)*yl^5)*sgn(-xu + yl))*ym + (((xm - xu)*sgn(-xl + yl) - (xl - xu)*sgn(-xm + yl) + (xl - xm)*sgn(-xu + yl))*ym^5 - 5*((xl - xm)*xu*sgn(-xu + yl) + (xl*xm - xl*xu)*sgn(-xl + yl) - (xl*xm - xm*xu)*sgn(-xm + yl))*ym^4 + 20*(xl^2*xm - xl*xm^2 + (xl - xm)*xu^2 - (xl^2 - xm^2)*xu)*yl^3 + 10*((xl - xm)*xu^2*sgn(-xu + yl) + 4*xl^2*xm - 4*xl*xm^2 + 4*(xl - xm)*xu^2 - 4*(xl^2 - xm^2)*xu + (xl^2*xm - xl^2*xu)*sgn(-xl + yl) - (xl*xm^2 - xm^2*xu)*sgn(-xm + yl))*ym^3 - 20*(xl^3*xm - xl*xm^3 + (xl - xm)*xu^3 - (xl^3 - xm^3)*xu)*yl^2 - 10*((xl - xm)*xu^3*sgn(-xu + yl) + 2*xl^3*xm - 2*xl*xm^3 + 2*(xl - xm)*xu^3 - 2*(xl^3 - xm^3)*xu + 6*(xl^2*xm - xl*xm^2 + (xl - xm)*xu^2 - (xl^2 - xm^2)*xu)*yl + (xl^3*xm - xl^3*xu)*sgn(-xl + yl) - (xl*xm^3 - xm^3*xu)*sgn(-xm + yl))*ym^2 + 5*((xl - xm)*xu^4*sgn(-xu + yl) + 8*(xl^3*xm - xl*xm^3 + (xl - xm)*xu^3 - (xl^3 - xm^3)*xu)*yl + (xl^4*xm - xl^4*xu)*sgn(-xl + yl) - (xl*xm^4 - xm^4*xu)*sgn(-xm + yl))*ym + (xl^5*xm - xl^5*xu - 2*(xm - xu)*yl^5 + 10*(xl*xm - xl*xu)*yl^4 - 20*(xl^2*xm - xl^2*xu)*yl^3 + 20*(xl^3*xm - xl^3*xu)*yl^2 - 10*(xl^4*xm - xl^4*xu)*yl)*sgn(-xl + yl) - ((xm - xu)*ym^5*sgn(-xl + yl) - 5*(xl*xm - xl*xu)*ym^4*sgn(-xl + yl) + 10*(xl^2*xm - xl^2*xu)*ym^3*sgn(-xl + yl) - 10*(xl^3*xm - xl^3*xu)*ym^2*sgn(-xl + yl) + 5*(xl^4*xm - xl^4*xu)*ym*sgn(-xl + yl) - (xl^5*xm - xl^5*xu)*sgn(-xl + yl))*sgn(-xl + ym) - (xl*xm^5 - xm^5*xu - 2*(xl - xu)*yl^5 + 10*(xl*xm - xm*xu)*yl^4 - 20*(xl*xm^2 - xm^2*xu)*yl^3 + 20*(xl*xm^3 - xm^3*xu)*yl^2 - 10*(xl*xm^4 - xm^4*xu)*yl)*sgn(-xm + yl) + ((xl - xu)*ym^5*sgn(-xm + yl) - 5*(xl*xm - xm*xu)*ym^4*sgn(-xm + yl) + 10*(xl*xm^2 - xm^2*xu)*ym^3*sgn(-xm + yl) - 10*(xl*xm^3 - xm^3*xu)*ym^2*sgn(-xm + yl) + 5*(xl*xm^4 - xm^4*xu)*ym*sgn(-xm + yl) - (xl*xm^5 - xm^5*xu)*sgn(-xm + yl))*sgn(-xm + ym) + ((xl - xm)*xu^5 - 10*(xl - xm)*xu^4*yl + 20*(xl - xm)*xu^3*yl^2 - 20*(xl - xm)*xu^2*yl^3 + 10*(xl - xm)*xu*yl^4 - 2*(xl - xm)*yl^5)*sgn(-xu + yl) + ((xl - xm)*xu^5*sgn(-xu + yl) - 5*(xl - xm)*xu^4*ym*sgn(-xu + yl) + 10*(xl - xm)*xu^3*ym^2*sgn(-xu + yl) - 10*(xl - xm)*xu^2*ym^3*sgn(-xu + yl) + 5*(xl - xm)*xu*ym^4*sgn(-xu + yl) - (xl - xm)*ym^5*sgn(-xu + yl))*sgn(-xu + ym))*yu + ((xm - xu)*yl*ym^5*sgn(-xl + yl) - 5*(xl*xm - xl*xu)*yl*ym^4*sgn(-xl + yl) + 10*(xl^2*xm - xl^2*xu)*yl*ym^3*sgn(-xl + yl) - 10*(xl^3*xm - xl^3*xu)*yl*ym^2*sgn(-xl + yl) + 5*(xl^4*xm - xl^4*xu)*yl*ym*sgn(-xl + yl) - (xl^5*xm - xl^5*xu)*yl*sgn(-xl + yl))*sgn(-xl + ym) - ((xl - xu)*yl*ym^5*sgn(-xm + yl) - 5*(xl*xm - xm*xu)*yl*ym^4*sgn(-xm + yl) + 10*(xl*xm^2 - xm^2*xu)*yl*ym^3*sgn(-xm + yl) - 10*(xl*xm^3 - xm^3*xu)*yl*ym^2*sgn(-xm + yl) + 5*(xl*xm^4 - xm^4*xu)*yl*ym*sgn(-xm + yl) - (xl*xm^5 - xm^5*xu)*yl*sgn(-xm + yl))*sgn(-xm + ym) - ((xl - xm)*xu^5*yl*sgn(-xu + yl) - 5*(xl - xm)*xu^4*yl*ym*sgn(-xu + yl) + 10*(xl - xm)*xu^3*yl*ym^2*sgn(-xu + yl) - 10*(xl - xm)*xu^2*yl*ym^3*sgn(-xu + yl) + 5*(xl - xm)*xu*yl*ym^4*sgn(-xu + yl) - (xl - xm)*yl*ym^5*sgn(-xu + yl))*sgn(-xu + ym))/((xl^2*xm - xl*xm^2 + (xl - xm)*xu^2 - (xl^2 - xm^2)*xu)*yl^2*ym - (xl^2*xm - xl*xm^2 + (xl - xm)*xu^2 - (xl^2 - xm^2)*xu)*yl*ym^2 + ((xl^2*xm - xl*xm^2 + (xl - xm)*xu^2 - (xl^2 - xm^2)*xu)*yl - (xl^2*xm - xl*xm^2 + (xl - xm)*xu^2 - (xl^2 - xm^2)*xu)*ym)*yu^2 - ((xl^2*xm - xl*xm^2 + (xl - xm)*xu^2 - (xl^2 - xm^2)*xu)*yl^2 - (xl^2*xm - xl*xm^2 + (xl - xm)*xu^2 - (xl^2 - xm^2)*xu)*ym^2)*yu) + 1/120*(40*(xl^2*xm - xl*xm^2 + (xl - xm)*xu^2 - (xl^2 - xm^2)*xu)*ym^4 + 20*(xl^3*xm - xl*xm^3 + (xl - xm)*xu^3 - (xl^3 - xm^3)*xu)*yl*ym^2 - 20*(xl^3*xm - xl*xm^3 + (xl - xm)*xu^3 - (xl^3 - xm^3)*xu + 2*(xl^2*xm - xl*xm^2 + (xl - xm)*xu^2 - (xl^2 - xm^2)*xu)*yl)*ym^3 - 20*((xl^2*xm - xl*xm^2 + (xl - xm)*xu^2 - (xl^2 - xm^2)*xu)*yl - (xl^2*xm - xl*xm^2 + (xl - xm)*xu^2 - (xl^2 - xm^2)*xu)*ym)*yu^3 + 20*((xl^3*xm - xl*xm^3 + (xl - xm)*xu^3 - (xl^3 - xm^3)*xu)*yl - (xl^3*xm - xl*xm^3 + (xl - xm)*xu^3 - (xl^3 - xm^3)*xu)*ym)*yu^2 - 20*(3*(xl^2*xm - xl*xm^2 + (xl - xm)*xu^2 - (xl^2 - xm^2)*xu)*ym^3 + 2*(xl^3*xm - xl*xm^3 + (xl - xm)*xu^3 - (xl^3 - xm^3)*xu)*yl*ym - (2*xl^3*xm - 2*xl*xm^3 + 2*(xl - xm)*xu^3 - 2*(xl^3 - xm^3)*xu + 3*(xl^2*xm - xl*xm^2 + (xl - xm)*xu^2 - (xl^2 - xm^2)*xu)*yl)*ym^2)*yu - ((xm - xu)*yl*ym^5 - 5*(xl*xm - xl*xu)*yl*ym^4 - 2*((xm - xu)*yl - (xm - xu)*ym)*yu^5 + 10*(xl^2*xm - xl^2*xu)*yl*ym^3 + 10*((xl*xm - xl*xu)*yl - (xl*xm - xl*xu)*ym)*yu^4 - 10*(xl^3*xm - xl^3*xu)*yl*ym^2 - 20*((xl^2*xm - xl^2*xu)*yl - (xl^2*xm - xl^2*xu)*ym)*yu^3 + 20*((xl^3*xm - xl^3*xu)*yl - (xl^3*xm - xl^3*xu)*ym)*yu^2 + (xl^5*xm - xl^5*xu)*yl - (2*xl^5*xm - 2*xl^5*xu - 5*(xl^4*xm - xl^4*xu)*yl)*ym + (xl^5*xm - xl^5*xu - (xm - xu)*ym^5 + 5*(xl*xm - xl*xu)*ym^4 - 10*(xl^2*xm - xl^2*xu)*ym^3 + 10*(xl^3*xm - xl^3*xu)*ym^2 - 10*(xl^4*xm - xl^4*xu)*yl + 5*(xl^4*xm - xl^4*xu)*ym + (xl^5*xm - xl^5*xu - (xm - xu)*ym^5 + 5*(xl*xm - xl*xu)*ym^4 - 10*(xl^2*xm - xl^2*xu)*ym^3 + 10*(xl^3*xm - xl^3*xu)*ym^2 - 5*(xl^4*xm - xl^4*xu)*ym)*sgn(-xl + ym))*yu + ((xm - xu)*yl*ym^5 - 5*(xl*xm - xl*xu)*yl*ym^4 + 10*(xl^2*xm - xl^2*xu)*yl*ym^3 - 10*(xl^3*xm - xl^3*xu)*yl*ym^2 + 5*(xl^4*xm - xl^4*xu)*yl*ym - (xl^5*xm - xl^5*xu)*yl)*sgn(-xl + ym))*sgn(-xl + yu) + ((xl - xu)*yl*ym^5 - 5*(xl*xm - xm*xu)*yl*ym^4 - 2*((xl - xu)*yl - (xl - xu)*ym)*yu^5 + 10*(xl*xm^2 - xm^2*xu)*yl*ym^3 + 10*((xl*xm - xm*xu)*yl - (xl*xm - xm*xu)*ym)*yu^4 - 10*(xl*xm^3 - xm^3*xu)*yl*ym^2 - 20*((xl*xm^2 - xm^2*xu)*yl - (xl*xm^2 - xm^2*xu)*ym)*yu^3 + 20*((xl*xm^3 - xm^3*xu)*yl - (xl*xm^3 - xm^3*xu)*ym)*yu^2 + (xl*xm^5 - xm^5*xu)*yl - (2*xl*xm^5 - 2*xm^5*xu - 5*(xl*xm^4 - xm^4*xu)*yl)*ym + (xl*xm^5 - xm^5*xu - (xl - xu)*ym^5 + 5*(xl*xm - xm*xu)*ym^4 - 10*(xl*xm^2 - xm^2*xu)*ym^3 + 10*(xl*xm^3 - xm^3*xu)*ym^2 - 10*(xl*xm^4 - xm^4*xu)*yl + 5*(xl*xm^4 - xm^4*xu)*ym + (xl*xm^5 - xm^5*xu - (xl - xu)*ym^5 + 5*(xl*xm - xm*xu)*ym^4 - 10*(xl*xm^2 - xm^2*xu)*ym^3 + 10*(xl*xm^3 - xm^3*xu)*ym^2 - 5*(xl*xm^4 - xm^4*xu)*ym)*sgn(-xm + ym))*yu + ((xl - xu)*yl*ym^5 - 5*(xl*xm - xm*xu)*yl*ym^4 + 10*(xl*xm^2 - xm^2*xu)*yl*ym^3 - 10*(xl*xm^3 - xm^3*xu)*yl*ym^2 + 5*(xl*xm^4 - xm^4*xu)*yl*ym - (xl*xm^5 - xm^5*xu)*yl)*sgn(-xm + ym))*sgn(-xm + yu) - ((xl - xm)*xu^5*yl - 10*(xl - xm)*xu^3*yl*ym^2 + 10*(xl - xm)*xu^2*yl*ym^3 - 5*(xl - xm)*xu*yl*ym^4 + (xl - xm)*yl*ym^5 - 2*((xl - xm)*yl - (xl - xm)*ym)*yu^5 + 10*((xl - xm)*xu*yl - (xl - xm)*xu*ym)*yu^4 - 20*((xl - xm)*xu^2*yl - (xl - xm)*xu^2*ym)*yu^3 + 20*((xl - xm)*xu^3*yl - (xl - xm)*xu^3*ym)*yu^2 - (2*(xl - xm)*xu^5 - 5*(xl - xm)*xu^4*yl)*ym + ((xl - xm)*xu^5 - 10*(xl - xm)*xu^4*yl + 5*(xl - xm)*xu^4*ym + 10*(xl - xm)*xu^3*ym^2 - 10*(xl - xm)*xu^2*ym^3 + 5*(xl - xm)*xu*ym^4 - (xl - xm)*ym^5 + ((xl - xm)*xu^5 - 5*(xl - xm)*xu^4*ym + 10*(xl - xm)*xu^3*ym^2 - 10*(xl - xm)*xu^2*ym^3 + 5*(xl - xm)*xu*ym^4 - (xl - xm)*ym^5)*sgn(-xu + ym))*yu - ((xl - xm)*xu^5*yl - 5*(xl - xm)*xu^4*yl*ym + 10*(xl - xm)*xu^3*yl*ym^2 - 10*(xl - xm)*xu^2*yl*ym^3 + 5*(xl - xm)*xu*yl*ym^4 - (xl - xm)*yl*ym^5)*sgn(-xu + ym))*sgn(-xu + yu))/((xl^2*xm - xl*xm^2 + (xl - xm)*xu^2 - (xl^2 - xm^2)*xu)*yl^2*ym - (xl^2*xm - xl*xm^2 + (xl - xm)*xu^2 - (xl^2 - xm^2)*xu)*yl*ym^2 + ((xl^2*xm - xl*xm^2 + (xl - xm)*xu^2 - (xl^2 - xm^2)*xu)*yl - (xl^2*xm - xl*xm^2 + (xl - xm)*xu^2 - (xl^2 - xm^2)*xu)*ym)*yu^2 - ((xl^2*xm - xl*xm^2 + (xl - xm)*xu^2 - (xl^2 - xm^2)*xu)*yl^2 - (xl^2*xm - xl*xm^2 + (xl - xm)*xu^2 - (xl^2 - xm^2)*xu)*ym^2)*yu)

Answer 1

Let me first make a little change to the parameters and put them according to the sketch below.

Then let me introduce the Iverson bracket $$ \left[ P \right] = \left\{ {\begin{array}{*{20}c} 1 & {P = TRUE} \\ 0 & {P = FALSE} \\ \end{array} } \right. $$ which in this type of problems is of profitable use, and which is easy to program.

So, the function $f(x)$ can be written as $$ \eqalign{ & f(x) = {2 \over {\left( {b + c} \right)}}\left( {{{x - a} \over b}\;\left[ {a \le x < a + b} \right] - \left( {1 - {{x - a - b} \over c}} \right)\;\left[ {a + b \le x < a + b + c} \right]} \right) = \cr & = {2 \over {\left( {b + c} \right)}}\left( {{{x - a} \over b}\;\left[ {a \le x < a + b} \right] - \left( {{{a + b + c - x} \over c}} \right)\;\left[ {a + b \le x < a + b + c} \right]} \right) \cr} $$

The integrand function is $$ \eqalign{ & g(x,y) = \left( {x - y} \right)\left[ {y < x} \right]f(x)f(y) = \cr & = {4 \over {\left( {b_{\,x} + c_{\,x} } \right)\left( {b_{\,y} + c_{\,y} } \right)}}\left( {x - y} \right)\left[ {y < x} \right]\; \cdot \cr & \cdot \;\left( {{{\left( {x - a_{\,x} } \right)} \over {b_{\,x} }}\;\left[ {a_{\,x} \le x < a_{\,x} + b_{\,x} } \right] - \left( {{{a_{\,x} + b_{\,x} + c_{\,x} - x} \over {c_{\,x} }}} \right)\;\left[ {a_{\,x} + b_{\,x} \le x < a_{\,x} + b_{\,x} + c_{\,x} } \right]} \right)\; \cdot \cr & \cdot \;\left( {{{\left( {y - a_{\,y} } \right)} \over {b_{\,y} }}\;\left[ {a_{\,y} \le y < a_{\,y} + b_{\,y} } \right] - \left( {{{a_{\,y} + b_{\,y} + c_{\,y} - y} \over {c_{\,y} }}} \right)\;\left[ {a_{\,y} + b_{\,y} \le y < a_{\,y} + b_{\,y} + c_{\,y} } \right]} \right) \cr} $$

Now, the function $g(x,y)$ is intrinsically null outside of the established domain of integration, so that we can write the integral as $$ \eqalign{ & I = \int_{x = - \infty }^\infty {\int_{y = - \infty }^\infty {g(x,y)dxdy} } = \cr & = \int\!\!\!\int\limits_{\left( {x,y} \right)\, \in \,\left[ {\min \left( {a_{\,x} ,a_{\,y} } \right)\,,\;\max \left( {a_{\,x} + b_{\,x} + c_{\,x} ,a_{\,y} + b_{\,y} + c_{\,y} } \right)} \right]^{\,2} } {g(x,y)dxdy} = \cr & = \int_{x = a_{\,x} }^{a_{\,x} + b_{\,x} + c_{\,x} } {dx\int_{y = a_{\,y} }^{a_{\,y} + b_{\,y} + c_{\,y} } {g(x,y)dy} } \cr} $$

Analyzing further the integral, we see that it splits into the sum of four integrals $$ \eqalign{ & I_{\,n} = I\;\mathop /\limits_{} \;{4 \over {\left( {b_{\,x} + c_{\,x} } \right)\left( {b_{\,y} + c_{\,y} } \right)}} = \int_{x = - \infty }^\infty {\int_{y = - \infty }^\infty {g(x,y)dxdy} } = \int\!\!\!\int\limits_{\left( {x,y} \right)\, \in \,U^{\,2} } {g(x,y)dxdy} = \cr & = \int\!\!\!\int\limits_{\left( {x,y} \right)\, \in \,U^{\,2} } {h_{\,1} (x,y)\left[ {c_{\,1} (x,y)} \right]dxdy} + \cdots + \int\!\!\!\int\limits_{\left( {x,y} \right)\, \in \,U^{\,2} } {h_{\,4} (x,y)\left[ {c_{\,4} (x,y)} \right]dxdy} = \cr & = \int\!\!\!\int\limits_{\left( {x,y} \right)\;:\;c\,_1 (x,y) = TRUE} {h_1 (x,y)\,dxdy} + \cdots + \int\!\!\!\int\limits_{\left( {x,y} \right)\;:\;c_{\,4} (x,y) = TRUE} {h_4 (x,y)\,dxdy} \cr} $$ where $$ U = \left[ {\min \left( {a_{\,x} ,a_{\,y} } \right)\,,\;\max \left( {a_{\,x} + b_{\,x} + c_{\,x} ,a_{\,y} + b_{\,y} + c_{\,y} } \right)} \right] $$ which has mainly a computational utility, in avoiding the infinite set,
and where $$ \eqalign{ & h_{\,1} (x,y) = \left( {x - y} \right){{\left( {x - a_{\,x} } \right)} \over {b_{\,x} }}{{\left( {y - a_{\,y} } \right)} \over {b_{\,y} }} \cr & \left[ {c_{\,1} (x,y)} \right] = \left[ {y < x} \right]\left[ {a_{\,x} \le x < a_{\,x} + b_{\,x} } \right]\left[ {a_{\,y} \le y < a_{\,y} + b_{\,y} } \right] \cr} $$ the other clearly following from the product of the terms in $g(x,y)$.

$\left[ {c_{\,1} (x,y)} \right] $ has the role of the charactestic function of the integration interval, and can be represented in various alternative ways $$ \eqalign{ & \left[ {c_{\,1} (x,y)} \right] = \left[ {y < x} \right]\left[ {a_{\,x} \le x < a_{\,x} + b_{\,x} } \right]\left[ {a_{\,y} \le y < a_{\,y} + b_{\,y} } \right] = \cr & = \left[ {\left( {y < x} \right) \wedge \left( {a_{\,x} \le x < a_{\,x} + b_{\,x} } \right) \wedge \left( {a_{\,y} \le y < a_{\,y} + b_{\,y} } \right)} \right] \cr & c_{\,1} (x,y) = TRUE\quad \Leftrightarrow \quad \left( {x,y} \right):\;\left\{ \matrix{ y < x \hfill \cr a_{\,x} \le x < a_{\,x} + b_{\,x} \hfill \cr a_{\,y} \le y < a_{\,y} + b_{\,y} \hfill \cr} \right. \cr} $$

That means that you can compute each integral within the same bounds $x \in U, y \in U$ and keeping $c_k(x,y)$ in the integrand,
or solve the set of inequalities for each and obtain therefrom the integration bounds.

Example

$$ \left\{ \matrix{ a_{\,x} = 0,\,b_{\,x} = 2,\,c_{\,x} = 3 \hfill \cr a_{\,y} = 1,\,b_{\,y} = 2,\,c_{\,y} = 2 \hfill \cr} \right.\quad \Rightarrow \quad U = \left[ {0,5} \right] $$ then $$ \eqalign{ & h_{\,1} (x,y) = {1 \over 4}x\left( {x - y} \right)\left( {y - 1} \right) \cr & \left[ {c_{\,1} (x,y)} \right] = \left[ {y < x} \right]\left[ {0 \le x < 2} \right]\left[ {1 \le y < 3} \right] = \left[ {1 \le y < x} \right]\left[ {1 < x < 2} \right] \cr & h_{\,2} (x,y) = {1 \over 4}x\left( {x - y} \right)\left( {y - 5} \right) \cr & \left[ {c_{\,2} (x,y)} \right] = \left[ {y < x} \right]\left[ {0 \le x < 2} \right]\left[ {3 \le y < 5} \right] = 0 \cr & h_{\,3} (x,y) = {1 \over 6}\left( {x - y} \right)\left( {x - 5} \right)\left( {y - 1} \right) \cr & \left[ {c_{\,3} (x,y)} \right] = \left[ {y < x} \right]\left[ {2 \le x < 5} \right]\left[ {1 \le y < 3} \right] = \cr & = \left[ {y < x} \right]\left[ {2 \le x < 5} \right]\left( {\left[ {1 \le y < 2} \right] + \left[ {2 \le y < 3} \right]} \right) = \cr & = \left[ {2 \le x < 5} \right]\left[ {1 \le y < 2} \right] + \left[ {2 \le y < 3} \right]\left[ {y \le x < 5} \right] \cr & h_{\,4} (x,y) = {1 \over 6}\left( {x - y} \right)\left( {5 - x} \right)\left( {5 - y} \right) \cr & \left[ {c_{\,4} (x,y)} \right] = \left[ {y < x} \right]\left[ {2 \le x < 5} \right]\left[ {3 \le y < 5} \right] = \left[ {3 \le y < x} \right]\left[ {3 \le x < 5} \right] \cr} $$

therefore, the first integral is $$ \eqalign{ & I_{\,n\,1} = \int_{x = 0}^{\,5} {dx\int_{y = 0}^{\,5} {h_{\,1} (x,y)\left[ {c_{\,1} (x,y)} \right]dy} } = \cr & = {1 \over 4}\int_{x = 0}^{\,5} {x\left[ {1 < x < 2} \right]dx\int_{y = 0}^{\,5} {\left( {x - y} \right)\left( {y - 1} \right)\left[ {1 \le y < x} \right]dy} } = \cr & = {1 \over 4}\int_{x = 1}^{\,2} {xdx\int_{y = 1}^{\,x} {\left( {x - y} \right)\left( {y - 1} \right)dy} } = \cr & = {1 \over 4}\int_{u = 0}^{\,1} {udx\int_{t = 0}^{\,u} {\left( {u - t} \right)tdt} } = {1 \over {24}}\int_{u = 0}^{\,1} {u^{\,4} dx} = {1 \over {120}} \cr} $$ and surely you can continue for the others.