I'm not sure, but you can think of something like the following. (I did not go through the Wikipedia code)
For each coordinate (x, y, z), find the sum of all elements from (0,0,0) in this element. Call it S (0,0,0 to x, y, z) or S0 (x, y, z).
This can be easily built by going through a three-dimensional matrix.
S0( x,y,z ) = value[x,y,z] + S0(x-1,y-1,z-1) + S0( x,y,z-1 ) + S0( x, y-1, z ) + S0( x-1, y, z ) - S0( x-1, y-1, z ) - S0( x, y-1, z-1 ) - S0( x-1, y, z-1 )
(basically the value of the element + S0 (x-1, y-1, z-1) + 3 faces (xy, yz, zx))
Now for each query (x1, y1, z1) - (x2, y2, z2), first transform the coordinates so that x1, y1, z1 is the angle of the cuboid closest to the origin, and x2, y2, z2 is the angle that is further from the start.
S( (x1,y1,z1) to (x2,y2,z2) ) = S0( x2,y2,z2 ) - S0( x2, y2, z1 ) - S0( x2, y1, z2 ) - S0( x1, y2, z2 ) + S0( x1, y1, z2 ) + S0( x1, y2, z1 ) + S0( x2, y1, z1 ) - S0( x1, y1, z1 )
(as amended)