Python 3.8 (pre-release),  99   79   74   73   70   68  67 bytes

-20 bytes by removing the auxiliary t lambda, since it was only used once.
-5 bytes thanks to Lucenaposition.
-1 byte thanks to Arnauld.
-3 bytes thanks to l4m2.
-2 bytes by using the splat operator again.
-1 byte thanks to Lucenaposition.

lambda*a,c,d:sum(x-y*~y/2for x in range(*a)for y in range(c+x,d+x))

Takes input like f(inclusiveTopY, exclusiveBottomY, inclusiveLeftX=..., exclusiveRightX=...).
Try it online!

Pyth, 37 bytes

JEKE=GE=TEVrJKVr+GN+TN aY+/*H+1H2N;sY

Port of squareroot12621's answer. Had to work around the binary NOT.

Try it here

JEKE=GE=TE                  Assign variables J, K, G, T
VrJK                        For N in range(J, K)
    Vr+GN+TN                For H in range(G+N, T+N)
            aY+/*H+1H2N     Y.append(N + (H * (H + 1) / 2))
;sY                         Print sum(Y)

Jelly,  11  9 bytes

-2 thanks to Unrelated String (use of the table hyper, þ, to avoid separately taking the Cartesian product with Œp).

r+R;ɗþ/FS

A dyadic Link that accepts the top-left coordinate on the left ([row, column]) and the bottom-right coordinate on the right ([row, column]) and yields the sum.

Try it online! Or see the test-suite.

How?

The value at a given coordinate is the triangle number of the sum of the row and column (the minimum value on its running diagonal) plus its column (the number of steps taken along that diagonal).

r+R;ɗþ/FS - Link: [Top, Left], [Bottom, Right]
r         - {[Top, Left]} inclusive range {[Bottom, Right]} (vectorises)
              -> [[Top..Bottom], [Left..Right]]
      /   -   reduce by - f([Top..Bottom], [Left..Right]):
     þ    -     table of:
    ɗ     -       last three links as a dyad - f(R, C):
 +        -         {R} add {C} -> R+C
  R       -         range -> [1..R+C]
   ;      -         concatenate {C} -> [1,2,...,R+C,C]
       F  - flatten 
        S - sum

Haskell, 53 48 47 bytes

• -5 bytes thanks to xnor
• -1 byte thanks to Wheat Wizard ♦

f a b c d=sum$do i<-[a..c];j<-[b..d];j:[1..i+j]

Try it online!

The code above is equivalent to:

f a b c d = sum[h i j | i<-[a..c], j<-[b..d]]
h r c = sum[0..r+c] + c

• the helper function h computes the value for row r and column c
  • The value in the first column is equal to the number of elements in the upper triangle, i.e. the (anti-)diagonals above. Each diagonal has one more element than the preceding: $$h(r,0) =0+1+2+3+4+...+r=\sum_{i=0}^ri$$
  • To compute a general element we can compute the leftmost value of its (anti)diagonal (i.e. h(r+c,0)) and then go up by c steps (i.e. sum c): $$h(r,c)=h(r+c,0)+c$$
• [h i j | i<-[a..c], j<-[b..d]] applies h to all pairs (a,b), (a,b+1), ..., (a,d), (a+1,b), (a+1, b+1,), ..., (c,d)

Python, 61 59 bytes

f=lambda a,b,c,d:c<d and((c+b)**3-b)/6+b*c-f(b,a,c+(b<a),d)

Attempt This Online!

Takes top incl, bottom excl, left incl, right excl.

How?

Just seeing that @Jonathan Allan scooped me on this insight.

Based on the observation that the first column are just triangle numbers. Its partial sums are tetrahedral numbers and we can compute the sum of the first column as \$(\begin{smallmatrix}{\rm bottom}\\3\end{smallmatrix})-(\begin{smallmatrix}{\rm top}\\3\end{smallmatrix})\$. As for other columns, they are just shifted up and have a constant offset added, for example, the second column is the same as the first, except each value is shifted upwards one grid unit and incremented by one.

Desmos, 49 bytes

f(a,b,c,d)=∑_{x=a}^b∑_{y=c}^d((x+y)^2+3x+y)/2

Order: inclusiveLeftX, inclusiveRightX, inclusiveTopY, inclusiveBottomY

Try it online

Charcoal, 26 23 bytes

≔…ＮＮθＩΣＥ…ＮＮΣＥθ⁺λΣ…·⁰⁺ιλ

Try it online! Link is to verbose version of code. Takes inputs in the order exclusive column range, exclusive row range e.g. the test case 2, 0, 6, 1 becomes 0 2 2 7. Explanation:

≔…ＮＮθ

Read in the exclusive column range.

ＩΣＥ…ＮＮΣＥθ⁺λΣ…·⁰⁺ιλ

Map over the exclusive row range, calculating the value of each cell, then double sum over the results.

Python 3, 73 bytes

lambda a,b,c,d:(e:=c-a)*(f:=d-b)*(e*e+f*f+3*(~a-b-c-d)**2-12*(a+c)-17)/24

Try it online!

I tried golfing a closed formula (independently of ovs) that I first found using sympy, even though this code is longer than squareroot12621's solution which uses list comprehensions.

Inputs are inclusiveTopY, inclusiveLeftX, exclusiveBottomY, exclusiveRightX. This means (a,b,c,d) corresponds to the rectangle range(a,c) x range(b,d).

These half-open ranges have the advantage that c-a and d-b are factors, since the rectangles have dimensions of zero when a==c and b==d, without need for an off-by-1 correction. The differences c-a and d-b appear again later in the formula, so the code aliases them to e and f using walruses.

The formula is a bit nicer to see in a more symmetrical form:

$$ f(a,b,c,d) = \frac{1}{24}(c-a)(d-b) \cdot \\ \left(3(a+b+c+d)^2 + (c-a)^2 + (d-b)^2 + 6(d-c+b-a) - 14 \right)$$

79 bytes

lambda a,b,c,d:(c-a)*(d-b)*(3*(a+b+c+d)**2+(c-a)**2+(d-b)**2+6*(d-c+b-a)-14)/24

Try it online!

APL+WIN, 36 bytes

Prompts for the 3rd and 4th inter as specified in the examples followed by the 1st and 2nd. Index origin = 0

+/+/⎕↓n↑(⍳m)⊖⊃(+\⍳m)+⍳¨1+⍳m←2×↑n←1+⎕

Try it online! Thanks to Dyalog Classic

Python 3, 75 bytes

Just wanted to see how far a closed formula can go. Uses the same inclusive bounds as in the question, in slightly different order.

lambda a,b,c,d:(a+~b)*(c+~d)*((a+b)*(b+3*(c+d)/2+5)+d*(2+d+c)-c*~c+a*~-a)/6

Try it online!

Started with this formula generated by WolframAlpha:

$$ \sum_{x=a}^b \sum_{y=c}^d{x + \frac{\left(x + y\right) (x + y + 1)}2} = \\ \frac{1}{12} (a - b - 1) (c - d - 1) (2 a^2 + a (2 b + 3 c + 3 d + 8) + 2 b^2 + b (3 c + 3 d + 10) + 2 (c^2 + c d + c + d^2 + 2 d)) $$

C (gcc), 60 bytes

f(x,y,X,Y){x=X?x<X?f(y,x,Y,0)+f(x+1,y,X,Y):0:x-(y+=x)*~y/2;}

Try it online!

Simple port of l4m2's answer.

05AB1E, 15 bytes

ŸŠŸâεηO`D±;*-}O

Port of @l4m2's comment on @squareroot12621's Python answer, so make sure to upvote that answer, as well as @l4m2's Javascript answer!

Takes four inputs in the order \$b,d,a,c\$, all four inclusive.

Try it online or verify all test cases.

Explanation:

Ÿ             # Create a [b,d]-ranged list of the first two (implicit) inputs
 Š            # Triple-swap to push the next two (implicit) inputs to the stack
  Ÿ           # Create an [a,c]-ranged list of those two inputs as well
   â          # Cartesian product to create all pairs of these two lists
    ε         # Map over each pair [x,y]:
     ηO       #  Cumulative sum it:
     η        #   Get the prefixes of this pair: [[x],[x,y]]
      O       #   Sum each inner list: [x,x+y]
       `      #  Pop and push those two values to the stack
        D     #  Duplicate the top one
         ±    #  Bitwise-NOT: -(x+y)-1
          ;   #  Halve: (-(x+y)-1)/2
           *  #  Multiply the top two: (x+y)*((-(x+y)-1)/2)
            - #  Subtract the top two: x-(x+y)*((-(x+y)-1)/2)
    }O        # After the map: sum this list
              # (after which the result is output implicitly)

JavaScript (Node.js), 55 bytes

f=(x,y,X,Y)=>X?x<X&&f(y,x,Y)+f(x+1,y,X,Y):x-(y+=x)*~y/2

Try it online!