IFS movie maker

Youtube video here – NB: select 720p HD quality

Each time it’s run, the program will make frames for a different movie; with current settings, 5000 frames for a 15fps movie.

Algorithm: Starting from a random point, at each step one of the 3 maps is randomly chosen.
++++Each map is:

    \begin{equation*} \left\{ \begin{align*} &&&\quad x_{n+1} = ax_n\cos{\theta} - by_n\sin{\theta} + c \\ &&&\quad y_{n+1} = dx_n\sin{\theta} + ey_n\cos{\theta} + f \end{align*} \right. \end{equation*}

where a, b, d and e range between 0.2 and 1.1. It’s identical to the basic Sierpinski triangle IFS except the 3 maps are randomly changing position (c,f), x and y shrink (a, b, d, e), and angle (\theta).
++++The maps are assigned the colours R, G, B. The pixel colour encodes the previous 8 maps that were chosen: The latest map is 128, the one before 64, the one before 32 etc. So that if e.g. the 8 last maps used were all the red one, the pixel will be coloured RGB=255,0,0. Up to 10 colours are stored for each pixel, which are averaged at the end to find the colour for the image.
If the point goes off screen, adjust variables and start again. Continue until either 1000 points have been calculated with nothing drawn (just revisiting points) or 10 million points plotted, then save image.

I made a 15 frames/sec, h264 video with an old version of ffmpeg, with crf 15. (crf 23 looked fine; crf 15 makes about a 50% larger file.)

Cython program IFSmoviev9.pyx:
On the Sage command line %attach IFSmoviev9.pyx then OK() to run it.

#cython: boundscheck=False, wraparound=False, cdivision=True, nonecheck=False
#     IFS movie maker......v9 - june 2017
#ffmpeg command to make 15 fps movie from images:
#ffmpeg -f image2 -r 15 -i ~/ifs/0%04d.png -vcodec h264 -crf 23 ifsv9.mp4
 
import numpy as np
cimport numpy as np
from scipy import misc
from libc.stdlib cimport rand
from libc.math cimport sin,cos,fabs,sqrt
 
DEF EQNS = 2 #number of equations i.e. currently x and y -> 2 dimensions
DEF COEFFS = 2 #i.e. a, b (=coef[0],coef[1].)
DEF COLOURS = 3 #i.e. R,G,B currently
DEF SCREENWID = 1280
DEF SCREENHEI = 720
DEF PI = 3.14159265
DEF MAXLOOPS = 10000000
DEF MARGIN = 30 #20 #50
DEF REV = -0.2 #reverse: multiply speed by this when edge of screen/min/max reached
DEF MAXCONSTV = 0.5 #the most the C terms are allowed to change each time.
DEF MAXCONSTA = 0.1 #max translation acceleration , i.e. velocity difference
DEF SHRINKSTEP = 0.001 #amount largest shrink is reduced by if image is too big.
DEF REBOUNDRATIO = 0.005 # as if a car travelling 100m/s hit a wall and rebounds at 0.5m/s
DEF MAXSHRINKV = 0.05
DEF MAXSHRINKA = 0.01
DEF MAXANGV = 2*PI/(10*15) # at 15fps, max turn speed = 360 degrees in 10 seconds
DEF MAXANGA = MAXANGV/5
DEF MAXANGJERK = MAXANGA/5 #measure in radians/second^3
DEF GAPINC = 0.3 #if there's a gap on any side, move closest map offset towards edge
 
# Constants to be adjusted more often than those above:
DEF MAPS = 3 #12 #6
DEF FRAMES = 5000 #8000
DEF MAXSHRINK = 1.1 #1.2 is good for 3 maps. but too high for 12, it seems.
DEF MINSHRINK = 0.2
DEF LOOPSSINCEANYTHINGDRAWNLIMIT = 1000 #maybe 100 is enough?
#frame finished when the last this-many points were already drawn on at least MAXPTSPERPIXEL times.
DEF MAXPTSPERPIXEL = 10 #max number of points averaged to find col of each pixel
DEF MOVEINDIVISOR = 500 #200
DEF PICFOLDER = "/Users/admin/ifs/" #NB create this folder before running program
DEF LOGFILE = "/Users/admin/ifs-log.txt"
#each map has 2 equations, basically this:
#eq[0]: x = ax+by+c : a=.coef[0],b=.coef[1],c=.const
#eq[1]: y = dx+ey+f : d=.coef[0],e=.coef[1],f=.const
#with an added rotation term, making :
#eq[0]: x = axcos(ang) - bysin(ang) + c
#eq[1]: y = dxsin(ang) + eycos(ang) + f
cdef:
    struct eqstruct:
        float const, constv, coef[COEFFS], v[COEFFS]
    struct mapstruct:
        eqstruct eq[EQNS]
        float angle, angularv, angulara
    mapstruct map[MAPS]
    struct dirs:
        int left, right, top, bottom 
    dirs mostPt #i.e. x-coord of leftmost point in image etc
    struct px:
        short rgb[3], pts # = num of points' data added together in that pixel
    px image[SCREENWID][SCREENHEI] 
    int loopsdone, loopsSinceAnythingDrawn
    short out[2], colpowers[8]
    float pt[2],p[MAPS],ptotal,pcumulative[MAPS]
    enum:
        ALLGOOD = 0, WENTOFFSCREEN = 1, PLUS = 1, MINUS = -1
        InXEqn = 0, InYEqn = 1, X = 0,Y = 1
colpowers[:] = [128,64,32,16,8,4,2,1]
txt = open(LOGFILE, "w")
 
cdef SetInitialValues():
    cdef short m,e,t
    cdef float me[COLOURS][2][4]
    for m in xrange(MAPS):
        map[m].angle = StepRand(PI) 
        map[m].angularv = StepRand(MAXANGV)
        map[m].angulara = StepRand(MAXANGA)
    #initial fractal= a rough tiny sierpinski tri, doesnt really matter what it is.
    me[0][0][:] = [-40,1,0.5,0]
    me[0][1][:] = [-20,-1,0,0.5]
    me[1][0][:] = [40,1,0.5,0]
    me[1][1][:] = [-20,-1,0,0.5]
    me[2][0][:] = [0,0,0.5,0]
    me[2][1][:] = [30,1,0,0.5]
    for m in xrange(MAPS):
        for e in xrange(EQNS):
            map[m].eq[e].const = me[m%COLOURS][e][0] #i.e. puts the same values in map 1,4,7,10 etc etc
            map[m].eq[e].constv = me[m%COLOURS][e][1]
            map[m].eq[e].coef[:] = [me[m%COLOURS][e][2],me[m%COLOURS][e][3]]
            map[m].eq[e].v[:] = [0.01,0.01]
 
cdef float StepRand(float max): #returns num between -max and max. step size=that range divided by 200
    return -max + (rand() % 401)*max/200
 
cdef float fsign(float num): #sign of floating pt var: 1, -1 or 0
    if num < 0:
        return -1
    else:
        return num > 0
 
cdef KeepVariablesWithinLimits():
    cdef:
        short m,e,t
        float s
        eqstruct *ME
    for m in xrange(MAPS):
        for e in xrange(EQNS):
            map[m].angle = map[m].angle%(2*PI) #keep angle value in -PI< <PI but dont change angle!!!
            s = fsign(map[m].angularv)
            if s*map[m].angularv > MAXANGV:
                map[m].angulara = -s*MAXANGA*REBOUNDRATIO
            ME = &map[m].eq[e]
            s = fsign(ME.constv)
            if s*ME.constv > MAXCONSTV:
                ME.constv = s*MAXCONSTV 
            for t in xrange(COEFFS): #0,1 , i.e. x and y terms, not c
                s=fsign(ME.coef[t])
                if s*ME.coef[t]>MAXSHRINK: #if >max or <-max
                    ME.coef[t] = s*MAXSHRINK
                    ME.v[t] = -s*REBOUNDRATIO
                if fabs(ME.v[t]) > MAXSHRINKV:
                    ME.v[t] = fsign(ME.v[t])*MAXSHRINKV 
                if s*ME.coef[t]<MINSHRINK: #if <min and >-min
                    ME.coef[t] = s*MINSHRINK
                    ME.v[t] = s*REBOUNDRATIO
 
cdef UpdateEquationValues():
    cdef short m,e,t
    for m in xrange(MAPS):
        map[m].angularv += map[m].angulara
        map[m].angle += map[m].angularv
        for e in xrange(EQNS):
            map[m].eq[e].const += map[m].eq[e].constv
            for t in xrange(COEFFS):
                map[m].eq[e].coef[t] += map[m].eq[e].v[t]
 
cdef CalculateMapChoiceWeights():
    global ptotal
    cdef short i
    ptotal = 0
    #the formula: using x=ax+by, y=dx+ey (ignoring constant and rotation)
    #the image points of the triangle (0,1), (0,0) and (1,0), are
    #(b,e), (0,0) and (a,d). The area of this triangle is |bd-ae|/2.
    #Since only the ratios between maps are important, just use |bd-ae|
    for i in xrange(MAPS):
        p[i] = fabs(map[i].eq[0].coef[1]*map[i].eq[1].coef[0]-map[i].eq[0].coef[0]*map[i].eq[1].coef[1])
        ptotal += p[i]
        pcumulative[i] = ptotal
     
cdef short ChooseAMap():
    cdef short i,ir = rand() % <int>(ptotal*10000+1)
    cdef float r = (<float>ir)/10000
    for i in xrange(MAPS): 
        if r <= pcumulative[i]:
            return i #i think but not sure that this ALWAYS works. otherwise:
    return MAPS-1
 
cdef CalculateNewPoint(short m):
    cdef float tempx, cosa = cos(map[m].angle), sina = sin(map[m].angle)
    tempx = pt[X]*cosa*map[m].eq[0].coef[0] - pt[Y]*sina*map[m].eq[0].coef[1] + map[m].eq[0].const 
    pt[Y] = pt[X]*sina*map[m].eq[1].coef[0] + pt[Y]*cosa*map[m].eq[1].coef[1] + map[m].eq[1].const
    pt[X] = tempx
 
cdef short TryMillionsOfLoops(): # it's a different num of loops every time.
    global loopsdone, loopsSinceAnythingDrawn
    cdef:
        short col,i,intpt[2],m, last8[8], RGB[3]
        int count
    loopsSinceAnythingDrawn = 0
    while loopsSinceAnythingDrawn < LOOPSSINCEANYTHINGDRAWNLIMIT: 
        count += 1
        if count > MAXLOOPS: 
            break
        m = ChooseAMap()
        CalculateNewPoint(m)
        intpt[:] = [<short> pt[X] + SCREENWID/2 , <short> pt[Y] + SCREENHEI/2]
        for i from 6 >= i >= 0: 
            last8[i+1] = last8[i]
        last8[0] = m%COLOURS
        if count < 10: 
            continue #dont draw first 10 values
        if intpt[X] >= SCREENWID or intpt[X] < 0 or intpt[Y] >= SCREENHEI or intpt[Y] < 0:
            out[:] = intpt
            return WENTOFFSCREEN #things arent All Good.
        if image[intpt[X]][intpt[Y]].pts < MAXPTSPERPIXEL:
            RGB[:] = [0,0,0]
            for i in xrange(8):
                RGB[last8[i]] += colpowers[i] #0->128,1->64,2->32,3->16,4->8,5->4,6->2,7->1
            for col in xrange(3):
                image[intpt[X]][intpt[Y]].rgb[col] += RGB[col]
            image[intpt[X]][intpt[Y]].pts += 1
            loopsSinceAnythingDrawn = 0 #reset counter as point was drawn to
        else: loopsSinceAnythingDrawn += 1
        if intpt[X] > mostPt.right: mostPt.right = intpt[X]
        if intpt[X] < mostPt.left: mostPt.left = intpt[X]
        if intpt[Y] > mostPt.top: mostPt.top = intpt[Y]
        if intpt[Y] < mostPt.bottom: mostPt.bottom = intpt[Y]
    loopsdone = count #save count to write to logfile later
    return ALLGOOD
 
cdef short MapWithGreatestConst(short equ): #in equations equ (i.e. in X eqs or Y eqs)
    cdef short m,greatest = 0
    cdef float maxval = map[0].eq[equ].const #assume map[0] is greatest..
    for m in xrange(1,MAPS):
        if map[m].eq[equ].const > maxval:
            greatest = m
            maxval = map[m].eq[equ].const
    return greatest
 
cdef short MapWithSmallestConst(short equ):
    cdef short m,smallest = 0
    cdef float minval = map[0].eq[equ].const #assume map[0] is least..
    for m in xrange(1,MAPS):
        if map[m].eq[equ].const < minval:
            smallest = m
            minval = map[m].eq[equ].const
    return smallest
 
cdef AddToGreatestConst(float ch, short equ):
    map[MapWithGreatestConst(equ)].eq[equ].const += ch
 
cdef AddToSmallestConst(float ch,short equ):
    map[MapWithSmallestConst(equ)].eq[equ].const += ch
 
cdef MoveConstsInward(float xr, float yr): 
    cdef short m,e,t
    if out[X] >= SCREENWID: #i.e. if 100, move in by 0.5; if 400, move in by 2
        AddToGreatestConst(-xr/MOVEINDIVISOR,InXEqn) 
    if out[X] < 0:
        AddToSmallestConst(xr/MOVEINDIVISOR,InXEqn) 
    if out[Y] >= SCREENHEI:
        AddToGreatestConst(-yr/MOVEINDIVISOR,InYEqn)
    if out[Y] < 0:
        AddToSmallestConst(yr/MOVEINDIVISOR,InYEqn)
         
cdef ReduceLargestLinearCoefficient(): #find max shrink val..
    cdef float maxshrinkval = 0, sign
    cdef short maxterm,m,e,t,maxm,maxe
    for m in xrange(MAPS):
        for e in xrange(EQNS):
            for t in xrange(COEFFS):
                if fabs(map[m].eq[e].coef[t]) > maxshrinkval:
                    maxm = m
                    maxe = e
                    maxshrinkval = fabs(map[m].eq[e].coef[t])
                    maxterm = t
    sign = fsign(map[maxm].eq[maxe].coef[maxterm])
    map[maxm].eq[maxe].coef[maxterm] -= sign*SHRINKSTEP
 
cdef float ConstRange(short equ):
    return map[MapWithGreatestConst(equ)].eq[equ].const-map[MapWithSmallestConst(equ)].eq[equ].const
 
cdef MakeFractalSmaller(): 
    cdef short ConstsAreDecentlySpreadOut
    cdef float xrange = ConstRange(X), yrange = ConstRange(Y)
    ConstsAreDecentlySpreadOut = xrange*yrange > 10000 #hmm 100x100, not much..
    print "xrange*yrange: %f x %f = %f ? %d" % (xrange,yrange,xrange*yrange,ConstsAreDecentlySpreadOut)
    if ConstsAreDecentlySpreadOut:
        if rand()%2: # 50% chance
            MoveConstsInward(xrange,yrange)
            return
    ReduceLargestLinearCoefficient()
 
cdef ShiftAllXorYConstants(short XorY, float sign):
    cdef short i
    for i in xrange(MAPS):  
        map[i].eq[XorY].const += sign*GAPINC #lil bit
     
cdef MakeFractalMoreCentred(): 
    cdef short i, LeftGap = mostPt.left, RightGap = SCREENWID-mostPt.right
    cdef short TopGap = SCREENHEI-mostPt.top, BottomGap = mostPt.bottom
    #if a bigger (by > MARGIN) gap on an edge, move nearest map const closer to that edge
    if LeftGap > MARGIN+RightGap:
        ShiftAllXorYConstants(InXEqn,MINUS)
    elif RightGap > MARGIN+LeftGap:
        ShiftAllXorYConstants(InXEqn,PLUS)
    if TopGap > MARGIN+BottomGap:
        ShiftAllXorYConstants(InYEqn,PLUS)
    elif BottomGap > MARGIN+TopGap:
        ShiftAllXorYConstants(InYEqn,MINUS)
 
cdef RandomlyAlterOneValueInEachEquation():
    cdef short m,e,t
    for m in xrange(MAPS):
        map[m].angulara += StepRand(MAXANGJERK) 
        for e in xrange(EQNS):
            t = rand()%(COEFFS+1)
            if t == COEFFS: #ie 2
                map[m].eq[e].constv += StepRand(MAXCONSTA) 
                continue #else if t == 0 or t == 1:
            map[m].eq[e].v[t] += StepRand(MAXSHRINKA)
 
cdef WriteToLogFile(short n):
    cdef short m
    txt.write("Frame no. %d - loops: %d\n" % (n,loopsdone))
    for m in xrange(MAPS):
        txt.write("Map %d        a          b          c\n" % m)
        txt.write("x: %f %f %f v: %f %f %f\n" % (map[m].eq[0].coef[0],
        map[m].eq[0].coef[1],map[m].eq[0].const,map[m].eq[0].v[0],map[m].eq[0].v[1],map[m].eq[0].constv))
        txt.write("y: %f %f %f v: %f %f %f\n" % (map[m].eq[1].coef[0],
        map[m].eq[1].coef[1],map[m].eq[1].const,map[m].eq[1].v[0],map[m].eq[1].v[1],map[m].eq[1].constv))
        txt.write("angle: %f angularv: %f angulara: %f\n" % (map[m].angle,map[m].angularv,map[m].angulara))
    txt.write("mostPt.left: %d R: %d T: %d B: %d\n" % (mostPt.left,mostPt.right,mostPt.top,mostPt.bottom))
    txt.write("ptotal=%f\n" % ptotal)
    for m in xrange(MAPS):
        txt.write("p[%d]=%f p[%d]/ptotal=%f\n" % (m,p[m],m,p[m]/ptotal))
    txt.write("-------------------------------------\n")
 
cdef DisplayOutOfBoundsData(short picnum):
    cdef short m,e,t
    print "\nOUT of BOUNDS frame not saved = ",picnum
    for m in xrange(MAPS):
        for e in xrange(EQNS):
            print "map[%d].eq[%d].const: %f v: %f" % (m,e,map[m].eq[e].const,map[m].eq[e].constv)
            print "coef[0]: %f v: %f coef[1]: %f v: %f" % (map[m].eq[e].coef[0],map[m].eq[e].v[0],map[m].eq[e].coef[1],map[m].eq[e].v[1])
     
cdef SetupMainLoopVars():
    cdef short i, j
    for i in xrange(SCREENWID):
        for j in xrange(SCREENHEI):
            image[i][j].pts = 0
            image[i][j].rgb[:] = [0,0,0]
    pt[:] = [SCREENWID/2,SCREENHEI/2] #initial point. doesnt matter what it is. 
    mostPt.left = SCREENWID
    mostPt.right = 0
    mostPt.top = 0
    mostPt.bottom = SCREENHEI
    CalculateMapChoiceWeights()
 
cpdef OK():
    cdef:
        np.ndarray[np.uint8_t,ndim = 3] pic = np.empty((SCREENHEI,SCREENWID,3),dtype = np.uint8)
        short i,j,col,frame = 0, picnum = 0
    SetInitialValues()
    while True: #done like this because not every loop makes an image.
        SetupMainLoopVars()
        if TryMillionsOfLoops() == WENTOFFSCREEN: #then dont save
            MakeFractalSmaller()
            DisplayOutOfBoundsData(picnum)
            continue #start again; try another IFS with new vals
        for i in xrange(SCREENWID):
            for j in xrange(SCREENHEI):
                if image[i][j].pts > 0: #number of pts stored 'in' that pixel
                    for col in xrange(3):
                        pic[SCREENHEI-1-j,i,col] = image[i][j].rgb[col]/image[i][j].pts #swaps x and y..
                else: #not used, make black.
                    for col in xrange(3):
                        pic[SCREENHEI-1-j,i,col] = 0   
        name = '{}{:05}.png'.format(PICFOLDER,picnum)
        misc.imsave(name,pic)
        print picnum
        WriteToLogFile(picnum)
        picnum += 1
        if picnum > FRAMES:
            break
        MakeFractalMoreCentred() 
        RandomlyAlterOneValueInEachEquation() 
        UpdateEquationValues()
        KeepVariablesWithinLimits()
    txt.close()

Program notes

1. None of the ‘shrink’ values (a,b,d,e) are ever negative – making a flipped-over fractal – which is a shame. I couldn’t think of a good way of doing that, as they can’t get too close to 0. Except starting some of them off negative, maybe a good idea. I think the program will correctly deal with neg values as is, i.e. they will stay negative, between -0.2 and -1.1 with current settings.
2. Things are added to the angular acceleration, but to the velocity of the shrink (a,b,d,e) and location (c and f). There’s no reason for that, I just didn’t get around to adding jerk (rate of change of acceleration) to the others.
3. I tried with more than 3 maps, 12 also works well. (A minimum shrink of 0.1 seems better with large numbers of maps) Any number can be used >=3. The colours R,G,B are re-used for colours mod 3. Maybe 2 maps works also, I didn’t try that.
4. Things move too fast if anything; it should be slowed down a bit. Slowed down a lot more, 25 fps movies would be better than 15 I guess, but 15 looks fine.
p.s. I was reading Refactoring and Clean Code (both wonderful programming books) while working on this; I hope it shows. 🙂