2-dimensional Fourier transform interpretation I have this image of a duck:

Taking the FFT of this image, I can generate two new images. A magnitude spectrum (right) and phase spectrum (bottom left).
I am trying to re-construct the image of the duck MANUALLY by adding together the 65,535 sine waves encoded in this information. I am running into a problem where my duck is reconstructed but looks a little funny.
If the 2D Fourier transform is as follows:
$f(x,y) = \sum\sum F[h,k]e^{-2\pi i(hx+ky+\phi)}$
I was hoping to be schooled in my interpretation of this. My sine wave amplitude is $F[h,k]$. My phase is $\phi$. My indices for the wave (the number of time it 'cuts' each respective axis are $h, k$.
So, what the magnitude spectrum is showing is:


*

*The center is $h,k = 0,0$. The zero frequency component intensity represents the average intensity of the entire image.

*Each coordinate, for example $(2, 2)$ from the center represents a wave with $(h, k) = (2, 2)$. The intensity of this pixel is an indication of the amplitude of the wave. 

*Similarly, in the phase spectrum, the values range between $(0, \pi)$ and represent the offset of this sine wave from the origin.


The code I am using to doing the following is below. I SHOULD be able to perfectly reconstruct the image, but I keep getting some weird looking duck. The following is the duck I get after adding together 65,535 waves:

I was hoping somebody could point out my error. Thank you!
EDIT
If I only add the first 128 rows of my 2D FFT arrays, I get a duck that looks closer to the original. I feel like it has something to do with the fact that you only need 1/2 the information from these spectra, but I am adding together some of the wrong info. Like I need to cut each spectrum diagonally (about the symmetrical axis) and add only those waves...
EDIT
Final product...

import numpy as np
import matplotlib.pyplot as plt
import cv2

#Wave generator
def waveGenerator(h, k, a, phi, res):
    mesh = np.fromfunction(lambda x, y: a*np.sin(2*np.pi*h*x/res + 2*np.pi*k*y/res + phi), shape = (res,res))
    return mesh

# Import some image...
duck = cv2.imread("data/fourier_duck.png", 0)

#Take the fft of the duck
duck_fft = np.fft.fft2(duck)
duck_fshift = np.fft.fftshift(duck_fft)
magnitude_duck = np.log(np.abs(duck_fshift) + 1)
unshifted = np.abs(duck_fft)
phases = np.angle(duck_fft)

#THE RECONSTRUCTED IMAGE (MANUALLY)
recon_img = np.full((256, 256), unshifted[0][0])

for h in range(len(unshifted)):
    for k in range(len(unshifted[h])):
        recon_img = np.add(recon_img, waveGenerator(h, k, unshifted[h][k], phases[h][k], 256))
        print(str(k) + ',' + str(h))

plt.imshow(recon_img)
plt.show()

 A: Try changing:
a*np.sin
to 
a*np.cos  
in your waveGenerator function.  When the phase is 0, the real part of the wavefunction should be maximum which corresponds to a cosine: i.e. cos(0)=1.  The rest seems to be okay, other than scaling issues I believe.  I hope this helps.
Simplified test program below adapted from the question:
import numpy as np
import matplotlib.pyplot as plt
#import cv2

#Wave generator
def waveGenerator(h, k, a, phi, res):
    mesh = np.fromfunction(lambda x, y: a*np.cos(2*np.pi*h*x/res + 2*np.pi*k*y/res + phi), shape = (res,res))
    return mesh

# Import some image...
#duck = cv2.imread("data/fourier_duck.png", 0)
n_size = 64
duck = np.zeros((n_size,n_size))
duck[n_size//2][n_size//2]=1.0
duck[(n_size+10)//2][(n_size+5)//4]=2.0
plt.imshow(duck)
plt.show()

#Take the fft of the duck
duck_fft = np.fft.fft2(duck)
duck_fshift = np.fft.fftshift(duck_fft)
magnitude_duck = np.log(np.abs(duck_fshift) + 1)
unshifted = np.abs(duck_fft)
phases = np.angle(duck_fft)

#THE RECONSTRUCTED IMAGE (MANUALLY)
recon_img = np.full((n_size, n_size), unshifted[0][0])

for h in range(len(unshifted)):
    for k in range(len(unshifted[h])):
        recon_img = np.add(recon_img, waveGenerator(h, k, unshifted[h][k], phases[h][k], n_size))

plt.imshow(recon_img)
plt.show()

