blender_fourier_epicycloids.py

import bpy
import math
import numpy as np
import random

rawdata = ['978.1466642174169,416.0576927689958',
 '1006.8180157843148,417.5706213803142',
 '1036.7920721872522,416.6868738547563',
 '1066.747012249827,416.38719890540835',
 '1059.6128569614289,409.33416276234385',
 '1030.4718652642766,409.628854704408',
 '1008.5275919584764,420.25884165882167',
 '1036.4203684584786,429.65852805336385',
 '1064.6695258001453,421.36596761061674',
 '1080.8564753191515,425.9346107654944',
 '1090.9984963426,400.82704344706326',
 '1078.6661209658157,375.1664727697521',
 '1064.0144585072057,350.25253770164574',
 '1059.8475661455668,323.93839257157595',
 '1087.1871874100893,332.6397822206011',
 '1105.1059792999029,316.14430017846826',
 '1077.175046723591,318.1169795123767',
 '1053.703865537524,306.48574862718584',
 '1058.225717806081,278.0421617637622',
 '1084.452139150267,264.66398768225235',
 '1108.6607852363586,279.80636041975026',
 '1109.9401065979005,302.74444876098636',
 '1085.2504394566724,289.73662123329007',
 '1058.0552730576014,298.2836563703446',
 '1068.7866844735145,308.06452259159084',
 '1088.8364804506116,301.6318118687812',
 '1115.253388076255,315.59369279529704',
 '1128.1572909822555,339.52934132504276',
 '1146.5404981807806,350.0051269094656',
 '1162.269709075033,328.3426062877566',
 '1155.3540038117635,299.19646220953473',
 '1173.4451232142446,278.82681166744226',
 '1161.044548750053,278.2542499944765',
 '1167.7861322234614,260.84601954238457',
 '1172.5363431703317,231.60312845227486',
 '1159.3293471098475,213.628257404038',
 '1164.9530392292977,185.1047985041461',
 '1153.4034589106652,157.70456697254528',
 '1129.515820828623,144.09878957001018',
 '1127.726126945944,123.30631653741972',
 '1110.4286695885544,99.9447931205676',
 '1087.3814076555084,90.25549114093583',
 '1069.4235962551095,68.37482075941926',
 '1041.1459798879305,74.45901309162323',
 '1031.5595713938603,58.062273905340504',
 '1020.2467463603627,49.0748725717776',
 '1004.201848598177,72.78956949752504',
 '977.0443885996532,63.36270830252842',
 '949.7206862542814,52.19057912174605',
 '923.3768584249019,63.76739849805824',
 '896.1574551459275,64.52551336386742',
 '877.3427119991331,86.07753105109177',
 '849.9608565829661,95.01036628087789',
 '831.5292999222204,117.6655084491073',
 '819.3654740803756,142.10357366468435',
 '802.6985884536655,166.64025882287245',
 '799.9778804291485,196.3460626153945',
 '789.433805785656,222.04187707757941',
 '789.7298593412536,250.6460222195745',
 '792.1914666050973,271.5166639492116',
 '787.5311748166084,298.6722774066924',
 '805.3605669977735,322.30568913559836',
 '806.4908128194212,350.7519532231538',
 '825.1992916277546,372.827607452347',
 '843.7861779785155,392.9475413513182',
 '822.0227207198853,409.9479972100096',
 '815.3921139655857,438.64238140839325',
 '808.9462128678133,467.9378576473391',
 '801.2797156945089,496.93750625711397',
 '787.0451242002493,520.1851481285263',
 '768.3939242641702,543.0677247453955',
 '749.2845157683714,562.8526668527126',
 '725.1938734938783,580.6647380489387',
 '703.0184260047147,600.8583110166996',
 '682.2296494620856,622.4382418626228',
 '658.9400968081057,641.3547778175473',
 '637.0055493552186,661.7467868168475',
 '620.5959846064281,686.833495412328',
 '606.539741474095,713.3337581527118',
 '594.4100699803955,740.7466408093295',
 '584.8749225317647,769.1883434039999',
 '579.7286985848407,798.6938797748791',
 '578.9505205461227,828.672042019606',
 '578.1262376951071,858.673749615512',
 '577.2986430512902,888.6623793582785',
 '576.6551753537049,918.6391299156733',
 '587.473943231686,941.8244633377129',
 '616.7358682844088,948.3659452351558',
 '645.1367618423905,957.9961906272499',
 '673.2429771552366,968.4951939393227',
 '701.5552665954497,978.4135360149999',
 '730.5167496062887,986.1516108193387',
 '758.9768510453033,995.6039048511695',
 '787.2731046504543,1005.5697421670098',
 '811.9073972324151,999.4039262395524',
 '815.4398243729845,969.719536238946',
 '812.6366230269657,939.9231227894702',
 '807.7807668165304,910.3236615216472',
 '806.2068280918155,880.3667098514865',
 '809.548171088636,850.6621160412429',
 '816.2070622788435,821.4146102250675',
 '809.5090475548581,792.8678302171899',
 '812.0068851659096,787.1557432869096',
 '821.9285871819856,815.1191455785614',
 '822.0920631469473,845.0683471306743',
 '818.1123924767971,874.7958788846777',
 '817.719119416261,904.6817120974295',
 '819.4755451943281,934.6196118792706',
 '819.9447129281169,964.607347139099',
 '817.8424557461474,994.5254579437762',
 '822.0385615969265,1022.6489192781614',
 '851.0883332872866,1029.8073886880286',
 '880.83610328016,1033.3853046555696',
 '910.8143468565403,1034.5825355757568',
 '940.7980449034123,1035.6732564662057',
 '970.7380048568921,1037.6322580375224',
 '1000.7002438731487,1038.999413535164',
 '1030.6860939541125,1039.0656487669314',
 '1060.6541777200214,1037.5345140136637',
 '1090.202863656865,1032.6581278302303',
 '1096.718303540081,1006.1130075883602',
 '1096.253889907558,976.1524816335886',
 '1094.6543214214657,946.1933338666717',
 '1083.0500259485245,926.0973418602944',
 '1062.806431609136,945.1728705599402',
 '1037.6354104146885,935.8521463449972',
 '1042.3724077964575,908.3497085786387',
 '1069.314235396582,907.5681466480419',
 '1089.5546414288763,909.7714638534262',
 '1083.9358777827863,880.3031091487036',
 '1073.8667670495056,852.0687241924352',
 '1064.0969534725589,823.7464389971014',
 '1055.5756616163994,795.1847895495652',
 '1033.6704805926722,806.7664175894868',
 '1007.140743088504,812.2888757390776',
 '999.8622103799681,786.7617697774056',
 '1024.777975067139,775.0743997535707',
 '1047.2114994659426,778.7749465827944',
 '1036.3260611606438,751.0335618384856',
 '1023.610728301276,723.8547129474405',
 '1011.1890396095237,696.5524062290501',
 '995.7350860939029,671.089434343815',
 '981.9142215976121,687.2306264415981',
 '960.0458049964982,703.7694255805484',
 '943.4485674180839,683.3749200854528',
 '961.6521818847659,664.3189311523439',
 '985.7644201353959,653.2568320689297',
 '973.1107958943585,628.667943126064',
 '953.7950753912776,607.0625902393867',
 '945.0653364784189,578.3874206582888',
 '934.4279547405738,550.3437132786948',
 '921.8161129478019,523.1546185652871',
 '915.4306235786083,494.0480717939737',
 '889.9994350284238,478.2249963410011',
 '863.4322606583538,464.3252034274193',
 '841.125740015357,444.6430143723211',
 '832.1597089843751,416.68297338867205',
 '852.7967170038555,426.24421353474315',
 '872.8988188172264,448.3251682426311',
 '896.5020073839762,466.8387807942623',
 '919.0330250670156,486.3512391016084',
 '931.6856551446749,513.2482250151021',
 '951.3150272307547,535.9312506855878',
 '967.4292668319283,560.8602516792349',
 '976.7682332191469,589.3589217681887',
 '992.1735054567459,614.6914310812891',
 '1020.1706085609908,619.6494757937121',
 '1049.8740074390132,615.5467077661186',
 '1079.4050889415741,610.428851085663',
 '1087.6306106960706,584.3212140987808',
 '1084.8516908767651,554.45766665606',
 '1089.7600298885145,525.0700531432435',
 '1102.7138797182292,498.19748682585936',
 '1126.4459982838493,497.7994484351752',
 '1136.97622858693,525.2803775467726',
 '1141.3641615484898,554.8386718904536',
 '1149.9205216711598,583.5994019862134',
 '1157.3284258973524,612.6062397837567',
 '1163.301085467548,641.1195811544902',
 '1186.4427666210972,660.1735367151558',
 '1199.2060536511592,683.3807101490389',
 '1201.0866751629255,713.3246200883885',
 '1203.4756900589052,743.2225451550562',
 '1206.0667395629885,773.1112700500491',
 '1206.0058566746018,803.0597020789261',
 '1190.6628765297366,828.5143399579799',
 '1177.6309453532513,855.3517855548356',
 '1165.2051141925958,882.6630793922985',
 '1149.4845230351168,907.8603601462388',
 '1142.2740067132768,929.3315868414513',
 '1146.8354295628142,958.5297778352514',
 '1147.7873535919039,988.5207965548711',
 '1149.6752216796872,1018.4436347656253',
 '1171.4956373021212,1024.1400102822386',
 '1201.0585756764306,1019.0994134618817',
 '1230.1244313858747,1011.7661314844046',
 '1255.7956550998583,996.6240467788534',
 '1278.156718051425,976.6506111989333',
 '1299.3593118839642,955.4293065375681',
 '1322.9221728135228,937.0434633495938',
 '1343.4269501210638,915.664359041385',
 '1341.7949476937279,885.8054167082837',
 '1339.937637906324,855.8652814729327',
 '1338.1814581971144,825.9082622935084',
 '1335.257822617054,796.1167087140684',
 '1326.697456235868,767.3579839656383',
 '1315.3306622459254,739.6368086221582',
 '1297.3401760448858,715.9639300073924',
 '1282.4496765641816,689.9137346776718',
 '1268.1762919061382,663.5237485612813',
 '1253.932543774509,637.1300977317942',
 '1232.5350355323858,617.5281329956492',
 '1205.3635748155766,604.8125498488647',
 '1185.4678808186086,585.2842299715778',
 '1177.1723202812861,556.4485389826566',
 '1169.93178134922,527.3502407133574',
 '1160.0513274098323,500.3919462873576',
 '1135.46267647171,483.2097750968937',
 '1116.5530622079002,462.31061909971197',
 '1120.2826106155287,432.5497459239016',
 '1124.347633358909,402.82601385348767',
 '1127.9490568699782,373.0491393193195',
 '1124.0950897049142,390.4452132600552',
 '1117.224843999718,419.62816485432495',
 '1107.6098012219545,448.0404886451371',
 '1093.8969296102096,474.6720703444021',
 '1073.2518813593665,496.1344716451654',
 '1044.7638437627907,504.7891077939371',
 '1014.9621953216631,502.8689839324396',
 '986.9377571645072,492.407779566458',
 '961.303095608816,476.85404246994347',
 '935.3950598051047,461.7165069712374',
 '908.869373304675,447.70333859870885',
 '884.0187543550259,431.2291019680282',
 '874.9110931530446,403.6151221847808',
 '853.3104599395023,385.92837971987615',
 '835.8706017584964,361.6104135534689',
 '829.1969073039652,332.7435779668267',
 '850.4395607740348,327.43331758895545',
 '869.1531639092989,350.19257788805646',
 '880.2245954682978,345.896941257972',
 '888.198931383133,318.84026803016695',
 '917.3517304993569,311.90744393345744',
 '941.882932342957,297.3319863973027',
 '968.7575844986653,284.861512179537',
 '996.3267089991571,294.0739526677134',
 '988.7841366462037,317.18418587139274',
 '959.5753813782707,322.41108186734493',
 '943.4542715250524,314.04343632434615',
 '966.2426726869605,313.5361781800978',
 '987.5143134577274,306.7644904957192',
 '1005.5581814824025,296.7602005959546',
 '1006.22673846215,267.21622724056255',
 '976.9750383105352,263.7596260981189',
 '947.4503158981419,268.3470186547712',
 '920.7241653507629,281.53994285641534',
 '897.0797798722263,300.0018779252265',
 '894.2692139347633,284.5552731928911',
 '881.5430700073242,260.7740968322755',
 '872.9221121573019,232.71451288425104',
 '899.2275884654484,239.7897354509329',
 '914.4626858632707,235.12938166865453',
 '890.3336597223654,218.68960929897247',
 '904.4116647284192,216.22734997960006',
 '925.4876499567381,221.2659989079192',
 '908.996038482441,198.72253104602098',
 '928.8908680238359,176.25435178803016',
 '938.7973320875528,197.08057683531632',
 '966.656043010268,201.805375450866',
 '994.2398179731936,190.32715160444397',
 '1019.0339273331006,176.07428654910922',
 '1044.8507263183592,190.4707411499025',
 '1059.6564597088916,216.08903150109032',
 '1063.461979368081,223.63641920040726',
 '1076.318572396533,217.93066658985197',
 '1091.1897776823903,243.17447449002813',
 '1118.6781110025315,249.26156518060822',
 '1121.8660714800972,279.0835877821946',
 '1112.2895728067344,273.30200048876105',
 '1087.2438842948534,259.28878028982706',
 '1058.3004502611157,263.8314424943925',
 '1045.3938454331155,290.1723545527459',
 '1051.6618344203762,319.205366518584',
 '1060.7545923218975,347.743856915568',
 '1069.771471676087,376.22300437558823',
 '1060.3004452514645,382.2845714454651',
 '1041.007152740478,390.7117742309571',
 '1025.084280820727,377.995883095026',
 '1008.95425615833,372.42729715893233',
 '1003.6128487057681,367.8281517781615',
 '982.3058090782057,388.66547887369984',
 '970.1585135695633,416.0076674030535']

def get_grease_pencil(gpencil_obj_name='GPencil') -> bpy.types.GreasePencil:
    """
    Return the grease-pencil object with the given name. Initialize one if not already present.
    :param gpencil_obj_name: name/key of the grease pencil object in the scene
    """

    # If not present already, create grease pencil object
    if gpencil_obj_name not in bpy.context.scene.objects:
        bpy.ops.object.gpencil_add(location=(0, 0, 0), type='EMPTY')
        # rename grease pencil
        bpy.context.scene.objects[-1].name = gpencil_obj_name

    # Get grease pencil object
    gpencil = bpy.context.scene.objects[gpencil_obj_name]

    return gpencil


def get_grease_pencil_layer(gpencil: bpy.types.GreasePencil, gpencil_layer_name='GP_Layer',
                            clear_layer=False) -> bpy.types.GPencilLayer:
    """
    Return the grease-pencil layer with the given name. Create one if not already present.
    :param gpencil: grease-pencil object for the layer data
    :param gpencil_layer_name: name/key of the grease pencil layer
    :param clear_layer: whether to clear all previous layer data
    """

    # Get grease pencil layer or create one if none exists
    if gpencil.data.layers and gpencil_layer_name in gpencil.data.layers:
        gpencil_layer = gpencil.data.layers[gpencil_layer_name]
    else:
        gpencil_layer = gpencil.data.layers.new(gpencil_layer_name, set_active=True)

    if clear_layer:
        gpencil_layer.clear()  # clear all previous layer data

    # bpy.ops.gpencil.paintmode_toggle()  # need to trigger otherwise there is no frame

    return gpencil_layer


# Util for default behavior merging previous two methods
def init_grease_pencil(gpencil_obj_name='GPencil', gpencil_layer_name='GP_Layer',
                       clear_layer=True) -> bpy.types.GPencilLayer:
    gpencil = get_grease_pencil(gpencil_obj_name)
    gpencil_layer = get_grease_pencil_layer(gpencil, gpencil_layer_name, clear_layer=clear_layer)
    return gpencil_layer

def draw_line(gp_frame, p0: tuple, p1: tuple):
    # Init new stroke
    gp_stroke = gp_frame.strokes.new()
    gp_stroke.display_mode = '3DSPACE'  # allows for editing
    gp_stroke.draw_cyclic = False        # closes the stroke
    gp_stroke.line_width = linethickness
    gp_stroke.material_index = 0
    
    gp_stroke.points.add(count=2)
    gp_stroke.points[0].co = p0
    gp_stroke.points[1].co = p1
    return gp_stroke

def draw_circle(gp_frame, center: tuple, radius: float, segments: int):
    # Init new stroke
    gp_stroke = gp_frame.strokes.new()
    gp_stroke.display_mode = '3DSPACE'  # allows for editing
    gp_stroke.draw_cyclic = True        # closes the stroke
    
    gp_stroke.line_width = circlethickness
    gp_stroke.material_index = 1

    # Define stroke geometry
    angle = 2*math.pi/segments  # angle in radians
    gp_stroke.points.add(count=segments)
    for i in range(segments):
        x = center[0] + radius*math.cos(angle*i)
        y = center[1] + radius*math.sin(angle*i)
        z = center[2]
        gp_stroke.points[i].co = (x, y, z) 

    return gp_stroke




def dft(x,scale):
    X = []
    N = len(x)
    for k in range(0,N):
        re = 0
        im = 0
        for n in range(0,N):
            phi = (2*np.pi*k*n)/N
            re += x[n]*np.cos(phi)
            im -= x[n]*np.sin(phi)
        re = re*scale/N
        im = im*scale/N
        freq = k
        mag = math.sqrt(re*re + im*im)
        phase = np.arctan(im/re)
        X.append([re,im,freq,mag,phase])
    return X

def fft_mag_phase(fourier_fft,freqs):
    X = []
    for i in range(len(fourier_fft)):
        im = fourier_fft[i].imag
        re = fourier_fft[i].real
        mag = abs(fourier_fft[i])
        phase = np.arctan(im/re)
        freq = freqs[i]
        X.append([re,im,freq,mag,phase])
    return X


def epicycles(x,y,fourier,phase_offset=0):
    for i in range(len(fourier)):
        pre_x = x
        pre_y = y
        re,im,freq,mag,phase = fourier[i]
        N = freq
        radius = mag
        if(radius > 1): #dont draw if its too small
            draw_circle(gp_frame , (pre_x,pre_y,0) , radius , 20)
        x += radius  * np.cos(N * time + phase+phase_offset + np.pi/2)
        y += radius  * np.sin(N * time + phase+phase_offset + np.pi/2)
        if(radius > 1): #dont draw if its too small
            draw_circle(gp_frame , (x,y,0), radius*0.1 , 10)
        draw_line(gp_frame,(pre_x,pre_y,0),(x,y,0))
    return x,y




polisheddata = []
cleaned_data = [i.rsplit(',') for i in rawdata]
for i in cleaned_data:
    try:
        polisheddata.append([float(i[0]),float(i[1])])
    except ValueError:
        continue
    
offsetx,offsety = 100,100
scalex,scaley = 0.1,0.1
xvals,yvals = [i[0]*scalex+offsetx for i in polisheddata],[i[1]*(-1)*scaley+offsety for i in polisheddata]    

angle = 0

gp_layer = init_grease_pencil()
wavex = []
wavey = []
m = 3
y=0 
draw_hold = (400,0,0)
trace_length = 400

#xvals = [10*np.cos(2*np.pi*i/100) for i in range(0,100)]
#yvals = [10*np.sin(2*np.pi*i/100) for i in range(0,100)]

linethickness = 3
circlethickness = 2

freqsx = np.fft.fftfreq(len(xvals))
freqsy = np.fft.fftfreq(len(yvals))

fourierY = fft_mag_phase(np.fft.fft(yvals)/len(yvals),freqsy)
fourierX = fft_mag_phase(np.fft.fft(xvals)/len(xvals),freqsx)
prev_wave = draw_hold

wave_shift = 0

time_scale = 100*2*np.pi/len(fourierY)

for frame in range(1000):
    x,y=0,0
    gp_frame = gp_layer.frames.new(frame)
    
    time = frame*time_scale
    
    #epicycle construction 
    vy=epicycles(x,y,fourierY )
    hx=epicycles(0,0,fourierX,np.pi/2)

    wavey.insert(0,vy[1])
    wavex.insert(0,hx[0])
    draw_line(gp_frame , (vy[0],vy[1],0) , (wavex[0],wavey[0],0))
    draw_line(gp_frame , (hx[0],hx[1],0) , (wavex[0],wavey[0],0))
    draw_circle(gp_frame , (wavex[0],wavey[0],0) , 1,10)
    for i in range(0,len(wavex)):
        if(not prev_wave == draw_hold or prev_wave == (0,0,0)):
            draw_line(gp_frame , prev_wave , (wavex[i] + wave_shift , wavey[i] , 0))
        prev_wave = (wavex[i] + wave_shift  , wavey[i],0)
    if(len(wavex) > trace_length):
        wavex.pop()
        wavey.pop()
    
    
    
    
    '''    
    hx,hy=200,100
    for i in range(1,m):
        hpre_x = hx
        hpre_y = hy
        hre,him,hfreq,hmag,hphase = fourierX[i]
        hN = hfreq
        hradius = hmag
        draw_circle(gp_frame , (hpre_x,hpre_y,0) , hradius , 20)
        hx += hradius  * np.cos(hN * time + hphase + np.pi/2)
        hy += hradius  * np.sin(hN * time + hphase + np.pi/2)
        draw_circle(gp_frame , (hx,hy,0), hradius*0.1 , 10)
        draw_line(gp_frame,(hpre_x,hpre_y,0),(hx,hy,0))
    wavex.insert(0,hx)
    for i in range(len(wavey)):
        draw_line(gp_frame , prev_wave , (i + wave_shift,wavey[i],0))

        prev_wave = (i + wave_shift , wavey[i] ,0)
    if(len(wavey) > trace_length):
        wavey.pop()
    
    
   '''  
   
    
    
    #wavelet drawing
    
'''
    
    
'''