1##
2##  pullDemo.py: sample wrappers for pull API
3##  Copyright (C) 2010, 2009, Gordon Kindlmann
4##
5##  Permission is hereby granted, free of charge, to any person obtaining
6##  a copy of this software and associated documentation files (the
7##  "Software"), to deal in the Software without restriction, including
8##  without limitation the rights to use, copy, modify, merge, publish,
9##  distribute, sublicense, and/or sell copies of the Software, and to
10##  permit persons to whom the Software is furnished to do so, subject to
11##  the following conditions:
12##
13##  The above copyright notice and this permission notice shall be
14##  included in all copies or substantial portions of the Software.
15##
16##  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
17##  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
18##  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
19##  NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
20##  LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
21##  OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
22##  WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
23##
24
25'''
26This is really hastily written, and doesn't conform in the least
27to good conventions of Python coding. The worst part is using exit()
28instead of raising exceptions when tings go wrong.
29
30The "API" here (this python interface to the pull library) is not well
31thought-out, and won't be helpful for integrating particles in other
32applications.  Better designed interfaces should come in the future,
33as more and more of Teem develops idiomatic python interfaces.
34
35The code here is mainly to serve as a demo of using the pull library
36in Teem (whose API should be quite stable), and to simplify
37documenting using pull to produce typical results.
38'''
39
40import teem
41import ctypes
42import sys
43
44##
45## volLoad: loads in volumes
46##
47def volLoad(vlist, cachePath, kssBlurStr, verbose):
48    kSSblur = teem.nrrdKernelSpecNew()
49    if (teem.nrrdKernelSpecParse(kSSblur, kssBlurStr)):
50        estr = teem.biffGetDone('nrrd')
51        print "problem with kernel:\n%s" % estr
52        sys.exit(1)
53    vol = (ctypes.POINTER(teem.meetPullVol) * len(vlist))()
54    E = 0
55    for i in range(len(vlist)):
56        vol[i] = teem.meetPullVolNew()
57        if (not E):
58            E += teem.meetPullVolParse(vol[i], vlist[i])
59    if (not E):
60        E += teem.meetPullVolLoadMulti(vol, len(vol),
61                                       cachePath, kSSblur,
62                                       teem.nrrdBoundaryWrap, 0.0,
63                                       verbose)
64    if (E):
65        estr = teem.biffGetDone('meet')
66        print "problem with volumes:\n%s" % estr
67        sys.exit(1)
68    teem.nrrdKernelSpecNix(kSSblur)
69    return vol
70
71def volFree(vol):
72    for i in range(len(vol)):
73        vol[i] = teem.meetPullVolNix(vol[i])
74
75##
76## especParse: returns a length-4 list of the arguments to pass
77## to pullInterEnergySet. Arguments are:
78## type: from pullIterType* enum
79## r: string definition for energy function along space (especR)
80## s: string definition for energy function along scale (especS)
81## w: string definition for windowing energy function (especW)
82##
83def energyParse(**kwargs):
84    type = kwargs['type']
85    if (teem.pullInterTypeJustR == type
86        or teem.pullInterTypeUnivariate == type):
87        espec = [teem.pullEnergySpecNew(), None, None]
88    elif (teem.pullInterTypeSeparable == type):
89        espec = [teem.pullEnergySpecNew(),
90                 teem.pullEnergySpecNew(), None]
91    elif (teem.pullInterTypeAdditive == type):
92        espec = [teem.pullEnergySpecNew(),
93                 teem.pullEnergySpecNew(),
94                 teem.pullEnergySpecNew()]
95    lett = ['r', 's', 'w']
96    E = 0
97    for i in range(3):
98        if (not E and espec[i]):
99            E += teem.pullEnergySpecParse(espec[i], kwargs[lett[i]])
100    if (E):
101        estr = teem.biffGetDone('pull')
102        print "problem with infos:\n%s" % estr
103        sys.exit(1)
104    ret = [type, espec[0], espec[1], espec[2]]
105    return ret
106
107def energyFree(espec):
108    espec[1] = teem.pullEnergySpecNix(espec[1])
109    espec[2] = teem.pullEnergySpecNix(espec[2])
110    espec[3] = teem.pullEnergySpecNix(espec[3])
111    return
112
113##
114## initSet: calls (and returns value from) the appropriate
115## pullInitMethod* function. Possibilities are:
116## [teem.pullInitMethodRandom, num]
117## [teem.pullInitMethodPointPerVoxel,
118##     ppv, zslcmin, zslcmax, alongscalenum, jitter]
119## [teem.pullInitMethodGivenPos, npos]
120##
121def initSet(pctx, info):
122    if (teem.pullInitMethodRandom == info[0]):
123        return teem.pullInitRandomSet(pctx, info[1])
124    elif (teem.pullInitMethodPointPerVoxel == info[0]):
125        return teem.pullInitPointPerVoxelSet(pctx, info[1],
126                                             info[2], info[3],
127                                             info[4], info[5])
128    elif (teem.pullInitMethodGivenPos == info[0]):
129        return teem.pullInitGivenPosSet(pctx, info[1])
130
131
132## place for storing pullContext->sysParm.gamma
133gamma = 0
134
135## way of getting arguments that allows for checking whether
136## there were any bogus options
137def a(args, a, default=None):
138    if default is None:
139        ret = args.get(a)
140    else:
141        ret = args.get(a, default)
142    if a in args:
143        del args[a]
144    return ret
145
146##
147## All the set-up and commands that are required to do a pullContext run
148##
149def run(nposOut, **args):
150    global gamma
151    if (not ('vol' in args and 'info' in args and 'efs' in args)):
152        print "run: didn't get vol, info, and efs args"
153        sys.exit(1)
154
155    # learn args, with defaults
156    vol = a(args, 'vol')
157    infoStr = a(args, 'info')
158    efs = a(args, 'efs')
159    init = a(args, 'init', [teem.pullInitMethodRandom, 100])
160    energyDict = a(args, 'energy', {'type':teem.pullInterTypeJustR,
161                                    'r':'qwell:0.68'})
162    verbose = a(args, 'verbose', 1)
163    rngSeed = a(args, 'rngSeed', 42)
164    nobin = a(args, 'nobin', False)
165    noadd = a(args, 'noadd', False)
166    ac3c = a(args, 'ac3c', False)
167    usa = a(args, 'usa', False)
168    nave = a(args, 'nave', True)
169    cbst = a(args, 'cbst', True)
170    ratb = a(args, 'ratb', True)
171    lti = a(args, 'lti', False)
172    npcwza = a(args, 'npcwza', False)
173    ubfgl = a(args, 'ubfgl', False)
174    iterMin = a(args, 'iterMin', 0)
175    iterMax = a(args, 'iterMax', 100)
176    snap = a(args, 'snap', 0)
177    pcp = a(args, 'pcp', 5)
178    iad = a(args, 'iad', 10)
179    radSpace = a(args, 'radSpace', 1)
180    radScale = a(args, 'radScale', 1)
181    alpha = a(args, 'alpha', 0.5)
182    beta = a(args, 'beta', 0.5)
183    gamma = a(args, 'gamma', 1)
184    stepInitial = a(args, 'step', 1)
185    ssBack = a(args, 'ssBack', 0.2)
186    ssOppor = a(args, 'ssOppor', 1.1)
187    pbm = a(args, 'pbm', 50)
188    k00Str = a(args, 'k00', 'c4h')
189    k11Str = a(args, 'k11', 'c4hd')
190    k22Str = a(args, 'k22', 'c4hdd')
191    kSSreconStr = a(args, 'kSSrecon', 'hermite')
192    eip = a(args, 'eip', 0.00005)
193    eiphl = a(args, 'eiphl', 0)
194    fnnm = a(args, 'fnnm', 0.25)
195    edmin = a(args, 'edmin', 0.00001)
196    edpcmin = a(args, 'edpcmin', 0.01)
197    maxci = a(args, 'maxci', 15)
198    if args:
199        print
200        print "sorry, got unrecognized arguments:"
201        print args
202        sys.exit(1)
203
204    # create pullContext and set up all its state.  Its kind of silly that
205    # the pullContext is created anew everytime we want to run it, but thats
206    # because currently the pullContext doesn't have the kind of internal
207    # network of flags (like the gageContext does) that allow the context
208    # to be modifyied and re-used indefinitely.  This will be fixed after
209    # the Teem v1.11 release.
210    pctx = teem.pullContextNew()
211    energyList = energyParse(**energyDict)
212    if (teem.pullVerboseSet(pctx, verbose) or
213        teem.pullRngSeedSet(pctx, rngSeed) or
214        teem.pullFlagSet(pctx, teem.pullFlagNixAtVolumeEdgeSpace, nave) or
215        teem.pullFlagSet(pctx, teem.pullFlagConstraintBeforeSeedThresh,
216                         cbst) or
217        teem.pullFlagSet(pctx, teem.pullFlagEnergyFromStrength, efs) or
218        teem.pullFlagSet(pctx, teem.pullFlagRestrictiveAddToBins, ratb) or
219        teem.pullFlagSet(pctx, teem.pullFlagNoPopCntlWithZeroAlpha, npcwza) or
220        teem.pullFlagSet(pctx, teem.pullFlagUseBetaForGammaLearn, ubfgl) or
221        teem.pullFlagSet(pctx, teem.pullFlagBinSingle, nobin) or
222        teem.pullFlagSet(pctx, teem.pullFlagNoAdd, noadd) or
223        teem.pullFlagSet(pctx,
224                         teem.pullFlagAllowCodimension3Constraints,
225                         ac3c) or
226        teem.pullInitUnequalShapesAllowSet(pctx, usa) or
227        teem.pullIterParmSet(pctx, teem.pullIterParmMin, iterMin) or
228        teem.pullIterParmSet(pctx, teem.pullIterParmMax, iterMax) or
229        teem.pullIterParmSet(pctx, teem.pullIterParmSnap, snap) or
230        teem.pullIterParmSet(pctx, teem.pullIterParmConstraintMax, maxci) or
231        teem.pullIterParmSet(pctx, teem.pullIterParmPopCntlPeriod, pcp) or
232        teem.pullIterParmSet(pctx, teem.pullIterParmAddDescent, iad) or
233        teem.pullIterParmSet(pctx,
234                             teem.pullIterParmEnergyIncreasePermitHalfLife,
235                             eiphl) or
236        teem.pullSysParmSet(pctx, teem.pullSysParmStepInitial, stepInitial) or
237        teem.pullSysParmSet(pctx, teem.pullSysParmRadiusSpace, radSpace) or
238        teem.pullSysParmSet(pctx, teem.pullSysParmRadiusScale, radScale) or
239        teem.pullSysParmSet(pctx, teem.pullSysParmAlpha, alpha) or
240        teem.pullSysParmSet(pctx, teem.pullSysParmBeta, beta) or
241        teem.pullSysParmSet(pctx, teem.pullSysParmGamma, gamma) or
242        teem.pullSysParmSet(pctx, teem.pullSysParmEnergyIncreasePermit, eip) or
243        teem.pullSysParmSet(pctx, teem.pullSysParmFracNeighNixedMax, fnnm) or
244        teem.pullSysParmSet(pctx, teem.pullSysParmEnergyDecreaseMin, edmin) or
245        teem.pullSysParmSet(pctx, teem.pullSysParmEnergyDecreasePopCntlMin,
246                            edpcmin) or
247        teem.pullSysParmSet(pctx, teem.pullSysParmBackStepScale, ssBack) or
248        teem.pullSysParmSet(pctx, teem.pullSysParmOpporStepScale, ssOppor) or
249        teem.pullProgressBinModSet(pctx, pbm) or
250        teem.pullInterEnergySet(pctx, energyList[0], energyList[1],
251                                energyList[2], energyList[3])):
252        estr = teem.biffGetDone('pull')
253        print "problem with set-up:\n%s" % estr
254        sys.exit(1)
255
256    # The meet library offers a slightly higher-level abstraction of
257    # the pullInfo; this code parses a bunch of info specifications
258    info = (ctypes.POINTER(teem.meetPullInfo) * len(infoStr))()
259    for i in range(len(infoStr)):
260        info[i] = teem.meetPullInfoNew()
261        if (teem.meetPullInfoParse(info[i], infoStr[i])):
262            estr = teem.biffGetDone('meet')
263            print "problem with infos:\n%s" % estr
264            sys.exit(1)
265
266    # Create all the kernels needed for reconstruction
267    k00 = teem.nrrdKernelSpecNew()
268    k11 = teem.nrrdKernelSpecNew()
269    k22 = teem.nrrdKernelSpecNew()
270    kSSrecon = teem.nrrdKernelSpecNew()
271    if (teem.nrrdKernelSpecParse(k00, k00Str) or
272        teem.nrrdKernelSpecParse(k11, k11Str) or
273        teem.nrrdKernelSpecParse(k22, k22Str) or
274        teem.nrrdKernelSpecParse(kSSrecon, kSSreconStr)):
275        estr = teem.biffGetDone('nrrd')
276        print "problem with kernels:\n%s" % estr
277        sys.exit(1)
278
279    # We do assume that the volumes came in loaded (as opposed to reloading
280    # them here for every one); here we add all the volumes and infos into
281    # the pullContext
282    if (teem.meetPullVolAddMulti(pctx, vol, len(vol),
283                                 k00, k11, k22, kSSrecon) or
284        teem.meetPullInfoAddMulti(pctx, info, len(info))):
285            estr = teem.biffGetDone('meet')
286            print "problem with infos:\n%s" % estr
287            sys.exit(1)
288
289    # Final setup and the actual "pullRun" call
290    if (initSet(pctx, init) or
291        teem.pullInitLiveThreshUseSet(pctx, lti) or
292        teem.pullStart(pctx) or
293        teem.pullRun(pctx) or
294        teem.pullOutputGet(nposOut, None, None, None, 0, pctx)):
295        estr = teem.biffGetDone('pull')
296        print "problem running system:\n%s" % estr
297        sys.exit(1)
298
299    # until there's a clean way to re-use a pullContext for different runs,
300    # there's no good way for the user to call pullGammaLearn (because pctx
301    # doesn't survive outside the call to "run".  So we call pullGammaLearn
302    # it whenever it makes sense to do so and save the result in the global
303    # "gamma"
304    if (teem.pullInterTypeAdditive == pctx.contents.interType and
305        teem.pullEnergyButterworthParabola.contents.name
306        == pctx.contents.energySpecS.contents.energy.contents.name):
307        if (teem.pullGammaLearn(pctx)):
308            estr = teem.biffGetDone('pull')
309            print "problem learning gamma:\n%s" % estr
310            sys.exit(1)
311        gamma = pctx.contents.sysParm.gamma
312        print "pullGammaLearn returned", gamma
313
314    teem.nrrdKernelSpecNix(k00)
315    teem.nrrdKernelSpecNix(k11)
316    teem.nrrdKernelSpecNix(k22)
317    teem.nrrdKernelSpecNix(kSSrecon)
318    energyFree(energyList)
319    teem.pullContextNix(pctx)
320    return
321