05AB1E, 75 bytes

-.øε0δ.ø©}ε®Ù.ø}D€øDø€ø««εÁÁεN3@iD0ÚPi1V}YiγÁεN2@id}}À˜}}}ÀÀ2V}€`€ü2€`ʒË≠}g

Dang, this was hard in 05AB1E.. But it works now. 05AB1E and matrices are already a bad combination, so add an additional dimension and it's a complete disaster, haha..

Try it online or verify all test cases.

Explanation:

Step 1: Surround the entire input 3D matrix with a layers of empty cells (0s) in each dimension:

-                    # Transform all values into 0s by subtracting the values in the
                     # (implicit) input 3D-matrix by the values in the (implicit) input
 .ø                  # Surround the (implicit) input-matrix with this 2D-matrix of 0s as
                     # both leading and trailing item
   ε                 # Map each 2D matrix of the 3D matrix to:
     δ               #  For each row of the 2D matrix:
    0 .ø             #   Surround it with a leading and trailing 0
        ©            #  Store the modified 2D matrix in variable `®` (without popping)
   }ε                # After the map: map over each 2D matrix in the 3D matrix again:
     ®Ù              #  Uniquify the last 2D matrix that was stored in `®`,
                     #  so we'll have a row of 0s wrapped inside a list
       .ø            #  Surround each 2D matrix with this row of 0s
    }                # And close this map as well

(Note: the z-axis actually contains two surrounding empty cells instead of one, but this doesn't really matter for the rest of the program.)

Step 2: Get a list of all strings of cells along the x, y, and z axes respectively:

D                    # Duplicate the current 3D-matrix, which of itself already contains
                     # all strings of cells along the x-axis
 €                   # Map each 2D matrix of the 3D matrix to:
  ø                  #  Zip/transpose; swapping rows/columns
D                    # Duplicate as well, which are the strings of cells along the y-axis
 ø                   # Zip/transpose; swapping rows/columns of this 3D matrix
  €                  # Map each 2D matrix of the 3D matrix to:
   ø                 #  Zip/transpose; swapping rows/columns
                     # And we now also have the strings of cells along the z-axis
««                   # Merge all three lists together

This will result in a 3D matrix with three inner 2D matrices (one for each dimension), which are each lists of strings of cells.

Step 3: Fill all inner bubbles with 1s:

ε                    # Map each 2D matrix of the 3D matrix to:
 ÁÁ                  #  Rotate the rows of the matrix twice towards the left
   ε                 #  Map each string of cells in the current 2D matrix to:
    N3@i             #   If the 0-based index is >= 3:
        D            #    Create a copy of the string of cells
         0Ú          #    Remove all leading and trailing empty cells
           Pi  }     #    If there are now only filled cells left:
             1V      #     Set flag `Y` to 1
         Yi          #    If flag `Y` is 1:
           γ         #     Split the string of cells into groups of equal adjacent values
            Á        #     Rotate these groups once towards the left
             ε       #     Map each group to:
              N2@i } #      If the 0-based index is >= 2:
                  d  #       Fill all empty cells (1 remains 1, 0 becomes 1)
             }À      #     After the map: rotate the groups back to the right
               ˜     #     And flatten it to a single string of cells again
   }}}ÀÀ             #  After the map: rotate the rows twice back towards the right
        2V           #  Reset flag `Y` back to 2 for the next iteration
}                    # Close the map

We basically skip the first and last strings of cells, since we know those are surrounding layers of empty cells we added in step 1. In addition, we also don't want to modify the second and second to last strings of cells, since those are the outer layers of the initial input 3D matrix. We do however want to start checking from the second string of cells onward until we find a solid string of filled cells (minus the surrounding empty cells). For all strings of cells after we've encountered such a solid string of filled cells, we want to transform them into solid strings of filled cells as well (minus the surrounding empty cells) to fill the bubble.

Step 4: Now that we've filled the bubbles, we want to get a list of all pairs of cells:

€`                   # Flatten the 3D matrix containing the three list of strings of
                     # cells one level down to a single list of strings of cells
  €                  # Map each string of cells to:
   ü2                #  Create overlapping pairs of cells
     €`              # And flatten this list of list of pairs one level down as well to a
                     # list of pairs

Step 5: Filter out any pairs of two empty or two filled cells, so we only have pairs containing one of each:

ʒ                    # Filter this list of paired cells by:
 Ë≠                  #  Check that both values in the pair are NOT the same
}                    # Close the filter

Step 6: Get the amount of pairs left containing both a filled and empty cell, and output it as result:

g                    # Pop and push the length of the filtered list
                     # (after which it is output implicitly as result)

Try it online with each of these steps output separately.

Python 3 with NumPy, 177 161 bytes

-16 bytes thanks to @fireflame241 !

f=lambda l:g(pad(pad(l,1)-2,1)+2,1,1,1)
def g(l,*i):l[i]+=2;return l[i]%2if l[i]-2else sum(g(l,*(t*d+i))for d in eye(3,3,0,int)for t in[1,-1])
from numpy import*

Try it online!

Big idea:

DFS over all outer empty cells. Every time an outer empty cell touches a cube, adds 1 to the counter.

Explanation:

• 0 denotes air (empty cell), odd positive numbers denote walls, and even positive numbers denotes paint.
• First, pad everything with a layer of 0 (air), so all outer air cells are connected: pad(l,1)
• Then, pad everything with a layer of 2 (paint), which acts as a barrier to prevent out of bound search later. To do this, subtracts 2 from all cells, pad everything with 0, then adds 2 back: pad(arr - 2, 1) + 2
• Start DFS at l[1,1,1], which is guaranteed to be an outer air cell.
• At each DFS step (function g):
  • If current cell is paint, stops recursion.
  • If current cell is wall, adds 1 to counter, and stops recursions.
  • If current cell is air, changes it to paint, and recurs on the 6 neighbors.

APL (Dyalog Unicode), ⁵⁵ 48 bytes

≢⍸↑2≠/¨⊢∘g3⍴⊂2=2(g⊣(⌈∧⊢)/,)⍣6⍣≡(1,g←⍉⍤2⍉∘⌽)⍣6~⎕

Try it online!

-7 bytes thanks to @ngn.

Improvements:

• 2 3 1⍉ → ⍉⍤2⍉: Replace "cycle the axes once" with "swap 1st and 3rd axes, then 2nd and 3rd"
• {⍵(g⍵)(g g⍵)} → ⊢∘g3⍴⊂: A scan that ignores the left arg and applies g to the right arg, so it works like this:

3⍴⊂x gives (x x x)
⊢∘g3⍴⊂x gives (x)(x ⊢∘g x)(x ⊢∘g x ⊢∘g x)
which is the same as (x)(g x)(g g x) because:
  x ⊢∘g x
→ x ⊢ g x
→ x ⊢ (g x)
→ g x

---

APL (Dyalog Unicode), 55 bytes

{≢⍸↑2≠/¨⍵(g⍵)(g g⍵)}2=2(g⊣(⌈∧⊢)/,)⍣6⍣≡(1,g←2 3 1⍉⌽)⍣6~⎕

Try it online!

A full program that takes a 3D array. Uses the flood fill already used here. Another key idea is g←2 3 1⍉⌽, which effectively cycles through all six sides when applied with ⍣6 (repeat six times).

How it works

{≢⍸↑2≠/¨⍵(g⍵)(g g⍵)}2=2(g⊣(⌈∧⊢)/,)⍣6⍣≡(1,g←2 3 1⍉⌽)⍣6~⎕

~⎕                ⍝ Logical negation of the input
(1,g←2 3 1⍉⌽)⍣6   ⍝ Pad with a layer of ones on all six sides
2(g⊣(⌈∧⊢)/,)⍣6⍣≡  ⍝ Flood fill from the outside, changing 1s to 2s:
2(        ,)      ⍝   Prepend 2 on the last axis
   ⊣(⌈∧⊢)/        ⍝   Pairwise lcm(max(x,y),y) over the last axis
                  ⍝   Effectively, propagate 2 to an adjacent 1 on the right
  g               ⍝   Cycle the orientation once
            ⍣6⍣≡  ⍝   Repeat 6 times until the flood fill is complete
2=                ⍝ Map 2s to 1s, and anything else to 0s
{⍵(g⍵)(g g⍵)}     ⍝ Construct 3 arrays so that each axis becomes the last axis
2≠/¨              ⍝ Extract faces (where 0 and 1 are adjacent) for each array
≢⍸↑               ⍝ Count ones in all arrays

JavaScript (ES6),  295  291 bytes

a=>a.map((s,z)=>s.map((r,y)=>r.map((v,x)=>v|!(g=(x,y,z,R=a[z]&&a[z][y])=>R&&1/R[x]?R[x]?0:R[x]++|[0,1,2,3,4,5].some(n=>(i=n&1||-1,g(n&6?x:x+i,n&2?y+i:y,n&4?z+i:z)))|--R[x]:1)(x,y,z)))).map((s,z,a)=>s.map((r,y)=>r.map((v,x)=>n+=v&&!r[x+1]+!((q=s[y+1])&&q[x])+!((q=a[z+1])&&q[y][x]))),n=0)|n*2

Try it online!

NB: This is a bit too slow to reliably complete the 6^th test case on TIO.

Python 2, 493 bytes

A=lambda*z:0<sum(abs(a-b)for a,b in zip(*z))<2
R=lambda c:reduce(lambda a,b:a|b,c)
def C(c,d,h,w):
 a=[[{(i/w/h,i/w%h,i%w)}for i in range(d*h*w)if c[i]-v]for v in[1,0]]
 for r in a:
	i=0
	for j in range(len(r)**2):i=j/len(r);c=[[f for f in r[i:]if any(A(j,k)for k in f for j in r[i])^j]for j in[0,1]];r[i:]=(c[0]and[R(c[0])])+c[1]
 a[0]=[s for s in a[0]if all(0<e[i]<[d,h,w][i]-1for i in[0,1,2]for e in s)]
 p,q=[sum(6-sum(A(x,y)for x in r)for y in r)for r in[k and R(k)for k in a]]
 print q-p

Try it online!

Takes input as a flattened array along with depth, height, and width.

How?

1. Find 6-connected components of 0s and of 1s
2. Remove the components of 0s which contain a 0 on the outer edge
3. Take 6 times the number of 1s minus the number of 1s that border each other to get the number of 1s that are exposed to any 0. This includes 0s on the inside (internal 0s / air pockets), so:
4. Subtract (6 times the number of internal 0s minus the number of internal 0s that border each other) to get the number of internal 0s that are exposed to any 1. This subtracts all of the faces on the inside.
5. Done!

# Are the arguments adjacent via 6-connectivity?
A=lambda *z:0<sum(abs(a-b)for a,b in zip(*z))<2

R=lambda c:reduce(lambda a,b:a|b,c)

def C(c,d,h,w):
    a=[
        [
            {(i/w/h,i/w%h,i%w)}
            for i in range(d*h*w)
            if c[i]-v
        ]
        for v in[1,0]

    ]
    # a[0]: set of coordinates of all 0s
    # a[1]: set of coordinates of all 1s
    # Find connected components:
    for r in a:
        i=0
        for j in range(len(r)**2):
            # for each index i
            i=j/len(r);
            # do len(r) passes:
            # c[0]: all components with index > i+1 that are adjacent to component i
            # c[1]: all components with index > i+1 that are not adjacent to component i
            c=[
                [f for f in r[i:]if any(A(j,k)for k in f for j in r[i])^j]
                for j in[0,1]
            ];
            # Replace components i and higher with:
            r[i:]=(
                # If c[0] is nonempty, then the union of c[0]
                c[0]and[R(c[0])]
            )+c[1] # append c[1]
    # a[0]: set of connected components of 0s
    # a[1]: set of connected components of 1s
    # Remove all of a[0] that border the outside:
    a[0]=[
        # Filter for:
        s for s in a[0]if
        all(
            # The coordinates along each axis are between 1 and that axis's length minus 2, inclusive
            0<e[i]<[d,h,w][i]-1
            for i in[0,1,2]
            # For all points
            for e in s
        )
    ]
    # a[0] now: set of connected components of 0s that do not border the outside
    p,q=[
        sum(
            6- # cube contributes 6 sides
            sum(A(x,y)for x in r) # minus the number of adjacent cells
            for y in r # for each cube
        )
        for r in # for each connected component
        [k and R(k)for k in a]
    ]
    print q-p

MATL, 37 36 35 bytes

e7BYa~6&1ZIt1)-tz6*yZybfX[hhtZPq~z-

Input is a row vector of zeros and ones, and a row vector of three integers with dimensions from inner nesting level to outer.

Try it online! Or verify all test cases.

Explanation

The code initially adds a frame of empty space around the 3D array. Any cell that is not space connected to that frame is filled. This has the effect of filling any holes in the original 3D shape.

The number of painted faces is the number of cubes in that filled shape times 6, minus the number of cubes that touch some other cube (two cubes touching means a face is not accessible to the paint; pairs are counted twice).

To detect which cubes touch, all pairwise distances between the cubes are computed, and two cubes touch if their distance is 1.

e      % Implicit inputs: vector or zeros and ones, and 3-element vector specifying
       % size along each dimension. Reshape the first according to the second. This
       % produces the 3D array
7BYa   % Pad the 3D array with a frame of zeros along the three dimensions
~      % Negate. This changes 0 to 1 and vice versa (*)
6&1ZI  % Label connected components using 6-connectivity in 3D (so no diagonals)
       % This considers the zeros in (*) as background, and nonzeros as foreground.
       % The foreground cells are labelled with a different integer according to
       % indicate the component. There will be an outer component that will include
       % the added frame and any space surrounding the shape, and potentially more
       % components if the shape has inner holes
t1)-   % Duplicate Subtract the upper-right-front entry from each entry. This
       % makes the outer connected component (originally the space surrounding the
       %  shape) equal to 0, and other components or brackground become nonzero. 
       % So now the shape plus any inner holes in it are nonzero (**)
tz6*   % Duplicate. Number of nonzeros times 6. This is the maximum number of faces
       % of unit cubes that could be painted (some won't actually get pointed,
       % namely any face that touches any other face) (***)
yZy    % Duplicate from below: pushes a copy of (**). Get its size as a length-3
       %  vector
bf     % Bubble up: moves the original copy of (**) to the top. Push linear indices
       % of its nonzero entries. Linear indices run down, then accros (left to
       % right), then front to bottom
X[     % Convert linear indices to a set of three indices assuming an array of the
       % specified size. Gives three column vectors
hh     % Concatenate veftically twice. This gives a 3-column matrix where each row
       % contains the coordinates of a cube in (**)
tZP    % Duplicate. Pairwise distances between rows of the 3-column matrix and
       % those of its copy
q~     % Subtract 1, negate. This gives 1 for distances equal to 1, and 0 otherwise
z      % Number of nonzeros
-      % Subtract from (***). Implicit display

Wolfram Language (Mathematica), 103 102 bytes

Count[x=ImageData@FillingTransform@Image3D@#~ArrayPad~1;x~Differences~#&/@{1,{0,1},{0,0,1}},1.|-1.,4]&

Try it online!

To use FillingTransform (replace all the inner 0 with 1), I have to convert the data to Image3D and convert it back. The rest is just to count the number of nonzero consecutive differences etc.