Domain distortion on a spherical texture (planetary surface)

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.imshow(R[-1], cmap = cm)
plt.show()

• I don't think wrapping an extra fbm deep will remove a discontinuity that's already present in the input. It looks like somewhere here, there's a 2D/non-wrapping input being used instead of a 3D wrap-around input, and seams in the input (where i wraps around from n to 0) lead to seams in the output. But I'm not familiar enough with this syntax to identify where that's happening. Jan 12, 2021 at 18:49
• @DMGregory Looking at it again, I agree, my original assessment was incorrect & I have removed it. Since the original texture wraps, could OP produce additional randomized instances of those & use them as the basis for a spherically wrapping displacement? Jan 13, 2021 at 2:22

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.
• that would be sphereWarp(<x, y, z>) = radius * normalize(warp(<x, y, z>)) Jan 13, 2021 at 16:43