IFFT 海洋呈圆形而不是尖状

Question

我使用 IFFT 生成的波相当上下，产生凹凸，而不是彼此分开然后回到一起，产生尖峰，如果这有意义的话。 我对 FFT/IFFT 没有太多经验，但据我所知，我的代码似乎正在做它应该做的事情。现在，由于我自己找不到任何圆形水的原因，这是我的问题：

有没有所有指南和论文都没有提到（或者我忽略了）的东西，使水变得尖刺而不是圆形？

我知道这是一个相当开放式的问题，但是由于我花了很多时间查看代码，没有发现任何缺陷，我只是希望有人可能有类似的经历或暗示可能会发生什么开。

我尝试实现 Jerry Tessendorf 介绍的 IFFT 算法来模拟海水。我基本上是从这篇paper中复制了代码。 我还将代码与本guide中的代码进行了比较。 经过几个小时的详细比较、一遍又一遍地阅读论文并调整参数后，我的代码生成的水仍然看起来是圆形/凹凸不平的，而不是尖刺的。 我在 python 中依次实现了它，并在带有计算着色器的虚幻引擎中实现了它。在这两种情况下，我的水看起来都不像我预期的那样。

论文提到 x 和 y 偏移（假设 z 是高度轴）很重要，但我已经实现了这一点。 我发现当我翻转 z 轴（将其乘以 -1）时，水看起来更符合我的预期，但是计算出的法线贴图将不再适合水，因为它仍然适合圆形水。

以下是一些参考图片：

圆形/凹凸不平的水：

尖水（高度乘以-1）：

编辑： 这是我用 python 编写的代码：

幅度部分：

def ph(k):
    L_ = (windspeed * windspeed) / g
    mag = np.linalg.norm(k)
    if(mag < 0.0001):
        mag = 0.0001
    magSq = mag * mag

    wave_exp = 6.0

    Ph = (A/(magSq*magSq) * pow(np.dot(normalize(k),normalize(wind_direction)), wave_exp) * math.exp(-(1.0/(magSq * L_ * L_))))

    return Ph * math.exp(-magSq * pow(l, 2.0))



def generate_h0k():
        for _x in range(N):
        for _y in range(N):
            x = np.array([_x, _y])  - N / 2.0
            
            k = x * 2 * math.pi / L
            
            L_ = windspeed * windspeed / g
            mag = np.linalg.norm(k)
            if mag < 0.00001:
                mag = 0.00001

            magsq = mag * mag
            
            h0k = math.sqrt((A / (magsq * magsq)) * pow(np.dot(normalize(k),normalize(wind_direction)), 6.0) * math.exp(-(1.0 / (magsq * L_ * L_))) * math.exp(-magsq * pow(L/2000.0, 2.0))) / math.sqrt(2.0)
            h0minusk = math.sqrt((A / (magsq * magsq)) * pow(np.dot(normalize(-k),normalize(wind_direction)), 6.0) * math.exp(-(1.0 / (magsq * L_ * L_))) * math.exp(-magsq * pow(L/2000.0, 2.0))) / math.sqrt(2.0)

            h0k = math.sqrt(ph(k)) / math.sqrt(2)
            h0minusk = math.sqrt(ph(-k)) / math.sqrt(2)

            h0k = clamp(h0k, -4000, 4000)
            h0minusk = clamp(h0minusk, -4000, 4000)


            gauss_random = np.array([random.gauss(0, 1), random.gauss(0, 1), random.gauss(0, 1), random.gauss(0, 1)])

            
            tilde_h0_k[_x, _y] = gauss_random[0:2] * h0k
            tilde_h0_minusk[_x, _y] = gauss_random[2:4] * h0minusk


def dispersion(k):
    w_0 = 2 * math.pi / 200
    return math.floor(math.sqrt(g * np.linalg.norm(k) ) / w_0) * w_0

def generate_hkt(t):
    for _x in range(N):
        for _y in range(N):
            x = np.array([_x, _y])  - N / 2

            k = np.array([2.0 * math.pi * x[0]/L, 2.0 * math.pi * x[1]/L])
            
            mag = np.linalg.norm(k)
            if mag < 0.00001:
                mag = 0.00001
            
            w = dispersion(k)

            h0k = complex(tilde_h0_k[_x, _y][0], tilde_h0_k[_x, _y][1])
            h0_minusk_conj = comp_conj(complex(tilde_h0_minusk[_x, _y][0], tilde_h0_minusk[_x, _y][1]))

            cos_wt = math.cos(w * t)
            sin_wt = math.sin(w * t)

            exp_iwt = complex(cos_wt, sin_wt)
            exp_iwt_inv = complex(cos_wt, -sin_wt)

            h_k_t_dy = comp_add(comp_mul(h0k, exp_iwt), comp_mul(h0_minusk_conj, exp_iwt_inv))

            dx = complex(0, -k[0]/ mag)
            h_k_t_dx = comp_mul(dx, h_k_t_dy)

            dy = complex(0, -k[1]/mag)
            h_k_t_dz = comp_mul(dy, h_k_t_dy)
            
            
            tilde_hkt_dy[_x, _y] = np.array([h_k_t_dy.real, h_k_t_dy.im])
            tilde_hkt_dx[_x, _y]  = np.array([h_k_t_dx.real, h_k_t_dx.im])
            tilde_hkt_dz[_x, _y]  = np.array([h_k_t_dz.real, h_k_t_dz.im])

FFT部分：

 def generate_init_bf_first_stage_indices():
    for i in range(N):
        i = np.uint8(i)
        b = int('{:08b}'.format(i)[::-1], 2) # bit reversing
        bf_first_stage_indices.append(b)

    def generate_init_bf_texture():
        generate_init_bf_first_stage_indices()
    
        for _x in range(bf_height):
            for _y in range(N):
                
                k = (_y * (float(N) / pow(2, _x+1))) % N
    
                twiddle = complex(math.cos(2.0 * math.pi * k / float(N)), math.sin(2.0 * math.pi * k / float(N)))
                butterflyspan = int(pow(2, _x))
    
                butterflywing = 0
                if _y % pow(2, _x+1) < pow(2, _x):
                    butterflywing = 1
    
                if _x == 0:
                    if(butterflywing == 1):
                        bf_twiddle_indices[_x, _y] = np.array([twiddle.real, twiddle.im, int(bf_first_stage_indices[int(_y)]), int(bf_first_stage_indices[int(_y+1)])])
                    else:
                        bf_twiddle_indices[_x, _y] = np.array([twiddle.real, twiddle.im, int(bf_first_stage_indices[int(_y-1)]), int(bf_first_stage_indices[int(_y)])])
                else:
                    if(butterflywing == 1):
                        bf_twiddle_indices[_x, _y] = np.array([twiddle.real, twiddle.im, int(_y), int(_y + butterflyspan)])
                    else:
                        bf_twiddle_indices[_x, _y] = np.array([twiddle.real, twiddle.im, int(_y - butterflyspan), int(_y)])

def switch_pingpong(_pingpong):
    if _pingpong == 1:
        return 0
    else:
        return 1

def compute_horizontal_butterflies():
    global pingpong_0
    global pingpong_1
    global current_pingpong
    for current_stage in range(bf_height):
        for _y in range(N):
            for _x in range(N):
                
                if current_pingpong == 0:
                    data = bf_twiddle_indices[current_stage, _x]
                    p_ = np.array([pingpong_0[int(data[2]), _y][0], pingpong_0[int(data[2]), _y][1]]) # red and green part seperately
                    q_ = np.array([pingpong_0[int(data[3]), _y][0], pingpong_0[int(data[3]), _y][1]])
                    w_ = np.array([data[0], data[1]])

                    p = complex(p_[0], p_[1])
                    q = complex(q_[0], q_[1])
                    w = complex(w_[0], w_[1])

                    h = comp_add(p, comp_mul(w, q))
                   
                    pingpong_1[_x, _y][0] = h.real
                    pingpong_1[_x, _y][1] = h.im

                else:
                    data = bf_twiddle_indices[current_stage, _x]
                    p_ = np.array([pingpong_1[int(data[2]), _y][0], pingpong_1[int(data[2]), _y][1]]) 
                    q_ = np.array([pingpong_1[int(data[3]), _y][0], pingpong_1[int(data[3]), _y][1]]) 
                    w_ = np.array([data[0], data[1]])

                    p = complex(p_[0], p_[1])
                    q = complex(q_[0], q_[1])
                    w = complex(w_[0], w_[1])

                    h = comp_add(p, comp_mul(w, q))

                    pingpong_0[_x, _y][0] = h.real
                    pingpong_0[_x, _y][1] = h.im
                

        current_pingpong = switch_pingpong(current_pingpong)
        


def compute_vertical_butterflies():
    global pingpong_0
    global pingpong_1
    global current_pingpong

    for current_stage in range(bf_height):
        for _x in range(N):
            for _y in range(N):
                if current_pingpong == 0:
                    data = bf_twiddle_indices[current_stage, _y]
                    p_ = np.array([pingpong_0[_x, int(data[2])][0], pingpong_0[_x, int(data[2])][1]]) 
                    q_ = np.array([pingpong_0[_x, int(data[3])][0], pingpong_0[_x, int(data[3])][1]])
                    w_ = np.array([data[0], data[1]])

                    p = complex(p_[0], p_[1])
                    q = complex(q_[0], q_[1])
                    w = complex(w_[0], w_[1])

                    h = comp_add(p, comp_mul(w, q))

                    pingpong_1[_x, _y][0] = h.real
                    pingpong_1[_x, _y][1] = h.im
                else:
                    data = bf_twiddle_indices[current_stage, _y]
                    p_ = np.array([pingpong_1[_x, int(data[2])][0], pingpong_1[_x, int(data[2])][1]])
                    q_ = np.array([pingpong_1[_x, int(data[3])][0], pingpong_1[_x, int(data[3])][1]])
                    w_ = np.array([data[0], data[1]])

                    p = complex(p_[0], p_[1])
                    q = complex(q_[0], q_[1])
                    w = complex(w_[0], w_[1])

                    h = comp_add(p, comp_mul(w, q))

                    pingpong_0[_x, _y][0] = h.real
                    pingpong_0[_x, _y][1] = h.im
                    
        current_pingpong = switch_pingpong(current_pingpong)

def compute_fourier_transforms(amplitudes):
    global current_pingpong
    pingpong_0[:,:, 0] = amplitudes[:, :, 0]
    pingpong_0[:, :, 1] = amplitudes[:, :, 1]


    current_pingpong = 0
    compute_horizontal_butterflies()
    compute_vertical_butterflies()


    return copy.deepcopy(pingpong_0)


def generate_displacement(amplitudes, multiplier=1):
    new_displacement = np.zeros(shape=(N, N))
    for _x in range(N):
        for _y in range(N):
            perms = [1.0, -1.0]
            index = int(int(_x + _y)) % 2
            perm = perms[index]

            h = amplitudes[_x, _y][0]
            
            new_displacement[_x, _y] = multiplier * perm * float(h)/float(N*N)
            
    return new_displacement

一切如何执行：

def init():
    random.seed(1234)
    generate_h0k()

    generate_init_bf_texture()


if __name__ == "__main__":
    init()

    duration_seconds = 5
    deltatime = 1/30
    frame_amount = int(duration_seconds // deltatime)
    
    for i in range(frame_amount):
        time = i * deltatime

        generate_hkt(time)        
        
        height_amplitude = compute_fourier_transforms(tilde_hkt_dy)
        x_amplitude = compute_fourier_transforms(tilde_hkt_dx)
        z_amplitude = compute_fourier_transforms(tilde_hkt_dz)

        height_disp = generate_displacement(height_amplitude)
        x_disp = generate_displacement(x_amplitude)
        z_disp = generate_displacement(z_amplitude)

参数：

N = 256
L = 1000
l = 0.5
windspeed = 26
wind_direction = [1, 1]
g = 9.81
A = 4

bf_height = int(math.log2(N))

数据结构的初始化：

bf_height = int(math.log2(N))
bf_first_stage_indices = []
bf_twiddle_indices = np.zeros(shape=(bf_height, N, 4))#, dtype=ArrayType)

tilde_h0_k = np.zeros(shape=(N,N, 2))
tilde_h0_minusk = np.zeros(shape=(N, N, 2))

tilde_hkt_dx = np.zeros(shape=(N, N, 2))
tilde_hkt_dy = np.zeros(shape=(N, N, 2))
tilde_hkt_dz = np.zeros(shape=(N, N, 2))

pingpong_0 = np.zeros(shape=(N, N, 2))
pingpong_1 = np.zeros(shape=(N, N, 2))
current_pingpong = 0

实用工具：

class complex():
    def __init__(self, _real, _im):
        self.real = _real
        self.im = _im

def comp_conj(_complex):
    return complex(_complex.real, -_complex.im)

def comp_mul(c0, c1):
    new_real = c0.real * c1.real - c0.im * c1.im
    new_im = c0.real * c1.im + c0.im *c1.real

    return complex(new_real, new_im)

def comp_add(c0, c1):
    new_real = c0.real + c1.real
    new_im = c0.im + c1.im
    return complex(new_real, new_im)

def clamp(num, min_value, max_value):
    num = max(min(num, max_value), min_value)
    return num


def normalize(in_vec):
    if(np.linalg.norm(in_vec) == 0):
        return in_vec
    return in_vec / np.linalg.norm(in_vec)