Question

Como gerar o Vicsek fractal usando contrações?

Answer 1

Seguindo a ideia do Triângulo de Sierpinski descrita aqui, o Vicsek fractal pode ser gerado através de uma contração \(V\), que pode ser definida como a união dos seguintes operadores:

\[V_1(x,y)=\frac{1}{3}\left[\begin{array}{c} x \\ y \\ \end{array}\right]\]

\[V_2(x,y)=\frac{1}{3}\left[\begin{array}{c} x \\ y \\ \end{array}\right]+\left[\begin{array}{c} 0 \\ \frac{2}{3} \\ \end{array}\right]\]

\[V_3(x,y)=\frac{1}{3}\left[\begin{array}{c} x \\ y \\ \end{array}\right]+\left[\begin{array}{c} \frac{2}{3} \\ 0 \\ \end{array}\right]\]

\[V_4(x,y)=\frac{1}{3}\left[\begin{array}{c} x \\ y \\ \end{array}\right]+\left[\begin{array}{c} \frac{1}{3} \\ \frac{1}{3} \\ \end{array}\right]\]

\[V_5(x,y)=\frac{1}{3}\left[\begin{array}{c} x \\ y \\ \end{array}\right]+\left[\begin{array}{c} \frac{2}{3} \\ \frac{2}{3} \\ \end{array}\right]\]

Podemos gerar a imagem com o número de iterações desejadas através da seguinte implementação em python:

import numpy as np
import matplotlib.pyplot as plt

def generatePlane(N):
    plane=np.empty([N*N,2]) 
    count=-1
    for i in range(N):
        for j in range(N):
            count=count+1
            plane[count,0]=1.0*i/N
            plane[count,1]=1.0*j/N
    return plane

def contractionVicsek(plane):
    [m,n]=np.shape(plane)
    newRegion=np.empty([5*m,2])
    count=-1
    for i in range(m):
        count=count+1
        newRegion[count,0]=(1.0*plane[i,0])/3.0 
        newRegion[count,1]=(1.0*plane[i,1])/3.0 

    for i in range(m):
        count=count+1
        newRegion[count,0]=(1.0*plane[i,0])/3.0
        newRegion[count,1]=(1.0*plane[i,1])/3.0+2.0/3.0 

    for i in range(m):
        count=count+1
        newRegion[count,0]=(1.0*plane[i,0])/3.0+2.0/3.0
        newRegion[count,1]=(1.0*plane[i,1])/3.0

    for i in range(m):
        count=count+1
        newRegion[count,0]=(1.0*plane[i,0])/3.0+(1.0)/3.0
        newRegion[count,1]=(1.0*plane[i,1])/3.0+(1.0)/3.0

    for i in range(m):
        count=count+1
        newRegion[count,0]=(1.0*plane[i,0])/3.0+(2.0)/3.0
        newRegion[count,1]=(1.0*plane[i,1])/3.0+(2.0)/3.0
    return newRegion    

if __name__ == '__main__':

    N=100
    plane=generatePlane(N)

    #iteration 1
    region=contractionVicsek(plane)
    #iteration 2
    region=contractionVicsek(region)
    #iteration 3
    region=contractionVicsek(region)
    #iteration 4    
    region=contractionVicsek(region)  

    fig = plt.figure()
    ax = fig.add_subplot(111)  
    fig.hold()
    ax.plot(region[:,0],region[:,1],'r',markersize=10, markeredgewidth=0)
    ax.set_ylabel('y')
    ax.set_xlabel('x')

Após 4 iterações, temos: