from rsf.proj import *
import math

ax=15
ay=25
az=35

Flow('phase',None,
     '''
     spike n1=64 n2=64 n3=64 d1=1 d2=1 d3=1 unit1= unit2= unit3= |
     rtoc |
     fft3 axis=1 pad=1 |
     fft3 axis=2 pad=1 |
     fft3 axis=3 pad=1 |
     real |
     noise rep=y seed=2010 type=n |
     rtoc |
     math output="exp(%g*I*input)"
     ''' % (2*math.pi))

def model(method,ax=15,ay=25,az=35,hurst=0.1):
    if method==0:
        equation="phase*exp(-input/8)"
    elif method==1:
        equation="phase/(1+input)^%g" % (0.75+0.5*hurst)
    else:
        equation="phase/(1+input)"

    return '''
    math output="%g*x1*x1+%g*x2*x2+%g*x3*x3" |
    math phase=$SOURCE output="%s" |
    fft3 axis=3 inv=y |
    fft3 axis=2 inv=y |
    fft3 axis=1 inv=y |
    real
    ''' % (ax*ax,ay*ay,az*az,equation)

Flow('gauss','phase',model(0))
Flow('fractal','phase',model(1))
Flow('expo', 'phase',model(2))
Flow('layered','phase',model(2,5,80,80))
Flow('isotropic','phase',model(0,15,15,15))

for case in Split('gauss expo fractal layered isotropic'):
    Result(case,
           '''
           byte |
           grey3 flat=n label1=X label2=Y label3=Z wanttitle=n
           ''')

autocorr =  '''
fft1 | 
math output="input*conj(input)" |
fft1 inv=y |
stack axis=3 |
stack axis=2 |
window n1=50 |
scale axis=1
'''

for case in Split('layered isotropic'):
    bin = case[:3]+'-bin'
    Flow(bin,case,'mask min=0')
    Result(bin,
           '''
           dd type=float |
           byte allpos=y |
           grey3 flat=n label1=X label2=Y label3=Z wanttitle=n
           ''')

    auto = case[:3]+'-auto'
    Flow(auto,case,autocorr)

    autobin = auto+'bin'
    Flow(autobin,bin,'dd type=float | add add=-0.5 | ' + autocorr)

    Result(auto,[auto,autobin],
           '''
           cat axis=2 ${SOURCES[1]} |
           graph dash=0,1 title="Autocorrelation"
           label1=Lag label2="Normalized Autocorrelation"
           ''')

End()