I'm not sure but something like the following can be thought of. ( I haven't gone through the Wikipedia code )
For every coordinate (x,y,z) find the sum of all elements from (0,0,0) to this element.
Call it S(0,0,0 to x,y,z) or S0(x,y,z).
This can be easily built by traversing the 3D matrix once.
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 element value + S0(x-1,y-1,z-1) + 3 faces (xy,yz,zx) )
Now for every query (x1,y1,z1) to (x2,y2,z2), first convert the coordinates so that x1,y1,z1 is the corner of the cuboid closest to origin and x2,y2,z2 is the corner that is farthest from origin.
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 )
(subject to corrections)