Inverse fourier transform by using a random phase and an original magnitude (matlab)
You image is not an $N \times M$ matrix. It's actually an $N \times M \times 3$ array of numbers. I assume that's why you're using fftn instead of fft2 to take its Fourier transform.

You have to either 


convert your image to grayscale; or
process each component independently.


To work get grayscale, we can naïvely average over the RGB colors like so:

img = mean(double(img),3)


To work with all three components independently, your phase array needs to have an extra dimension:

rand_phase= 4*pi*(rand(s(1),s(2), s(3))-0.5);