Using an image as orbit trap for coloring a Julia set

Question

Searching the web I couldn't find any example on how to use a bitmap (a raster image) as orbit trap to color a Julia set or Mandelbrot set. Here Inigo Quilez is explaining the method he uses in only two lines, and displays the resulting Julia set. Also Malin Christersson displays the upper-left image, and on another page a Julia set based on Mona Lisa's image, without any comment or explanation. Hence I need help in undersanding how orbit trap method works when the trap is an image. How is placed the image in the complex plane to detect the orbit trapping? I succeeded only to plot the Julia set when I defined a meshgrid over a rectangular region in the complex plane, of the same resolution like the image, and iterated the points of the meshgrid until abs(z)>2 or the number of iterations exceeds a threshold. The Julia set gets at the position of the stating point the color of the image pixel, corresponding to the last point of the computed orbit. Here is what I got, using Julia lang:

using Images, Interpolations

f(z::Complex; c= 0.56666-0.5im, n=4) =  z^n + c
function iteratef(z, iter)
    k=0
    while abs(z) <2.0 && k < iter
        z=f(z;)
        k+=1
    end
    return z
end
function JuliasetImage(img::Union{Matrix{RGB{T}}, Matrix{Gray{T}}}; 
                rectangle=(remin=-1.5, remax=1.5, immin=-1.5, immax=1.5), iter=255) where T
  
   #define  a constant  image for fractal, inialized on white
    frimg = fill(RGB{N0f8}(1, 1, 1), size(img)...)
    nrow, ncol = size(img)
    a, b, c, d = rectangle
    @assert a<b && c<d || error("wrong boundaries for the rectangle {$a, $b}x[$c, $d]")
    for j in axes(img, 2)
        for i in  axes(img, 1)
            z = Complex(a+(b-a)*(j-1)/(ncol-1), d-(d-c)*(i-1)/(nrow-1)) 
            z = iteratef(z, iter)
            #get indices in the image of the last z
            idxr = floor(Int, 1.5+(nrow-1)*(d-imag(z))/(d-c))
            idxc = floor(Int, 1.5+(ncol-1)*(real(z)-a)/(b-a))
            if idxr < 1 || idxr > nrow
               idxr = mod(idxr, nrow) + 1  
            end
            if idxc < 1 || idxc > ncol
               idxc = mod(idxc, ncol) + 1
            end
          
            frimg[i,j] = img[idxr, idxc]
        end
   end
   return frimg
end

img =load("ai-boy.jpg")
frimg = JuliasetImage(img;)
fritp = interpolate(frimg, BSpline(Linear()), OnGrid())

But I want to understand how the above mentioned Julia sets have been generated based on an image.

Answer 1

Finally I succeeded myself to plot a Julia set with an image as orbit trap. The algorithm for detecting the trapped points of an orbit works as follows:

The Julia fractal is included in a rectangular region, [A, B] x [C,D], from the complex plane. Over this rectangle is defined a meshgrid of resolution Nrow x Ncol. The trapping set is another rectangle, [a,b] x [c,d], included in the first one, chosen by user through experimentation. On the latter rectangle one defines a meshgrid of the same resolution as the image resolution. Thus is ensured a one-to-ne correspondence between the points of the last meshgrid and the image pixels.

The algorithm works similar to the known one, for plotting a Julia set. But here is tested for each point of an orbit whether it entered the trap (the rectangle [a,b]x[c,d]). If this is a case, and (i,j) are indices of the nearest point to the trapped point, in the meshgrid over this rectangle, then the fractal image at the position, (k,l), of the orbit starting point, is colored with that pixel color. The code is written in a non-Julian way to be accessible to all those interested in this type of orbit trap, even if they do not use Julia lang. This code works both for images with background removed (i.e. png images with background pixels of color RGBA(0,0,0,0)) as well as jpg images with their original background. If the background is removed then, when the orbit is in fact trapped, but the point of the orbit in the trap [a,b] x[c,d] corresponds to a transparent pixel (i.e. it is not a pixel on the shape represented by the image) it is considered as untrapped (see the function istrapped). Here is an example with an image of transparent background, set via fotor.

using Images, Interpolations

function get_index(z::Complex, nrow::Int, ncol::Int; imgrectangle=(a=-0.25, b=0.7, c=-0.25, d=0.7))
    a, b, c, d = imgrectangle
    idxrow = floor(Int, 1.5+ (nrow-1)*(d-imag(z))/(d-c)) #1.5 =1.0+0.5; 1
    idxcol = floor(Int, 1.5+ (ncol-1)*(real(z)-a)/(b-a))
    return idxrow, idxcol
end  

function istrapped(img::Matrix{T}, i::Int, j::Int) where T
       1 <= i <=size(img, 1)  && 1 <= j <= size(img, 2) && (img[i, j] != RGBA(0,0,0,0))     
end 

f(z::Complex; c=-0.8+0.156im)=z^2+c 

function iteratef(z::Complex, img::Matrix{T}; maxiter=512) where T
    n=0
    i, j = get_index(z, size(img)...;)
    while (n < maxiter && abs(z) < 2) && (n < 2 || !istrapped(img, i, j ))
        z = f(z)
        n += 1
        i, j = get_index(z, size(img)...;)
        if istrapped(img, i, j)
            return img[i, j]
        end
    end
end    

function Juliafractal(img::Matrix{T}; 
                      fractrectangle=(xmin=-1.55, xmax=1.55, ymin=-1.4, ymax=1.4),
                      Nrow=1000, Ncol=1000,  maxiter=1024,
                      bgcolor=RGB{N0f8}(1.0,1.0, 1.0)) where T
    A, B, C, D = fractrectangle
    frimg = fill(bgcolor, 1:Nrow, 1:Ncol)
    for l in axes(frimg, 2)
        for k in  axes(frimg, 1)
            z = Complex(A+(B-A)*(l-1)/(Ncol-1), D-(D-C)*(k-1)/(Nrow-1)) 
            vr = iteratef(z, img; maxiter=maxiter) 
            if vr != nothing
                frimg[k, l] = vr
            end    
        end
    end
    return interpolate(frimg, BSpline(Linear()), OnGrid())
end    

img = load("imgs/ai-younglady.png") 
fritp = Juliafractal(img; bgcolor=RGB{N0f8}(1,0.953,0.91))