---

---

Voici une méthode versatile, qui se révèle également efficace contre les moirés, pour un coût plus faible. Je ne l'ai pas inventée.

Au fur et à mesure qu'on calcule les points z_n d'une orbite, on peut également calculer la dérivée d_n de la fonction qui au premier point z₀ associe le n+1-ième z_n. On a les relations de récurrence $$\begin{align} z_{n+1}&=f(z_n)\\ d_0&=1\\ d_{n+1}&=f'(z_n)d_n \end{align}$$

Le piège est illustré ci-dessous.

(boucle d'itération)
  ...
  z = f(z)
  d = f_prime(z)*d # oups...
  ...  

Si z contenait z_n et d contenait d_n alors maintenant z contient f(z_n) et d contient f'(z_n+1)d_n au lieu de f'(z_n)d_n.

Solution 1 : changer l'ordre des lignes (ça marche ici)

(boucle d'itération)
  ...
  d = f_prime(z)*d
  z = f(z)
  ...  

Solution 2 : utiliser des variables temporaires (c'est plus général)

(boucle d'itération)
  ...
  new_z = f(z)
  new_d = f_prime(z)*d
  z=new_z
  d=new_d
  ...  

Note : on peut de même calculer les dérivées d'ordre supérieur au fur et à mesure des itérations mais on ne les utilisera pas ici.

L'image d'un disque de rayon r par fⁿ est approximativement un disque de rayon d_n×r. Cette approximation linéaire n'est en théorie bonne que pour d_n≠0 et r suffisamment petit. L'idée est ici de faire comme si c'était le cas tout le temps, quand r est de l'ordre de la taille d'un pixel.

Remarque : En fait en dynamique holomorphe en dimension 1, on est fréquemment en situation de distorsion bornée, donc cette procédure donne souvent le bon ordre de grandeur pour la taille de l'image d'un tel disque.

Dans l'algo ci-dessous je choisis r = la longueur du côté d'un pixel, ce qui fait que le disque déborde bien du pixel (son aire est ≈ 3× plus grande). J'itère et considère le disque de centre z_n et de rayon d_n×r. S'il contient un point de l'ensemble de Julia alors je considère le pixel initial est proche du Julia et je le colorie en conséquence. En effet, sauf cas exceptionnels, il y a une préimage de ce disque dont la taille est comparable au disque initial B(z₀,r). Celle-ci contient un point de l'ensemble de Julia car, comme souvent avec les ensembles définis dynamiquement :

Mais on a mieux :

Ce qui est encourageant : peut-être suffit-il de tester un seul point de J ? On procède donc comme suit:

• On trouve/choisit un point x ∈ J.
• Pour chaque pixel, correspondant à un certain z₀, on calcule z_n et d_n pour n=1,2,3,...
• Si pour un certain x ∈ B(z_n,d_n×r) alors on colorie le pixel en bleu.

Ce dernier test se rajoute aux tests précédents (seuil et bassin d'attraction).

Ci-dessous, on a choisi la pointe de l'ensemble de Julia, autrement dit le point fixe réel répulsif x=2.357... (j'ai calculé cette approximation à part), qui est aussi le point le plus à gauche de J sur l'axe réel. Le programme implémente la méthode décrite ci-dessus avec ce choix, les instructions correspondantes sont en rouge.

Il y a d'autres ajouts, signalés par d'autres couleurs. En bleu, un détecteur d'overflow par le mécanisme "try/catch" de capture d'exceptions. J'ai ajouté un compteur pour connaître le nombre de pixels qui ont subi un overflow et un autre, en vert, pour ceux pour lesquels le nombre max d'itérations est atteint. Enfin en magenta une variable sur pour gérer plus facilement le sur-échantillonnage.

from __future__ import division
import png
from array import array
import math
import cmath

center = 5 + 0j
xw = 7.0
sur = 1
width = sur*600
height = sur*400
max_iter = 50
thick = 1.0

tip=2.35767667
seuil_bassin=1.
seuil=500.

tau=2*math.pi

pix_size = xw/width

NONE=0 ; MAX=1 ; BASIN=2 ; SEUIL=3 ; OVERFLOW=4 ;  TIP=5

def N2(z) : return z.real**2 + z.imag**2

def compute(j,i) :
  z = center + (j-width/2+(height/2-i)*1.0j)*pix_size
  k = 0
  reason = NONE
  r = thick*pix_size
  r_sq = r*r
  global nb_overflow,nb_max
  
  while True:
    if k>=max_iter : reason = MAX; break
    
    z = z - tau*1j*math.floor(0.5+z.imag/tau)
    
    d_sq = N2(z-tip)
    if d_sq<r_sq : reason = TIP; break
    
    if z.real < seuil_bassin : reason = BASIN; break
    
    if z.real > seuil : reason = SEUIL; break;

    try : 
      z = cmath.exp(z-1.5)
      r_sq = r_sq * N2(z)
    except OverflowError: reason=OVERFLOW; break
    
    k += 1
    
  if reason==BASIN :    return [255,255,255]
  if reason==SEUIL :    return [255,127,255]
  if reason==TIP :      return [0,0,255]
  if reason==OVERFLOW : nb_overflow += 1 ; return [0,255,0]
  nb_max +=1 ; return [255,0,0]



pixels=[]

nb_overflow = 0
nb_max = 0
for i in xrange(height) :
  row = []
  for j in xrange(width) :
    row.extend(compute(j,i))
  pixels.append(row)
print 'nb_overflow',nb_overflow
print 'nb_max',nb_max


image = png.from_array(pixels,'RGB')
image.save('essai-3.png')

Comme on peut le voir j'ai décidé de travailler avec le carré du rayon plutôt que le rayon. Quand j'ai commencé la programmation il y a 30 ans, la fonction racine carrée était coûteuse en temps de calcul, j'ai donc pris le réflexe de l'éviter. J'ignore si c'est encore d'actualité. Notez aussi que j'ai dû abaisser le seuil. Dans ces conditions j'obtiens nb_overflow=0 et nb_max=0 sur tous les exemples.

Voici quelques exemples d'images produites

Dans les images suivantes, j'ai mis le mauve en bleu.

On voit une transition. Cela s'explique ainsi : nous représentons un épaississement de J. Pas exactement l'ensemble V_ε des points situés à epsilon de J, mais un ensemble qui contient un V_ε et est contenu dans un V_ε' avec ε et ε' pas trop loins. Quand les doigts d'une génération donnée sont suffisamment serrés, leurs ε-voisinnages respectifs se touchent et forment une surface en continuité. A partir de ce moment il n'y a plus de différence visible entre les différent niveaux de resserrement.

La transition n'est pas aussi nette que celle qui était due au seuil dans nos expériences précédentes, ce qui est une bonne chose à mon sens. Mais elle reste rapide. La couleur moyenne dévie du bleu moyen quand l'espacement inter digital devient comparable à ε. Sur la 2e génération de doigts, cet espacement = constante/exp(Re(z)). Il varie d'un facteur 10 pour un déplacement horizontal de ≈ 2.3, à comparer à la hauteur du doigt d'ordre 1, qui vaut ≈ 3.14. La variation s'accélère quand on va vers la droite et cela participe peut-être à la brusquerie. Ci-dessous une simulation de la couleur basée sur la formule ci-dessus, sur une bande de hauteur π :

C'est en effet un peu brusque, peut-être légèrement moins marqué que sur nos images de Julia.

Dans la dernière image on voit apparaître l'effet du seuil. Dans le programme ci-dessous, on traite le problème, voir les lignes en rouge.

[...]

NONE=0 ; MAX=1 ; BASIN=2 ; SEUIL=3 ; OVERFLOW=4 ; TIP=5 ; BANDS=6

[...]
  
  while True:
    if k>=max_iter : reason = MAX; break
        
    z = z - tau*1j*math.floor(0.5+z.imag/tau)
    
    if z.real>tip and math.log(r_sq)/2 > math.log(tau)+1.5-z.real : 
      d = z.imag / tau
      d = d-math.floor(d)
      d = math.fabs(d-0.5)
      d = max(0.0,0.25-d)
      d = d * tau
      d_sq = d*d
      if d_sq<r_sq : reason = BANDS; break
      else : reason = BASIN; break

    d_sq = N2(z-tip)
    if d_sq<r_sq : reason = TIP; break
    
[...]    
        
  if reason==BASIN :    return [255,255,255]
  if reason==SEUIL :    return [255,0,0]
  if reason==TIP :      return [0,0,255]
  if reason==BANDS :      return [0,0,255]
  if reason==OVERFLOW : nb_overflow += 1 ; return [0,255,0]
  nb_max +=1 ; return [255,0,0]

[...]

Pour ce test, on idéalise le Julia comme une alternance de bandes bleues et blanches de hauteur π. On regarde le disque B(z_n,d_n×r) et on teste s'il rencontre une bande bleue. Si oui, pixel←bleu, sinon pixel←blanc. Mais on ne fait le test que si la taille du disque est suffisamment grande pour que l'idéalisation soit raisonnable : ici j'ai pris que d_n×r doit être plus grand que la distance approximative entre les doigts de la génération 2. Résultat :

C'est intéressant : la transition se fait plus progressive qu'avec la détection des pré-images de la pointe. Par contre on récupère plus de moirés. Ci-dessous je les cache avec un "dirty trick".

---