https://github.com/thunil/ofblend
Revision aef880de2c500a7f7828fe93647299e8bc0170aa authored by Nils Thuerey on 06 December 2017, 08:13:20 UTC, committed by GitHub on 06 December 2017, 08:13:20 UTC
Update README.md
Tip revision: aef880de2c500a7f7828fe93647299e8bc0170aa authored by Nils Thuerey on 06 December 2017, 08:13:20 UTC
Merge pull request #1 from qinenergy/patch-1
Merge pull request #1 from qinenergy/patch-1
Tip revision: aef880d
flof.py
#
# FlOF example scene
#
import sys;
import os
from manta import *
from ofHelpers import *
paramUsed = [];
#
# quickstart guide:
#
# first gen, e.g.: "manta .../dataGen2Drop.py px 0 save4d 1 savehi 1 ; manta dataGen2Drop.py px 1 save4d 1 savehi 1 "
#
# then gen defo "manta .../flof.py dataid0 0 dataid1 1 mode 1 "
# check visually "manta .../flof.py dataid0 0 dataid1 1 mode 2 showgui 1 "
# test full ouput "manta .../flof.py dataid0 0 dataid1 1 mode 3 showgui 1 alpha 100 "
#
# gen other defo "manta .../flof.py dataid0 1 dataid1 0 mode 1 "
# test 2way defo "manta .../flof.py dataid0 1 dataid1 0 mode 2 showgui 1 alpha 100 "
#
# final 2way interpol "manta .../flof.py mode 3 showgui 1 twoway 1 alpha 50 writepngs 1 "
# (with png output)
#
setDebugLevel(1); # set to 2 for kernel msgs
# main prefix for outputs
screenPrefixOrg = "out";
screenPrefixOrg = getParam("prefix",screenPrefixOrg, paramUsed);
# read cmd line & os params
dataPath="./"
dataPathId = int(os.getenv('OFBLEND_PATH', 0));
dataPathId = int( getParam("datapathid",dataPathId, paramUsed) );
# win/mac paths
if dataPathId==1:
# default mac
dataPath="/Users/sinithue/dataOfBlend/"
elif dataPathId==2:
dataPath="I:/dataOfBlend/"
print("Using data path %s " % dataPath);
# optionally separate path for xl slice data
dataPathXl = int( getParam("datapathxl", 0, paramUsed) );
if dataPathXl==1:
# default mac
dataPathXl="/Users/sinithue/dataOfBlend/"
print("Using data path xl %s " % dataPathXl);
elif dataPathXl==2:
dataPathXl="I:/dataOfBlend/"
print("Using data path xl %s " % dataPathXl);
else:
dataPathXl=dataPath
# read params
# 1: source data id (can be interpreted differently)
# 2: target data id
# 3: output mode 1=run OF, 2=show defo, 3=load and gen hires, 4=hires 2way
dataId0 = int( getParam("dataid0", 0, paramUsed) );
dataId1 = int( getParam("dataid1", 1, paramUsed) );
mainMode = int( getParam("mode", 1, paramUsed) );
# main settings container for optical flow & data processing
ofsd = { }
ofsd["postProcMode"] = 0
ofsd["isosurfaceOff"] = 0
ofsd["isosurfaceMult"] = 1.
ofsd["fourdFactor"] = 0
ofsd["yFactor"] = 1.
ofsd["autoBorder"] = 0
ofsd["autoBrdRes"] = 0
ofsd["lsFactor"] = 1.
ofsd["maxDist"] = 0
ofsd["cgAccuracy"] = 0
ofsd["resBackup"] = 0
ofsd["repeatStartDist"] = 0
ofsd["postVelBlur"] = 0
ofsd["maxLSPDist"] = 0
ofsd["lspBlur"] = 0
# when loading hi-res SDF data, rescale values for new grid size
rescaleSdfSlices = False
sdfIsoOff = 0.
# output settings:
# writePngs png screen shot
# writeUni write volumes
# projPpm debugging, projections as ppm
# triMesh generate triangle triMesh
# meshPrefix filenames for blender
# sdfSmooth post processing for output, smooth sdf values
# output settings, note sync with checkerTest!
writePngs = 0
triMesh = 0
projPpm = 0
writeUni = 0
sdfSmooth = 0
sdfSmoothBord = 0
meshPrefix = "fluidsurface_final";
# remember whether we're dealing with SDFs or smoke data (from post-proc mode, which can be reset!)
haveLevelsetData = False;
# slight temporal filtering , good for liquids
doTimeSdfMerge = False
timeSdfIsFirst = True
setup = 1;
resHiresOut = 100;
# optional second interpolation dimension
doThirdLoad = dataId2 = 0;
doThirdLoad = int( getParam("thirdload", 0, paramUsed) );
# current data point positions (indexed by dataidX)
currp = [ vec3(0) , vec3(0) , vec3(0) , vec3(0) ];
if setup==1: # new two drop example
ofsd["fourdFactor"] = 1.5
resLoad = 40
resLoadXl = 2*resLoad
resDefo = 32 # could be anything...
ofsd["postProcMode"]= 1
haveLevelsetData = True
rescaleSdfSlices = True
sdfIsoOff = -0.5
sdfSmooth = 2
doTimeSdfMerge = True
dataPrefix = "out" # main data prefix
dataPrefixXl = "outxl" # hires slice data prefix
fileLoadPrefix = "%sdefo01_%03d_%03d_%03d" % ( dataPath, dataId0,dataId1,resDefo);
fileLoadInv = "%sdefo01_%03d_%03d_%03d" % ( dataPath, dataId1,dataId0,resDefo);
fileOutPrefix = fileLoadPrefix;
dataname_i0 = "%s/%s_r%03d_x%03d.uni" % (dataPathXl,dataPrefix, resLoad, dataId0 )
dataname_i1 = "%s/%s_r%03d_x%03d.uni" % (dataPathXl,dataPrefix, resLoad, dataId1 )
fileSliceLoadPrefix0= "%s/%s_r%03d_x%03d_%s.uni" % (dataPathXl, dataPrefixXl, resLoadXl, dataId0, "%04d" )
fileSliceLoadPrefix1= "%s/%s_r%03d_x%03d_%s.uni" % (dataPathXl, dataPrefixXl, resLoadXl, dataId1, "%04d" )
fileSliceIdxStart = 0;
fileSliceIdxEnd = resLoad*4.5;
else:
print("Error unkown setup %d " % (setup) ); exit(1);
# set other params more or less automatically
# main run mode
doRunOf = False;
showGui = 0;
doLoad = 1;
doSave = 0;
doCacheInputs = 0;
# main resolution , set later on depending on mode
res = -1;
# sliced loading: 0=off, 1=load single, 2=load uni sequence
doSlicedLoad = 0;
doOptLoadadv = 1; # use optimized slice load version?
#doOptLoadadv = 0; # NT_DEBUG off
# different modes for 2way interpol (0=off, 1=on, 2=naive blend)
doTwoWayBlend = 0;
# main blending param
blendAlpha = 1.0;
# faster load for debugging
debugSkipLoad = 9999; # = off
# for two dim interpol (optional, not for all setups)
thirdAlpha = 0.;
# interval for output
outputInterval = 1;
outputInterval = float( getParam("interval",outputInterval, paramUsed) );
if(outputInterval!=1.0): print("Note - changed outputInterval to %f " % outputInterval);
# get "partial load" fraction
partialLoad = -1.; # off
partialLoad = float( getParam("partialload",partialLoad, paramUsed) );
overrideToff = 0.; # accumulate offset
if mainMode==1:
# run OF and save OF data
res = resDefo;
doLoad = False;
doRunOf = True;
doSave = True;
# optional - enable caching
elif mainMode==2:
res = resLoad;
# allow main resolution override
res = int( getParam("resload",res, paramUsed) );
# apply defo from OF run to 4d data for visual check
doLoad = 1;
doSlicedLoad = 0;
doCacheInputs = 1;
doRunOf = False;
elif mainMode==3:
# load sliced hi res data and output deformed density files
doLoad = 0;
doSlicedLoad = 2;
ofsd["postProcMode"] = 0; # no post-processing! leave original data
if partialLoad>0.:
ofsd["postProcMode"] = -1; # also no rep. start frame
doRunOf = False;
#writeUni = 1; # uncomment to always write
# allow resHiresOut override
resHiresOut = int( getParam("hiresout",resHiresOut, paramUsed) );
res = resHiresOut;
else:
print("Unknown run mode %d " % mainMode); exit(1);
if mainMode>1:
blendAlpha = int( getParam("alpha", int(100.*blendAlpha), paramUsed) ) * 0.01;
if doThirdLoad:
thirdAlpha = int( getParam("thirdalpha", int(100.*thirdAlpha), paramUsed) ) * 0.01;
print("Blend alpha debug: %f %f %d "%(blendAlpha,thirdAlpha,doThirdLoad));
# calculate blend weights for 2dim 3point interpol using barycentric coords
def recalcUvw():
global uvw, blendAlpha, thirdAlpha;
# problem specific geometry
if setup==4:
r1 = Vec3(dataId0, depth0, 0);
r2 = Vec3(dataId1, depth1, 0);
r3 = Vec3(dataId2, depth2, 0);
elif setup==5:
r1 = currp[0]; r2 = currp[1]; r3 = currp[2]
else:
print("Missing coords"); exit(1);
pos = r1 + blendAlpha* (r2-r1) + thirdAlpha * (r3-r1);
# uvw = Vec3(triu,triv,triw); # "return" value
uvw = world2uvw( pos, r1,r2,r3 );
print("UVW for setup %f %f %f from alphas %f,%f" %(uvw.x,uvw.y,uvw.z, blendAlpha, thirdAlpha) );
#exit(1);
if doThirdLoad:
recalcUvw()
# copy to settings struct
ofsd["resBackup"] = res;
# override defaults
writePngs = int( getParam("writepngs",writePngs, paramUsed) );
writeUni = int( getParam("writeuni" ,writeUni , paramUsed) );
doTwoWayBlend = int( getParam("twoway" ,0 , paramUsed) );
doTWUnion = int( getParam("twunion" ,0 , paramUsed) );
if doTWUnion and not doTwoWayBlend:
doTwoWayBlend = 1
# output meshes? or ppms? (for debugging)
projPpm = int( getParam("projppm" ,projPpm , paramUsed) );
triMesh = int( getParam("trimesh" ,triMesh , paramUsed) );
debugSkipLoad = int( getParam("debugskip",debugSkipLoad, paramUsed) );
if(debugSkipLoad<9999): print("\nWARNING debug skip enabled %d!\n"%debugSkipLoad); # for quick tests
# note - defo "could" have different time dim (normally not)
ofsd["fourdFactorDefo"] = ofsd["fourdFactor"];
# for data load later on
ofsd["fourdFactorLoad"] = ofsd["fourdFactor"];
# data file names for caching
cacheprefix = "prep";
suffixThirdAlpha = ""
if doThirdLoad:
suffixThirdAlpha = ("b%03d" % (thirdAlpha*100.));
# adjust output path
screenPrefix = "%s/%s_f%01dt%01d_a%03d%s" % (dataPath,screenPrefixOrg,dataId0,dataId1, int(blendAlpha*100.), suffixThirdAlpha );
meshPrefix = "%s/%s" % (dataPath, meshPrefix);
# ========================
# main OF params
multiStep = 3;
ofsd["cgAccuracy"] = 0.01;
dim = 3; # not really, four in practice...
# projection params
ofsd["maxLSPDist"] = 10;
ofsd["lspBlur"] = 2.;
# misc parameters
ofsd["autoBorder"] = 0.1;
ofsd["autoBrdRes"] = int(res * ofsd["autoBorder"]);
ofsd["postVelBlur"] = 4;
sdfSmoothBord = ofsd["autoBrdRes"]; # just to pass value along...
ofsd["repeatStartDist"] = float(res) * ofsd["autoBorder"] * 1.; # keep float
ofsd["lsFactor"] = -0.1/20.0;
ofsd["maxDist"] = 40;
orgErr = 0.;
# solver setup
(gs,gs4d) = initRes( res, ofsd["yFactor"], ofsd["fourdFactorDefo"], 1.)
alloc4thDim = gs4d.t
if partialLoad>0.:
alloc4thDim = gs4d.t * partialLoad / 1.;
print("Partial sliced load active, alloc: %d " % (alloc4thDim) );
alloc4thDim = int(alloc4thDim)
solv = Solver(name='main', gridSize = gs, dim=dim, fourthDim=alloc4thDim );
solv.timestep = 1.0;
print("Main dimension "+vec4intString( Vec4(gs4d.x,gs4d.y,gs4d.z,alloc4thDim) ) +" , org length " +str(gs4d.t) );
timing = Timings();
# setup relics ...
flags = solv.create(FlagGrid)
flags.initDomain()
flags.fillGrid()
# allocate grid data
i0 = solv.create(Grid4Real)
if not doSlicedLoad:
i1 = solv.create(Grid4Real)
i0adv= solv.create(Grid4Real)
if doTwoWayBlend>0 or doThirdLoad>0:
i1adv = solv.create(Grid4Real)
velInv = solv.create(Grid4Vec4)
if doThirdLoad>0:
i2 = solv.create(Grid4Real)
thirdVel = solv.create(Grid4Vec4);
else:
if doTwoWayBlend>0 or doThirdLoad>0:
i1 = solv.create(Grid4Real)
i2 = solv.create(Grid4Real)
# preview intermediate step
i0adv_org=""
doIadvOrg = False;
if(doIadvOrg):
i0adv_org = solv.create(Grid4Real)
if not doSlicedLoad:
vel = solv.create(Grid4Vec4)
if doRunOf:
velTmp = solv.create(Grid4Vec4)
if doCacheInputs==1 and ofsd["postProcMode"]==3:
print("Turn off caching for smoke post-proc!"); exit(1);
# loader
(gsload,gs4dload) = initRes( resLoad, ofsd["yFactor"], ofsd["fourdFactorDefo"], 1.);
sload = Solver(name='loader', gridSize = gsload, dim=dim, fourthDim=resLoad*ofsd["fourdFactorLoad"]);
print("Data load dimension "+vec4intString(gs4dload)+" " );
# loading memory only required if preprocessede files are not available
loadMemRequired = True;
if doSlicedLoad==1:
prepname = "%s_%s%04d.uni" %(dataname_i0, cacheprefix, res);
if doCacheInputs and (os.path.isfile( prepname )):
loadMemRequired = False;
if loadMemRequired:
iload = sload.create(Grid4Real)
# effective time range should not be finer than org data (optionally full time range)
timerangeStart = int(ofsd["autoBrdRes"]*1.0) + ofsd["repeatStartDist"] + 1
timerangeStop = int( (res*ofsd["fourdFactor"] - ofsd["autoBrdRes"]-1)*1. );
timerangeStart = int( getParam("start",timerangeStart, paramUsed) );
timerangeStop = int( getParam("stop" ,timerangeStop , paramUsed) );
# optionally, allow as 0..1 factor
startfac = float( getParam("startfac", -1., paramUsed) );
stopfac = float( getParam("stopfac" , -1., paramUsed) );
if startfac>=0.:
timerangeStart = int(res*ofsd["fourdFactor"] * startfac);
if stopfac>=0.:
timerangeStop = int(res*ofsd["fourdFactor"] * stopfac);
# init offset
if timerangeStop>0. and partialLoad>0.:
alloc4thDim = gs4d.t * partialLoad / 1.;
# init first offset , use 1/2 of partial load window
overrideToff = -( timerangeStart - int( alloc4thDim * (1./2.) ) );
showGui = int( getParam("showgui",showGui, paramUsed) );
guiHideGrids = int( getParam("hidegrids", 1 , paramUsed) ); # on by default...
# ensure all cmd line params were used
checkUnusedParam(paramUsed);
# init UI
gui = 0;
if showGui and (GUI):
gui = Gui()
gui.show()
gui.setCamPos (0,0,-1.2);
if guiHideGrids:
gui.toggleHideGrids();
#gui.pause()
else:
writePngs = 0; # make sure it's off...
# output data
phiout = solv.create(LevelsetGrid)
meshout = solv.create(Mesh)
if(doIadvOrg):
phiout_org = solv.create(LevelsetGrid)
# time merge test for SDFs
if doTimeSdfMerge:
dstPrev1 = solv.create(LevelsetGrid)
dstPrev2 = solv.create(LevelsetGrid)
# allocate debugging grids
if not doSlicedLoad:
phiout1 = solv.create(LevelsetGrid)
phiout0 = solv.create(LevelsetGrid)
mesh1 = solv.create(Mesh)
mesh0 = solv.create(Mesh)
velout = solv.create(Vec3Grid)
veltout = solv.create(RealGrid)
velTmp_out = solv.create(Vec3Grid)
velTmp_tout = solv.create(RealGrid)
else:
# velocity debug
velout = solv.create(VecGrid)
veltout = solv.create(RealGrid)
if doTwoWayBlend==1:
phiout1 = solv.create(LevelsetGrid)
if doThirdLoad:
phiout1 = solv.create(LevelsetGrid)
phiout2 = solv.create(LevelsetGrid)
if doTWUnion==1:
phioutTmp = solv.create(LevelsetGrid); # tmp grid for union interpolation
if doTWUnion==1 and (doThirdLoad):
dstTmp = solv.create(LevelsetGrid); # another tmp grid needed...
# not always needed, but allocate solver for loading defo data
(gsDefoLoad,gs4dDefoLoad) = initRes( resDefo, ofsd["yFactor"], ofsd["fourdFactorDefo"], 1.)
sDefoload = Solver(name='defo_loader', gridSize = gsDefoLoad, dim=dim, fourthDim=int(gs4dDefoLoad.t) )
print("Defo (vel) load dimension "+vec4intString(gs4dDefoLoad)+" " )
# ======================================================================================================================
def my_range(start, end, step):
while start <= end:
yield start
start += step
def initOthers(dst0, dst1, dstvel, dst4d, dstvelTmp, dst4dTmp, t):
getSliceFrom4d(src=i0, dst=dst0, srct=t);
getSliceFrom4d(src=i1, dst=dst1, srct=t);
getSliceFrom4dVec(src=vel, dst=dstvel, dstt=dst4d, srct=t);
#getSliceFrom4dVec(src=velTmp, dst=dstvelTmp, dstt=dst4dTmp, srct=t);
# most parameters hidden in ofsd dictionary
def postProcessPhiData(iTarget, isSlice=False , loadOffset=vec4(0) , ofsd={} ):
printMinMax = True;
if(printMinMax): print("Min max before post proc: %f to %f" % (iTarget.getMin(), iTarget.getMax()) );
if ofsd["postProcMode"]==-1:
print("Post processing off, no rep. start. fr.");
return; # special for streaming load
# start frame , note - make sure it hits full data!
repeatFrame4d( iTarget, loadOffset.t + 1./float(ofsd["resBackup"]), ofsd["repeatStartDist"], bnd=0 );
if ofsd["postProcMode"]==0:
print("Post processing off");
elif ofsd["postProcMode"]==1:
# extrapolate levelset
print("Post processing as level set, offset %f, mult %f " %(ofsd["isosurfaceOff"],ofsd["isosurfaceMult"]) );
if ofsd["isosurfaceOff"] != 0.:
iTarget.addConst(ofsd["isosurfaceOff"]);
if ofsd["isosurfaceMult"] != 1.:
iTarget.multConst(ofsd["isosurfaceMult"]);
iTarget.setBound(0.1, ofsd["autoBrdRes"] ); # reset border of LS to outside
# extrapolate
extrap4dLsSimple(phi=iTarget, distance=ofsd["maxDist"])
extrap4dLsSimple(phi=iTarget, distance=ofsd["maxDist"], inside=True)
iTarget.multConst( ofsd["lsFactor"] );
elif ofsd["postProcMode"]==3:
print("Post processing with smokeWLT special!");
# convert from "LS" to smoke
iTarget.addConst (-0.125)
iTarget.multConst (-1.0 * ofsd["lsFactor"]); # undo mult later on!
iTarget.setBound (0. , ofsd["autoBrdRes"]);
# note - does special set bound later on!
if(printMinMax): print("Min max after post proc: %f to %f" % (iTarget.getMin(), iTarget.getMax()) );
# load data sets & registrations
# regular 4d data load (also handles caching)
def loadPhiData(iTarget, fname ):
# load & resize
global iload, loadOffset;
prepname = "%s_%s%04d.uni" %(fname,cacheprefix,res);
if doCacheInputs and (os.path.isfile( prepname )):
iTarget.load(prepname);
return;
iload.load( fname );
# special - reset sides before rescaling!
if ofsd["postProcMode"]==3:
iload.setBound(0.125,1); # neutralize offset
interpolateGrid4d(iTarget, iload, offset=loadOffset , scale=loadScale );
postProcessPhiData( iTarget, False, loadOffset, ofsd )
if doCacheInputs:
iTarget.save(prepname);
# sliced load from 3d sequence (for hi-res data)
def loadPhiDataSliced( phiDst, phiSliceName ):
# loop over given range and load
loadPlaceGrid4d(fname=phiSliceName , phi=phiDst, offset=loadOffset , scale=loadScale , \
fileIdxStart=fileSliceIdxStart, fileIdxEnd=fileSliceIdxEnd , debugSkipLoad=debugSkipLoad , spread=1. , \
overrideSize=gs4d , overrideTimeOff=overrideToff , overrideGoodRegion=0 , loadTimeScale=1. , rescaleSdfValues=rescaleSdfSlices, sdfIsoOff=sdfIsoOff, repeatStartFrame=ofsd["autoBorder"] );
postProcessPhiData( phiDst, True, loadOffset, ofsd )
# simple quadratic ease in/out , t = time, b = startvalue, c = change in value, d = duration
def easeInOut (t, b, c, d):
t = t / (d*0.5);
if t < 1.:
return( c/2.*t*t + b );
t = t-1;
return( -c/2. * (t*(t-2.) - 1.) + b);
# blend weight
def blendW (t, l):
w0 = ( 1.0+l - 2.0*t) / (1.0+l);
w1 = (-1.0+l + 2.0*t) / (1.0+l);
if(w0<0.): w0=0.;
if(w1<0.): w1=0.;
if(w0>1.): w0=1.;
if(w1>1.): w1=1.;
wu = 1. - w0 - w1;
print("Blend for %f: w0 %f, wu %f, w1 %f " % (t,w0,wu,w1) );
return (w0,wu,w1);
# ======================================================================================================================
# apply deformation to whole 4d volume (using global params)
def deform4d():
i0adv.copyFrom(i0);
if blendAlpha<(0.+1e-03) and not doThirdLoad:
print("Showing i0!\n"); # debug info, show input 1 , target
elif blendAlpha<(1.+1e-03) or doThirdLoad: # for third-load, always apply (dont skip to i1)
print("Applying defo with alpha=%f !\n" % blendAlpha);
if blendAlpha>(1.0-1e-03) and ofsd["postProcMode"]==0:
massBefore = debugGridAvg4d(i0adv,ofsd["autoBrdRes"]);
advect4d(vel=vel, grid=i0adv, dtFac=blendAlpha);
if doTwoWayBlend==1:
# reverse advect
i1adv.copyFrom(i1);
advect4d(vel=velInv, grid=i1adv, dtFac= (1.-blendAlpha) );
i0adv.multConst(1.-blendAlpha);
i0adv.addScaled(i1adv, blendAlpha);
elif doTwoWayBlend==2: # naive blend, for comparison only
print("\nShowing naive blend result!\n" );
i0adv.multConst(1.-blendAlpha);
i0adv.addScaled(i1, blendAlpha);
# only for full advect, recalc diff
if blendAlpha>(1.0-1e-03) and (orgErr>1e-03):
if ofsd["postProcMode"]==0:
massAfter = debugGridAvg4d(i0adv,ofsd["autoBrdRes"]);
print("Avg mass before/after adv %f,%f; factor %f" % (massBefore,massAfter,(massBefore/massAfter)) );
avgErr = calcLsDiff4d(i0=i0adv, i1=i1 );
print("Final ls diff = %f , org fac %f\n" % (avgErr, (avgErr/orgErr) ) );
else:
# debug, show input 1 , target
i0adv.copyFrom(i1); print("Showing i1!\n");
# second defo for two-dim data (note - copied from first defo)
if doThirdLoad:
if thirdAlpha<(0.+1e-03):
print("Not active: vel2!\n");
print("Applying defo2 with third-alpha=%f !\n" % thirdAlpha);
advect4d(vel=vel, grid=thirdVel, dtFac=blendAlpha ); # align
advect4d(vel=thirdVel, grid=i0adv, dtFac=thirdAlpha );
def dumpOut(iadv, outName, interval=1, nrOverride=-1):
meshCnt = 0;
for t in my_range( timerangeStart,timerangeStop, interval):
getSliceFrom4d(src=iadv, dst=phiout, srct=t);
phiout.createMesh(meshout);
if(doIadvOrg):
getSliceFrom4d(src=i0adv_org, dst=phiout_org, srct=t);
# show orgs , debugging!
if not triMesh:
getSliceFrom4d(src=i0 , dst=phiout0, srct=t);
getSliceFrom4d(src=i1 , dst=phiout1, srct=t);
phiout0.createMesh(mesh0); phiout1.createMesh(mesh1);
initOthers(phiout0, phiout1, velout,veltout, velTmp_out,velTmp_tout, t);
if nrOverride<0:
frameNum = int(t/interval);
if (frameNum != t):
print("NYI- only integer time interval supported here !"); exit(1);
writeResults(outName, frameNum, meshCnt, phiout,meshout,gui, writePngs , triMesh , projPpm , writeUni , sdfSmooth,sdfSmoothBord , meshPrefix );
meshCnt = meshCnt+1;
else:
print("Nr override active!");
writeResults(outName, nrOverride,nrOverride, phiout, meshout,gui, writePngs , triMesh , projPpm , writeUni , sdfSmooth,sdfSmoothBord , meshPrefix );
solv.step();
# similar to dumpOut, but more functionality (optimized loads, multiple defos)
isInitialized = False;
def loadSlice_dumpOut( dst, phiSrc, defoName0, outName, interval=1, phiSliceName="", \
defoName1="", dst1=0, phiSrc1=0, phiSliceName1="", \
defoName2="", dst2=0, phiSrc2=0, phiSliceName2="", \
doCleanup=True , nrOverride=-1 , \
defoName3="", dst3=0, phiSrc3=0, phiSliceName3="" , defoAniFac = 1.):
global overrideToff, isInitialized, timeSdfIsFirst;
meshCnt = 0;
lats1 = 7; # arbitrary IDs
lats2 = 13; # actual value doesnt matter
lats = [7, 13, 73, 11, 17 ]; # last one unused for now...
latName = [defoName0, defoName1, defoName2, defoName3, "-" ];
latSrc = [phiSrc, phiSrc1, phiSrc2, phiSrc3, 0 ];
latDst = [dst, dst1, dst2, dst3, 0 ];
latSlic = [phiSliceName, phiSliceName1, phiSliceName2, phiSliceName3, 0 ];
latNum = 0;
# offset by fraction of single frame
frameOffset = 0.5;
loadFunc = loadAdvectTimeSlice;
if doOptLoadadv: loadFunc = loadAdvectTimeSlice_OptRun; # optimized version
latNum = 1;
if doTwoWayBlend: latNum = 2;
if doThirdLoad : latNum = 3;
useDefoVols = (latNum>2)
if not isInitialized:
for i in range(0,latNum):
print("Advect init %d %s " % (i,latName[i]) );
loadAdvectTimeSlice_OptInit( lats[i], latName[i] , useDefoVols, False, partialLoad );
if useDefoVols:
loadAdvectTimeSlice_OptAdd( lats[i], latName[(i+1)%latNum] );
isInitialized = True
for t in my_range( timerangeStart,timerangeStop, interval):
print("Time %f, in range %f to %f, inter %f " % (t,timerangeStart,timerangeStop,interval) );
# activate partial loading , shift & load whenever necessary
print("Pos in slice %d , dim %d , toff %d " % ( (t+overrideToff) , alloc4thDim, overrideToff) );
keep = int(alloc4thDim-2);
keepCheck = int(alloc4thDim*0.5)+1
if partialLoad>0. and (t + overrideToff) > keepCheck:
shift = (alloc4thDim - keep); # shift amount (remainder)
overrideToff = overrideToff - shift;
print("Pos shift %d " % (shift) );
for i in range(0,latNum):
shiftForwGrid4d( phi=latSrc[i] , overrideGoodRegion=keep );
for i in range(0,latNum): loadPlaceGrid4d(fname=latSlic[i] , phi=latSrc[i], offset=loadOffset , scale=loadScale , \
fileIdxStart=fileSliceIdxStart, fileIdxEnd=fileSliceIdxEnd , debugSkipLoad=debugSkipLoad , spread=1. , \
overrideSize=gs4d , overrideTimeOff=overrideToff , overrideGoodRegion=keep, loadTimeScale=1. , rescaleSdfValues=rescaleSdfSlices ,sdfIsoOff=sdfIsoOff, repeatStartFrame=ofsd["autoBorder"] );
# interpolation go...
if not doThirdLoad:
# 2 data sets , 1d line interpolations
loadFunc ( lats[0], latName[0], latDst[0], latSrc[0], t+frameOffset*interval, blendAlpha , 1. , \
defoOffset, defoScale, defoFactor, overrideSize=gs4d, overrideTimeOff=overrideToff, zeroVel=False , thirdAlpha = thirdAlpha, bordSkip=ofsd["autoBrdRes"] , defoAniFac=defoAniFac );
alphaInv = 1. - blendAlpha;
if doTwoWayBlend==2: # recreate naive blend , alpha = 0 for second data
alphaInv = 0.;
if doTwoWayBlend==1 and doTWUnion:
loadFunc ( lats[1], latName[1], latDst[1], latSrc[1], t+frameOffset*interval, alphaInv , 1. , defoOffset, defoScale, defoFactor, gs4d, overrideToff, zeroVel=False, bordSkip=ofsd["autoBrdRes"] , defoAniFac=defoAniFac );
(w0,wu,w1) = blendW(blendAlpha, 0.1);
# blend union
phioutTmp.copyFrom(dst);
phioutTmp.join(latDst[1]);
dst.multConst( w0);
dst.addScaled(phioutTmp, wu);
dst.addScaled(latDst[1], w1);
print("Blended union with weights %f,%f,%f " % (w0,wu,w1) );
elif doTwoWayBlend==1 or doTwoWayBlend==2:
loadFunc ( lats[1], latName[1], latDst[1], latSrc[1], t+frameOffset*interval, alphaInv , 1. , defoOffset, defoScale, defoFactor, gs4d, overrideToff, zeroVel=False, bordSkip=ofsd["autoBrdRes"] , defoAniFac=defoAniFac );
# include additional blending factors
dst.multConst( (1.-blendAlpha) );
dst.addScaled(latDst[1], ( blendAlpha) );
elif doThirdLoad: # doThirdLoad , 3 data sets
if doTwoWayBlend==0: # simple interpol , 1 src data point only
loadFunc ( lats[0], latName[0], latDst[0], latSrc[0], t+frameOffset*interval, (uvw.y + uvw.z) , 1. ,defoOffset, defoScale, defoFactor, overrideSize=gs4d, overrideTimeOff=overrideToff, zeroVel=False , thirdAlpha = uvw.z , bordSkip=ofsd["autoBrdRes"] , defoAniFac=defoAniFac );
if doTwoWayBlend==1 and doTWUnion:
print("Advect stanton 2dim with uvws %f %f %f " % (uvw.x,uvw.y,uvw.z) );
loadFunc ( lats[0], latName[0], latDst[0], latSrc[0], t+frameOffset*interval, (uvw.y + uvw.z) , 1. ,defoOffset, defoScale, defoFactor, overrideSize=gs4d, overrideTimeOff=overrideToff, zeroVel=False , thirdAlpha = uvw.z , bordSkip=ofsd["autoBrdRes"] , defoAniFac=defoAniFac );
loadFunc ( lats[1], latName[1], latDst[1], latSrc[1], t+frameOffset*interval, (uvw.z + uvw.x) , 1. ,defoOffset, defoScale, defoFactor, overrideSize=gs4d, overrideTimeOff=overrideToff, zeroVel=False , thirdAlpha = uvw.x , bordSkip=ofsd["autoBrdRes"] , defoAniFac=defoAniFac );
loadFunc ( lats[2], latName[2], latDst[2], latSrc[2], t+frameOffset*interval, (uvw.x + uvw.y) , 1. ,defoOffset, defoScale, defoFactor, overrideSize=gs4d, overrideTimeOff=overrideToff, zeroVel=False , thirdAlpha = uvw.y , bordSkip=ofsd["autoBrdRes"] , defoAniFac=defoAniFac );
# calculate blendW
(w0,w1,w2,w_uv,w_vw,w_uw,w_uvw) = weights2dim(uvw)
sum =w0+w1+w2+w_uv+w_vw+w_uw+w_uvw;
print( "uvw-summ %4.3f,%4.3f,%4.3f w0=%4.3f, w1=%4.3f, w2=%4.3f, w_uv=%4.3f, w_vw=%4.3f, w_uw=%4.3f, w_uvw=%4.3f , sum %f" % (uvw.x,uvw.y,uvw.z, w0,w1,w2,w_uv,w_vw,w_uw,w_uvw, sum) );
# calculate unions w weight >0
dstTmp.setConst(0.);
if w_uv>1e-03:
phioutTmp.copyFrom (latDst[0]);
phioutTmp.join (latDst[1]);
dstTmp.addScaled(phioutTmp, w_uv );
if w_vw>1e-03:
phioutTmp.copyFrom (latDst[1]);
phioutTmp.join (latDst[2]);
dstTmp.addScaled(phioutTmp, w_vw );
if w_uw>1e-03:
phioutTmp.copyFrom (latDst[0]);
phioutTmp.join (latDst[2]);
dstTmp.addScaled(phioutTmp, w_uw );
if w_uvw>1e-03:
phioutTmp.copyFrom (latDst[0]);
phioutTmp.join (latDst[1]);
phioutTmp.join (latDst[2]);
dstTmp.addScaled(phioutTmp, w_uvw );
# last, add orgs
if w0>1e-03: dstTmp.addScaled(latDst[0], w0 );
if w1>1e-03: dstTmp.addScaled(latDst[1], w1 );
if w2>1e-03: dstTmp.addScaled(latDst[2], w2 );
dst.copyFrom(dstTmp); # finally copy over, overwrite deformed i0
elif doTwoWayBlend==1: # simple interpol
print("Advect 2dim with uvws %f %f %f " % (uvw.x,uvw.y,uvw.z) );
loadFunc ( lats[0], latName[0], latDst[0], latSrc[0], t+frameOffset*interval, (uvw.y + uvw.z) , 1. ,defoOffset, defoScale, defoFactor, overrideSize=gs4d, overrideTimeOff=overrideToff, zeroVel=False , thirdAlpha = uvw.z , bordSkip=ofsd["autoBrdRes"] , defoAniFac=defoAniFac );
loadFunc ( lats[1], latName[1], latDst[1], latSrc[1], t+frameOffset*interval, (uvw.z + uvw.x) , 1. ,defoOffset, defoScale, defoFactor, overrideSize=gs4d, overrideTimeOff=overrideToff, zeroVel=False , thirdAlpha = uvw.x , bordSkip=ofsd["autoBrdRes"] , defoAniFac=defoAniFac );
loadFunc ( lats[2], latName[2], latDst[2], latSrc[2], t+frameOffset*interval, (uvw.x + uvw.y) , 1. ,defoOffset, defoScale, defoFactor, overrideSize=gs4d, overrideTimeOff=overrideToff, zeroVel=False , thirdAlpha = uvw.y , bordSkip=ofsd["autoBrdRes"] , defoAniFac=defoAniFac );
if 0:
# NT DEBUG, only dst , leave unmodified
print("Note: No blend, only dst active!");
else:
dst.multConst( uvw.x );
dst.addScaled(latDst[1], uvw.y );
dst.addScaled(latDst[2], uvw.z );
if doTwoWayBlend==2: # not for naive blend!
print("No blend, naive on");
if doTimeSdfMerge: # "motion blur", temporal blur effect
if timeSdfIsFirst:
dstPrev2.copyFrom(dst)
dstPrev1.copyFrom(dst)
timeSdfIsFirst = False
dstPrev2.copyFrom(dstPrev1)
dstPrev1.copyFrom(dst)
dst.join(dstPrev2)
simpleBlurSpecial( dst, 1 , 0.);
print("Time merge sdf mode enabled!");
if haveLevelsetData == 1: # for levelset data reset border for LS
#print("Set bound %d "%(ofsd["autoBrdRes"]+50));
dst.setBound(0.5, ofsd["autoBrdRes"]+1);
if nrOverride<0:
frameNum = int(t/interval);
writeResults(outName, frameNum,meshCnt, dst, meshout,gui, writePngs , triMesh , projPpm , writeUni , sdfSmooth,sdfSmoothBord , meshPrefix );
meshCnt = meshCnt+1;
else:
print("Nr override active!");
writeResults(outName, nrOverride,nrOverride, dst, meshout,gui, writePngs , triMesh , projPpm , writeUni , sdfSmooth,sdfSmoothBord , meshPrefix );
#timing.display()
solv.step();
# cleanup globals
if doOptLoadadv and doCleanup:
for i in range(0,latNum):
loadAdvectTimeSlice_Finish( lats[i] );
# ======================================================================================================================
# main
# load input data
# first init offsets
(loadOffset,loadScale,defoOffset,defoScale,defoFactor) = initOffsets( gs4d,gs4dDefoLoad, ofsd["autoBorder"], ofsd["repeatStartDist"], 1.);
print("Data scaling factor "+str(loadScale)+", offset " +str(loadOffset)+" " );
# i0
if not doSlicedLoad==2:
loadPhiData(i0, dataname_i0);
else:
loadPhiDataSliced(i0, fileSliceLoadPrefix0 );
if doTwoWayBlend==1:
loadPhiDataSliced(i1, fileSliceLoadPrefix1 );
orgErr = 1.;
# i1
if not doSlicedLoad:
loadPhiData(i1, dataname_i1);
erri0i1 = calcLsDiff4d(i0=i0, i1=i1 );
orgErr = erri0i1;
print("Error between inputs %f " % (erri0i1) );
# i2
if doThirdLoad:
if not doSlicedLoad:
loadPhiData(i2, dataname_i2);
else:
loadPhiDataSliced(i2, fileSliceLoadPrefix2 );
# load previous defo data
if doLoad:
loadVelRescaled(vel, fileLoadPrefix, sDefoload, defoOffset, defoScale, defoFactor);
if doTwoWayBlend==1:
# load reverse defo
loadVelRescaled(velInv, fileLoadInv , sDefoload, defoOffset, defoScale, defoFactor);
print("Two-way (1dim) interpolation active" );
if doThirdLoad==1:
loadVelRescaled(thirdVel, fileLoadThird1 , sDefoload, defoOffset, defoScale, defoFactor);
print("Two-dim interpolation active" );
vel.setBoundNeumann(ofsd["autoBrdRes"]+0);
if doTwoWayBlend==1: velInv.setBoundNeumann(ofsd["autoBrdRes"]+0);
if doThirdLoad>0: thirdVel.setBoundNeumann(ofsd["autoBrdRes"]+0);
if loadMemRequired:
del iload;
# run of!
if doRunOf:
if ofsd["postProcMode"]>=2:
print("ERROR - smoke post-proc was active! Needs to be disabled for actual OF runs"); exit(1);
# run!
if doRunOf and 1:
opticalFlowMultiscale4d(vel=vel, i0=i0, i1=i1, wSmooth=0.001, wEnergy=0.0001 , cgAccuracy=ofsd["cgAccuracy"] , \
postVelBlur=ofsd["postVelBlur"], cfl=999, multiStep=multiStep , minGridSize=20,
doFinalProject=True, resetBndWidth=ofsd["autoBorder"] );
# main OF step done
solv.step()
if doSave:
vel.save( "%s_vel.uni" % fileOutPrefix );
# before applying defo, undo scaling for OF
if not doSlicedLoad:
i0.multConst( 1.0/ofsd["lsFactor"] );
i1.multConst( 1.0/ofsd["lsFactor"] );
if doTWUnion or doThirdLoad:
print("These modes are only supported with slices!"); exit(1);
# ======================================================================================================================
# get 3d slices and output result in some way
if ( showGui or triMesh or projPpm or writeUni or writePngs ):
doClean = True;
defoAniFac = 1.0
steps = 1;
# iterate for single step mode
nrOverride = -1;
for step in range(steps):
# apply deformation once
if step==0 and not doSlicedLoad:
deform4d()
# actual output modes
if not doSlicedLoad:
dumpOut(i0adv, screenPrefix, outputInterval, nrOverride);
else:
if doThirdLoad>0:
loadSlice_dumpOut( dst=phiout, phiSrc=i0, defoName0=("%s_vel.uni" % fileLoadPrefix), outName=screenPrefix, \
interval=outputInterval, phiSliceName=fileSliceLoadPrefix0 ,
defoName1=("%s_vel.uni" % fileLoadThird1), dst1=phiout1, phiSrc1=i1, phiSliceName1=fileSliceLoadPrefix1 ,
defoName2=("%s_vel.uni" % fileLoadThird2), dst2=phiout2, phiSrc2=i2, phiSliceName2=fileSliceLoadPrefix2 , nrOverride=nrOverride, doCleanup=doClean, defoAniFac=defoAniFac );
elif doTwoWayBlend>0:
loadSlice_dumpOut( dst=phiout, phiSrc=i0, defoName0=("%s_vel.uni" % fileLoadPrefix), outName=screenPrefix, \
interval=outputInterval, phiSliceName=fileSliceLoadPrefix0, \
defoName1=("%s_vel.uni" % fileLoadInv), dst1=phiout1, phiSrc1=i1, phiSliceName1=fileSliceLoadPrefix1 , nrOverride=nrOverride, doCleanup=doClean , defoAniFac=defoAniFac );
else:
loadSlice_dumpOut( dst=phiout, phiSrc=i0, defoName0=("%s_vel.uni" % fileLoadPrefix), outName=screenPrefix, \
interval=outputInterval, phiSliceName=fileSliceLoadPrefix0 , nrOverride=nrOverride, doCleanup=doClean , defoAniFac=defoAniFac );
else:
print("No frame output.");
# done!
solv.step()
Computing file changes ...