In this case you can multiple the rotation matrix of the region (region.rot) by the rotation matrices that perform the rotation around each axis.
| 1 0 0 |
X = | 0 A -B |
| 0 B A |
| C 0 -D |
Y = | 0 1 0 |
| D 0 C |
| E -F 0 |
Z = | F E 0 |
| 0 0 1 |
where
A, B - cos and sin of the rotation angle around X,
C, D - cos and sin of the rotation angle around Y,
E, F - cos and sin of the rotation angle around Z.
So for 90 degree rotation matrices will be the following:
xm = PhotoScan.Matrix( [[1,0,0],[0,0,-1],[0,1,0]] )
ym = PhotoScan.Matrix( [[0,0,-1],[0,1,0],[1,0,0]] )
zm = PhotoScan.Matrix( [[0,-1,0],[1,0,0],[0,0,1]] )
Hope it helps.