Domain distortion on a spherical texture (planetary surface)

Question

I am experimenting with procedural techniques to generate planetary terrain. I have close to zero knowledge so I started with something simple, using perlin noise and domain distortion. My goal is to wrap these textures on a sphere model (using openGL) so I tried to use spherical coordinates and it was partially successful (I learned from https://ronvalstar.nl/creating-tileable-noise-maps).

As long as I don't apply domain distortion I obtain the expected result with no discontinuity at the poles or anywhere:

But when I do apply the domain distortion, I get this:

You can see there are discontinuities at the poles, and also at a meridian.

I'm not sure to understand why the domain distortion is breaking my spherical mapping, I tried to apply a modified formula from here: https://iquilezles.org/www/articles/warp/warp.htm but it looks like the way I define the position vector is wrong in this case.

I do not know if spherical coordinates are adapted to this kind of problem. Can you please confirm this is the most efficient way to go?

If you have something else to suggest, please share. I have been reading a lot of things about cubemaps, (u, v) texture mapping, procedural texturing etc but did not really understand everything and how it might help me, also for the most part it looked quite sophisticated (intimidating for me).

Here is my code (the buggy part I guess is under the "# --- WARP (domain distortion) ---" line:

import numpy as np
from noise import pnoise3 as perlin3
from noise import pnoise2 as perlin2
import matplotlib.pyplot as plt
import noise
from PIL import Image
import time

# --- transformation X, Y (cartesian) -> a, b, c (spherical) ---
def spherical(x, y, naxis):

  fNx = (x + 0.5) / naxis
  fNy = (y + 0.5) / naxis

  fRdx= fNx * 2 * np.pi
  fRdy = fNy * np.pi
  fYsin = np.sin(fRdy + np.pi)

  a = naxis * np.sin(fRdx) * fYsin
  b = naxis * np.cos(fRdx) * fYsin
  c = naxis * np.cos(fRdy)

  return a, b, c

# --- MAIN ---
def procedural(naxis, WARP, octaves, persistence, lacunarity):

  t0 = time.time()
  M = []

  for i in range(len(WARP) + 1):
    M.append(np.zeros((naxis, naxis)))
  M = np.array(M)

  # --- perlin noise on spherical coordinates ---
  M[0] = np.array([[perlin3(float((spherical(i, j, naxis)[0]) / naxis),                          
                            float((spherical(i, j, naxis)[1]) / naxis),                          
                            float((spherical(i, j, naxis)[2]) / naxis),                           
                            persistence = persistence, lacunarity = lacunarity, octaves = octaves)    
                            for i in range(naxis)] for j in range(naxis)])

  # --- WARP (domain distortion) ---
  # f( p + f(p) )
  # p: spherical(i, j) -> a, b, c
  # components:  f( a + M[w-1] )
  #              f( b + M[w-1] )
  #              f( c + M[w-1] )

  for w in range(1, len(WARP) + 1):
    M[w] += np.array([[perlin3(float((M[w-1][i, j] * naxis * WARP[w-1] + spherical(i, j, naxis)[0]) / naxis), 
                               float((M[w-1][i, j] * naxis * WARP[w-1] + spherical(i, j, naxis)[1]) / naxis), 
                               float((M[w-1][i, j] * naxis * WARP[w-1] + spherical(i, j, naxis)[2]) / naxis),  
                               persistence = persistence, lacunarity = lacunarity, octaves = octaves)            
                               for i in range(naxis)] for j in range(naxis)])
 
  # --- normalization [0, 1] ---
  for n in range(len(WARP) + 1):
    M[n] = (M[n] - np.min(M[n])) / (np.max(M[n]) - np.min(M[n]))

  return M

cm = 'coolwarm'

naxis = 1024
WARP = [0.5, 1, 0.5]
#WARP = []
octaves = 24
persistence = 0.5
lacunarity = 2

R = procedural(naxis, WARP, octaves, persistence, lacunarity)

# --- save jpeg ---

im = Image.fromarray(R[-1] * 255)
im = im.convert('RGB')
im.save("your_file.jpeg")

fig = plt.figure()
ax = fig.add_subplot(111)
ax.imshow(R[-1], cmap = cm)
plt.show()

Answer 1

When applying the warp, do it on the 3D coordinates while you're performing the initial noise evaluation in the first place. I am not super familiar with the Numpy syntax here, but basically you would want to redefine your lambda function so that instead of generating noise(<x, y, z>) or noiseA(<x, y, z> + noiseB(<x, y, z>)) at a particular point, where <x, y, z> is converted from spherical coordinates, you want to go ahead and generate noise(warp(<x, y, z>)). Then perhaps warp(<x, y, z>) = warp(<x + noiseWarpX(<x, y, z>), y + noiseWarpY(<x, y, z>), z + noiseWarpZ(<x, y, z>)>). Adding noise to each individual coordinate is slightly biased, but it's a lot less biased than adding it to all three coordinates at once.

Speaking of bias, be wary of tutorials that recommend you to use "Perlin" noise. Actual Perlin noise, which you're using here, is an old function which produces a lot of 45 and 90 degree bias. Good simplex type noises are a newer approach, which typically produce less directionally aligned results. I would consider tutorials which discuss Perlin without at least any caveats or clarification on this matter, to be recirculating outdated information. IMO there's almost no reason to use actual Perlin noise nowadays unless you're doing something to address its grid bias.

I suggest using PyFastNoiseLite, with the OpenSimplex2 noise option and the FBm fractal option. For domain warping, FastNoiseLite has a dedicated, more efficient (and further bias reducing) domain warping function implemented. However it is not mirrored in PyFastNoiseLite, so you might either have to implement the wrapper yourself, or just use the three-noises technique.