"""denoise test in 3d model"""
from rsf.prog import RSFROOT
from rsf.proj import *

###
# To finishing fx-emdpf denoising part, it needs to install Matlab,
# and rebuild Madagascar with "API=Matlab" option.
###

### load data
Flow("data", "image3d", "dd form=native")
Result("realdata", "data",
       """
       byte gainpanel=e clip=1.5 |
       grey3 flat=y frame1=250 frame2=47 frame3=254
       point1=0.7 point2=0.45 title= color=g
       label2=X label3=Y unit2=km unit3=km
       """)

### f-x SPF
Flow("spf2", "data",
     """
     fxspfdenoise2 lambda1=280 lambda3=55 na1=4 na2=0 verb=y 
     """)
Result("realspf2", "spf2",
       """
       byte gainpanel=e clip=1.5 |
       grey3 flat=y frame1=250 frame2=47 frame3=254
       point1=0.7 point2=0.45 title= color=g
       label2=X label3=Y unit2=km unit3=km
       """)
Flow("errspf2", "data spf2", "add ${SOURCES[1]} scale=1,-1 ")
Result("realerrspf2", "errspf2",
       """
       byte gainpanel=e clip=1.5 |
       grey3 flat=y frame1=250 frame2=47 frame3=254
       point1=0.7 point2=0.45 title= color=g
       label2=X label3=Y unit2=km unit3=km
       """)

### f-x-y SPF
# padding data
length_pad = 100
Flow("data_pad", "data",
     """
     window n3=%d | reverse which=4 | put o3=0 d3=0.03
     """ % length_pad)
Flow("datas", "data_pad data", "cat ${SOURCES[1:2]} axis=3")
# denoise
Flow("spf3", "datas",
     """
     fxyspfdenoise3 lambda1=600 lambda2=600 lambda3=90 na1=2 na2=2 verb=y |
     window f3=%d | put o3=0
     """ % length_pad)
Result("realspf3", "spf3",
       """
       byte gainpanel=e clip=1.5 |
       grey3 flat=y frame1=250 frame2=47 frame3=254
       point1=0.7 point2=0.45 title= color=g
       label2=X label3=Y unit2=km unit3=km
       """)
Flow("errspf3", "data spf3", "add ${SOURCES[1]} scale=1,-1 ")
Result("realerrspf3", "errspf3",
       """
       byte gainpanel=e clip=1.5 |
       grey3 flat=y frame1=250 frame2=47 frame3=254
       point1=0.7 point2=0.45 title= color=g
       label2=X label3=Y unit2=km unit3=km
       """)

### f-x EMD (emdpf)
matlab = WhereIs("matlab")
matROOT = "../Matfun"
matfun = "Real"
matlabpath = os.environ.get("MATLABPATH", os.path.join(RSFROOT, "lib"))

if not matlab:
    sys.stderr.write("\nCannot find Matlab.\n")
    sys.exit(1)

fxemds = []
for n in range(0, 310):
    fxdat = "fxdat-%s" % n
    fxemd = "fxemd-%s" % n
    Flow(fxdat, "data", "window n3=1 f3=%d" % n)
    Flow(fxemd, [fxdat, os.path.join(matROOT, matfun + ".m")],
         """
         MATLABPATH=%(matlabpath)s %(matlab)s -nosplash -nojvm -r
         "addpath %(matROOT)s;%(matfun)s('${SOURCES[0]}','${TARGETS[0]}');quit"
         """ % vars(),
         stdin=0,
         stdout=-1)
    fxemds.append(fxemd)

Flow("fxemd", fxemds, "cat ${SOURCES[1:%d]} axis=3 " % len(fxemds))
Flow("errfxemd", "data fxemd", "add ${SOURCES[1]} scale=1,-1 ")
Result("realfxemd", "fxemd",
       """
       put o1=0 o2=0 o3=0 d1=0.001 d2=0.03 d3=0.03 |
       byte gainpanel=e clip=1.5 |
       grey3 flat=y frame1=250 frame2=47 frame3=254
       point1=0.7 point2=0.45 title= color=g
       label2=X label3=Y unit2=km unit3=km
       """)
Result("realerrfxemd", "errfxemd",
       """
       byte gainpanel=e clip=1.5 |
       grey3 flat=y frame1=250 frame2=47 frame3=254
       point1=0.7 point2=0.45 title= color=g
       label2=X label3=Y unit2=km unit3=km
       """)

### f-x-y RNA  (consume large memory space)
# Patch  f=73*2/5 (patching for reducing computational burden)
Flow('patch','data','patch w=140,266,310 p=10,1,1')
tpres = []
tpre2ds = []

for nwt in range(0,10):
    fd     = 'fd-%d' % nwt
    shiftsa= 'shiftsa-%d' % nwt
    sh1    = 'sh1-%d' % nwt
    shifts = 'shifts-%d' % nwt
    flt    = 'flt-%d' % nwt
    pre    = 'pre-%d' % nwt
    tpre   = 'tpre-%d' % nwt
    Flow(fd,'patch',
         '''
         window n4=1 f4=%d | fft1 | transp plane=13 memsize=1000 |
         transp plane=12 memsize=1000 
         '''  % nwt )

    Flow(shifts,fd,
        '''
        cshifts2 ns1=2 ns2=2 | transp plane=34 memsize=1000
        ''' )

    Flow([flt, pre],[shifts, fd],
         '''
         clpf match=${SOURCES[1]} pred=${TARGETS[1]} rect1=10 rect2=10 niter=10
         ''')
    Flow(tpre,pre,
         '''
         transp plane=13 memsize=1000 | transp plane=23 memsize=1000 | 
         fft1 inv=y 
         ''')
    tpres.append(tpre)
Flow('rna',tpres,
     'cat ${SOURCES[1:%d]} axis=4 | patch inv=y weight=y dim=3' % len(tpres))

Result("realrna", "rna",
       """
       byte gainpanel=e clip=1.5 |
       grey3 flat=y frame1=250 frame2=47 frame3=254
       point1=0.7 point2=0.45 title= color=g
       label2=X label3=Y unit2=km unit3=km
       """)
Flow("errrna", "data rna", "add ${SOURCES[1]} scale=1,-1 ")
Result("realerrrna", "errrna",
       """
       byte gainpanel=e clip=1.5 |
       grey3 flat=y frame1=250 frame2=47 frame3=254
       point1=0.7 point2=0.45 title= color=g
       label2=X label3=Y unit2=km unit3=km
       """)






End()