← Back to branch summary

~ubuntu-branches/ubuntu/raring/voxbo/raring

« back to all changes in this revision

Viewing changes to munge/vbimagemunge.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Michael Hanke
  • Date: 2010-06-06 11:33:11 UTC
  • Revision ID: james.westby@ubuntu.com-20100606113311-v3c13imdkkd5n7ae
Tags: upstream-1.8.5~svn1172
Import upstream version 1.8.5~svn1172

Show diffs side-by-side

added added

removed removed

Lines of Context:
munge/vbimagemunge.cpp
 
1
 
 
2
// vbimagemunge.cpp
 
3
// general-purpose munging util for voxbo
 
4
// Copyright (c) 2003-2009 by The VoxBo Development Team
 
5
 
 
6
// VoxBo is free software: you can redistribute it and/or modify it
 
7
// under the terms of the GNU General Public License as published by
 
8
// the Free Software Foundation, either version 3 of the License, or
 
9
// (at your option) any later version.
 
10
// 
 
11
// VoxBo is distributed in the hope that it will be useful, but
 
12
// WITHOUT ANY WARRANTY; without even the implied warranty of
 
13
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 
14
// General Public License for more details.
 
15
// 
 
16
// You should have received a copy of the GNU General Public License
 
17
// along with VoxBo.  If not, see <http://www.gnu.org/licenses/>.
 
18
// 
 
19
// For general information on VoxBo, including the latest complete
 
20
// source code and binary distributions, manual, and associated files,
 
21
// see the VoxBo home page at: http://www.voxbo.org/
 
22
//
 
23
// original version written by Dan Kimberg
 
24
 
 
25
#include <fstream>
 
26
#include <boost/format.hpp>
 
27
#include "vbutil.h"
 
28
#include "vbio.h"
 
29
#include "imageutils.h"
 
30
#include "vbversion.h"
 
31
#include "vbimagemunge.hlp.h"
 
32
 
 
33
void vbimagemunge_help();
 
34
void vbimagemunge_version();
 
35
 
 
36
int cube_combine(Cube &cb,tokenlist &args);
 
37
void selectcubes(list<Cube> &clist,size_t n);
 
38
bool comparedims(list<Cube> cubelist);
 
39
void cube_random01(Cube &cb);
 
40
void cube_printregioninfo(Cube &cb);
 
41
void reallyload(Cube &cube);
 
42
void build_oplist();
 
43
 
 
44
// ugly globals -- code needs some reorg
 
45
Cube mycube;
 
46
int multi_index;
 
47
 
 
48
class imageop {
 
49
public:
 
50
  string name;
 
51
  int nargs;
 
52
  tokenlist args;
 
53
  // an operator is combining if either initfn or finishfn is non-null
 
54
  vbreturn (*initfn)(list<Cube> &cubelist,tokenlist &args);
 
55
  vbreturn (*procfn)(Cube &cube,tokenlist &args);
 
56
  vbreturn (*finishfn)(list<Cube> &cubelist,tokenlist &args);
 
57
  // storage for combining operators
 
58
  // sole constructor requires name and nargs
 
59
  imageop();
 
60
  imageop(string xname,int xnargs,vbreturn (*xprocfn)(Cube &,tokenlist &));
 
61
  void init(string xname,int xnargs,vbreturn (*xprocfn)(Cube &,tokenlist &));
 
62
  // random methods
 
63
  bool iscombining() {return (initfn!=NULL || finishfn!=NULL);}
 
64
};
 
65
 
 
66
imageop::imageop()
 
67
{
 
68
  init("",0,NULL);
 
69
}
 
70
 
 
71
imageop::imageop(string xname,int xnargs,vbreturn (*xprocfn)(Cube &,tokenlist &))
 
72
{
 
73
  init(xname,xnargs,xprocfn);
 
74
}
 
75
                 
 
76
 
 
77
void imageop::init(string xname,int xnargs,vbreturn (*xprocfn)(Cube &,tokenlist &))
 
78
{
 
79
  name=xname;
 
80
  nargs=xnargs;
 
81
  args.clear();
 
82
  initfn=NULL;
 
83
  procfn=xprocfn;
 
84
  finishfn=NULL;
 
85
}
 
86
 
 
87
// three operator lists that are handled differently
 
88
list<imageop> phase1;      // do the cubes one at a time.  if last operator is combining, feed.
 
89
list<imageop> phase2;      // do each operator for all cubes
 
90
typedef list<imageop>::iterator OPI;
 
91
typedef list<Cube>::iterator CUBI;
 
92
 
 
93
map<string,imageop> oplist;
 
94
 
 
95
int
 
96
main(int argc,char *argv[])
 
97
{
 
98
  if (argc==1) {
 
99
    vbimagemunge_help();
 
100
    exit(0);
 
101
  }
 
102
  build_oplist();
 
103
  tokenlist args;
 
104
  list<string> filelist;
 
105
  // string outfile;
 
106
  args.Transfer(argc-1,argv+1);
 
107
  bool inops=0;
 
108
  int phase=1;
 
109
  for (int i=0; i<args.size(); i++) {
 
110
    if (inops) {
 
111
      if (oplist.count(args[i])==0) {
 
112
        printf("[E] vbimagemunge: unknown operator %s\n",args(i));
 
113
        exit(202);
 
114
      }
 
115
      imageop myop=oplist[args[i]];
 
116
      if (i+myop.nargs>=args.size()) {
 
117
        printf("[E] vbimagemunge: operator %s requires %d arguments\n",
 
118
               myop.name.c_str(),myop.nargs);
 
119
        exit(203);
 
120
      }
 
121
      myop.args.Add(args[i]);  // first arg is the operator name
 
122
      for (int j=0; j<myop.nargs; j++)
 
123
        myop.args.Add(args[++i]);
 
124
      // now decide which list gets the op
 
125
      if (phase==2)
 
126
        phase2.push_back(myop);
 
127
      // if we're in phase 1 and it's parallel, just put it on
 
128
      else if (!myop.iscombining())
 
129
        phase1.push_back(myop);
 
130
      // if we're in phase 1 and it's combining, see if we have any parallels first
 
131
      else {
 
132
        if (phase1.size())
 
133
          phase1.push_back(myop);
 
134
        else
 
135
          phase2.push_back(myop);
 
136
        phase=2;
 
137
      }
 
138
    }
 
139
    else {
 
140
      if (args[i]=="--") {
 
141
        inops=1;
 
142
      }
 
143
      else if (args[i]=="-h") {
 
144
        vbimagemunge_help();
 
145
        exit(0);
 
146
      }
 
147
      else if (args[i]=="-v") {
 
148
        vbimagemunge_version();
 
149
        exit(0);
 
150
      }
 
151
      else
 
152
        filelist.push_back(args[i]);
 
153
    }
 
154
  }
 
155
 
 
156
  list<Cube> cubelist;
 
157
  Cube master;
 
158
 
 
159
  // pre-load all the files.  we should figure out how not to load all
 
160
  // the image data from 4D files, maybe.
 
161
 
 
162
  for (list<string>::iterator ff=filelist.begin(); ff!=filelist.end(); ff++) {
 
163
    printf("[I] vbimagemunge: reading file %s\n",ff->c_str());
 
164
    Cube cb;
 
165
    Tes ts;
 
166
    if (cb.ReadHeader(*ff)==0) {
 
167
      cubelist.push_back(cb);
 
168
      continue;
 
169
    }
 
170
    if (ts.ReadFile(*ff)) {
 
171
      printf("[E] vbimagemunge: couldn't read file %s, continuing anyway\n",ff->c_str());
 
172
      continue;
 
173
    }
 
174
    for (int i=0; i<ts.dimt; i++) {
 
175
      Cube tmpc;
 
176
      cubelist.push_back(tmpc);
 
177
      ts.getCube(i,cubelist.back());
 
178
    }
 
179
  }
 
180
  if (cubelist.size()==0) {
 
181
    printf("[E] vbimagemunge: no valid input cubes found\n");
 
182
    exit(11);
 
183
  }
 
184
  
 
185
  // first do phase 1
 
186
  int index=0;
 
187
  if (phase1.size()) {
 
188
    if (phase1.back().iscombining() && phase1.back().initfn)
 
189
      phase1.back().initfn(cubelist,phase1.back().args);
 
190
    for (CUBI cc=cubelist.begin(); cc!=cubelist.end(); cc++) {
 
191
      for (OPI oo=phase1.begin(); oo!=phase1.end(); oo++) {
 
192
        if (oo->procfn)
 
193
          oo->procfn(*cc,oo->args);
 
194
      }
 
195
      index++;
 
196
    }
 
197
    if (phase1.back().iscombining() && phase1.back().finishfn)
 
198
      phase1.back().finishfn(cubelist,phase1.back().args);
 
199
  }
 
200
  // if we have a phase2 block, now do that
 
201
  for (OPI oo=phase2.begin(); oo!=phase2.end(); oo++) {
 
202
    if (oo->initfn)
 
203
      oo->initfn(cubelist,oo->args);
 
204
    index=0;
 
205
    for (CUBI cc=cubelist.begin(); cc!=cubelist.end(); cc++) {
 
206
      if (oo->procfn)
 
207
        oo->procfn(*cc,oo->args);
 
208
      index++;
 
209
    }
 
210
    if (oo->finishfn)
 
211
      oo->finishfn(cubelist,oo->args);
 
212
  }
 
213
  
 
214
  exit(0);
 
215
}
 
216
 
 
217
// for "combine" the arguments are:
 
218
// combine fixed/moving x y z op value
 
219
 
 
220
template<class T>
 
221
int
 
222
cube_combine_fixed(Cube &cb,tokenlist &args)
 
223
{
 
224
  int xx=strtol(args[2]);
 
225
  int yy=strtol(args[3]);
 
226
  int zz=strtol(args[4]);
 
227
  string method=args[5];
 
228
  double value=strtod(args[6]);
 
229
  int f_needval=0;
 
230
  if (method=="sum" || method=="average" || method=="sumthresh" || method=="averagethresh")
 
231
    f_needval=1;
 
232
 
 
233
  if (xx<1 || yy<1 || zz<1) {
 
234
    printf("[E] vbimagemunge: invalid subcube dimensions\n");
 
235
    return 101;
 
236
  }
 
237
 
 
238
  Cube cb2(cb);
 
239
  // initial subcube range 0:xx-1, etc.
 
240
  int x1=0,y1=0,z1=0;
 
241
  int xn=xx-1,yn=yy-1,zn=zz-1;
 
242
 
 
243
  int i,j,k,count;
 
244
  double total;
 
245
  T val;
 
246
  while(1) {
 
247
    count=0;
 
248
    total=0.0;
 
249
    if (x1>cb.dimx-1 || y1>cb.dimy-1 || z1>cb.dimz-1)
 
250
      break;
 
251
    if (xn>cb.dimx-1) xn=cb.dimx-1;
 
252
    if (yn>cb.dimy-1) yn=cb.dimy-1;
 
253
    if (zn>cb.dimz-1) zn=cb.dimz-1;
 
254
    for (i=x1; i<=xn; i++) {
 
255
      for (j=y1; j<=yn; j++) {
 
256
        for (k=z1; k<=zn; k++) {
 
257
          if (f_needval) {
 
258
            val=cb.getValue<T>(i,j,k);
 
259
            if (val>0.0)
 
260
              count++;
 
261
            total+=val;
 
262
          }
 
263
          else if (cb.testValue(i,j,k))
 
264
            count++;
 
265
        }
 
266
      }
 
267
    }
 
268
    // calculate the replacement value -- sum, count, average,
 
269
    // pct, any sumthresh, averagethresh
 
270
    T newvalue=0;
 
271
    int rsize=(xn-x1+1)*(yn-y1+1)*(zn-z1+1);
 
272
    if (method=="count")
 
273
      newvalue=count;
 
274
    else if (method=="sum")
 
275
      newvalue=total;
 
276
    else if (method=="pct")
 
277
      newvalue=(count*100)/rsize;
 
278
    else if (method=="average")
 
279
      newvalue=total/rsize;
 
280
    else if (method=="any") {
 
281
      if (count) newvalue=1;
 
282
    }
 
283
    else if (method=="sumthresh") {
 
284
      if (total>value) newvalue=1;
 
285
    }
 
286
    else if (method=="averagethresh") {
 
287
      if ((total/rsize)>value) newvalue=1;
 
288
    }
 
289
    else
 
290
      newvalue=0;
 
291
    // all voxels in the subcube get the identical new value
 
292
    for (i=x1; i<=xn; i++) {
 
293
      for (j=y1; j<=yn; j++) {
 
294
        for (k=z1; k<=zn; k++) {
 
295
          cb2.setValue<T>(i,j,k,newvalue);
 
296
        }
 
297
      }
 
298
    }
 
299
    // next subcube
 
300
    if (xn<cb.dimx-1) {
 
301
      x1+=xx;
 
302
      xn+=xx;
 
303
    }
 
304
    else if (yn<cb.dimy-1) {
 
305
      x1=0;
 
306
      xn=xx-1;
 
307
      y1+=yy;
 
308
      yn+=yy;
 
309
    }
 
310
    else if (zn<cb.dimz-1) {
 
311
      x1=y1=0;
 
312
      xn=xx-1;
 
313
      yn=yy-1;
 
314
      z1+=zz;
 
315
      zn+=zz;
 
316
    }
 
317
    else
 
318
      break;
 
319
  }
 
320
  cb=cb2;
 
321
  return 0;
 
322
}
 
323
 
 
324
template<class T>
 
325
int
 
326
cube_combine_moving(Cube &cb,tokenlist &args)
 
327
{
 
328
  int xx=strtol(args[2]);
 
329
  int yy=strtol(args[3]);
 
330
  int zz=strtol(args[4]);
 
331
  string method=args[5];
 
332
  double value=strtod(args[6]);
 
333
  int f_needval=0;
 
334
  if (method=="sum" || method=="average" || method=="sumthresh" || method=="averagethresh")
 
335
    f_needval=1;
 
336
 
 
337
  if (xx<1 || yy<1 || zz<1) {
 
338
    printf("[E] vbimagemunge: invalid subcube dimensions\n");
 
339
    return 101;
 
340
  }
 
341
 
 
342
  Cube cb2(cb);
 
343
  // initial subcube range 0:xx-1, etc.
 
344
  int x1=0,y1=0,z1=0;
 
345
  int xn=xx-1,yn=yy-1,zn=zz-1;
 
346
 
 
347
  int i,j,k,count;
 
348
  double total;
 
349
  T val;
 
350
 
 
351
  for (int xxx=0; xxx<=cb.dimx; xxx++) {
 
352
    for (int yyy=0; yyy<=cb.dimy; yyy++) {
 
353
      for (int zzz=0; zzz<=cb.dimz; zzz++) {
 
354
        count=0;
 
355
        total=0.0;
 
356
        // build a neighborhood around the voxel
 
357
        x1=xxx-xx/2;
 
358
        xn=x1+xx-1;
 
359
        y1=yyy-yy/2;
 
360
        yn=y1+yy-1;
 
361
        z1=zzz-zz/2;
 
362
        zn=z1+zz-1;
 
363
        // clip it to cube boundaries
 
364
        if (x1<0) x1=0;
 
365
        if (y1<0) y1=0;
 
366
        if (z1<0) z1=0;
 
367
        if (xn>cb.dimx-1) xn=cb.dimx-1;
 
368
        if (yn>cb.dimy-1) yn=cb.dimy-1;
 
369
        if (zn>cb.dimz-1) zn=cb.dimz-1;
 
370
        for (i=x1; i<=xn; i++) {
 
371
          for (j=y1; j<=yn; j++) {
 
372
            for (k=z1; k<=zn; k++) {
 
373
              if (f_needval) {
 
374
                val=cb.getValue<T>(i,j,k);
 
375
                if (val>0.0)
 
376
                  count++;
 
377
                total+=val;
 
378
              }
 
379
              else if (cb.testValue(i,j,k))
 
380
                count++;
 
381
            }
 
382
          }
 
383
        }
 
384
        // calculate the replacement value -- sum, count, average,
 
385
        // pct, any sumthresh, averagethresh
 
386
        T newvalue=0;
 
387
        int rsize=(xn-x1+1)*(yn-y1+1)*(zn-z1+1);
 
388
        if (method=="count")
 
389
          newvalue=count;
 
390
        else if (method=="sum")
 
391
          newvalue=total;
 
392
        else if (method=="pct")
 
393
          newvalue=(count*100)/rsize;
 
394
        else if (method=="average")
 
395
          newvalue=total/rsize;
 
396
        else if (method=="any") {
 
397
          if (count) newvalue=1;
 
398
        }
 
399
        else if (method=="sumthresh") {
 
400
          if (total>value) newvalue=1;
 
401
        }
 
402
        else if (method=="averagethresh") {
 
403
          if ((total/rsize)>value) newvalue=1;
 
404
        }
 
405
        else
 
406
          newvalue=0;
 
407
        // set just the original voxel
 
408
        cb2.setValue<T>(xxx,yyy,zzz,newvalue);
 
409
      }
 
410
    }
 
411
  }
 
412
  cb=cb2;
 
413
  return 0;
 
414
}
 
415
 
 
416
int
 
417
cube_combine(Cube &cb,tokenlist &args)
 
418
{
 
419
  if (args[1]=="fixed" || args[1]=="f") {
 
420
    switch(cb.datatype) {
 
421
    case vb_byte:
 
422
      return cube_combine_fixed<char>(cb,args);
 
423
      break;
 
424
    case vb_short:
 
425
      return cube_combine_fixed<int16>(cb,args);
 
426
      break;
 
427
    case vb_long:
 
428
      return cube_combine_fixed<int32>(cb,args);
 
429
      break;
 
430
    case vb_float:
 
431
      return cube_combine_fixed<float>(cb,args);
 
432
      break;
 
433
    case vb_double:
 
434
      return cube_combine_fixed<double>(cb,args);
 
435
      break;
 
436
    }
 
437
  }
 
438
  else if (args[1]=="moving" || args[1]=="m") {
 
439
    switch(cb.datatype) {
 
440
    case vb_byte:
 
441
      return cube_combine_moving<char>(cb,args);
 
442
      break;
 
443
    case vb_short:
 
444
      return cube_combine_moving<int16>(cb,args);
 
445
      break;
 
446
    case vb_long:
 
447
      return cube_combine_moving<int32>(cb,args);
 
448
      break;
 
449
    case vb_float:
 
450
      return cube_combine_moving<float>(cb,args);
 
451
      break;
 
452
    case vb_double:
 
453
      return cube_combine_moving<double>(cb,args);
 
454
      break;
 
455
    }
 
456
  }
 
457
 
 
458
  return 101;  // shouldn't happen
 
459
}
 
460
 
 
461
void
 
462
cube_printregioninfo(Cube &cb)
 
463
{
 
464
  vector<VBRegion> rlist=findregions(cb,vb_ne,0.0);
 
465
  int totalvoxels=0;
 
466
  for (size_t i=0; i<rlist.size(); i++)
 
467
    totalvoxels+=rlist[i].size();
 
468
  printf("[I] vbimagemunge: %d voxels in %d regions found in volume %s\n",
 
469
         totalvoxels,(int)(rlist.size()),cb.GetFileName().c_str());
 
470
  for (size_t i=0; i<rlist.size(); i++) {
 
471
    // calculate centers of mass
 
472
    double x1=0,y1=0,z1=0,x2=0,y2=0,z2=0,totalmass=0;
 
473
    for (VI myvox=rlist[i].begin(); myvox!=rlist[i].end(); myvox++) {
 
474
      x1+=myvox->second.x;
 
475
      y1+=myvox->second.y;
 
476
      z1+=myvox->second.z;
 
477
      x2+=myvox->second.x*myvox->second.val;
 
478
      y2+=myvox->second.y*myvox->second.val;
 
479
      z2+=myvox->second.z*myvox->second.val;
 
480
      totalmass+=myvox->second.val;
 
481
    }
 
482
    x1/=rlist[i].size();
 
483
    y1/=rlist[i].size();
 
484
    z1/=rlist[i].size();
 
485
    x2/=totalmass;
 
486
    y2/=totalmass;
 
487
    z2/=totalmass;
 
488
    
 
489
    printf("[I] vbimagemunge:   region %d, %d voxels,",(int)i,(int)(rlist[i].size()));
 
490
    printf(" center of mass: %g,%g,%g (weighted: %g,%g,%g)\n",x1,y1,z1,x2,y2,z2);
 
491
    // totalvoxels+=rlist[i].size();    
 
492
  }
 
493
  // printf("[I] vbimagemunge:   %d total voxels in volume %s\n",totalvoxels,cb.GetFileName().c_str());
 
494
}
 
495
 
 
496
void
 
497
cube_random01(Cube &cb)
 
498
{
 
499
  uint32 rr=0;
 
500
  int pos=32;
 
501
  cb.zero();
 
502
  for (int i=0; i<cb.dimx; i++) {
 
503
    for (int j=0; j<cb.dimy; j++) {
 
504
      for (int k=0; k<cb.dimz; k++) {
 
505
        if (pos>31) {
 
506
          rr=VBRandom();
 
507
          pos=0;
 
508
        }
 
509
        if (rr & 1<<pos)
 
510
          cb.SetValue(i,j,k,1.0);
 
511
        pos++;
 
512
      }
 
513
    }
 
514
  }
 
515
}
 
516
 
 
517
void
 
518
selectcubes(list<Cube> &clist,size_t n)
 
519
{
 
520
  if (n>=clist.size() || n<1)
 
521
    return;
 
522
  map<uint32,bool> randlist;
 
523
  while (randlist.size()<n)
 
524
    randlist[VBRandom()]=1;
 
525
  while (randlist.size()<clist.size())
 
526
    randlist[VBRandom()]=0;
 
527
  map<uint32,bool>::iterator mm=randlist.begin();
 
528
  list<Cube>::iterator cc=clist.begin();
 
529
  int oldn=clist.size();
 
530
  for (int i=0; i<oldn; i++) {
 
531
    if (mm->second)
 
532
      cc++;
 
533
    else
 
534
      cc=clist.erase(cc);
 
535
    mm++;
 
536
  }
 
537
}
 
538
 
 
539
bool
 
540
comparedims(list<Cube> cubelist)
 
541
{
 
542
  if (cubelist.size()<2)
 
543
    return 1;
 
544
  list<Cube>::iterator cc=cubelist.begin();
 
545
  int dimx=cc->dimx;
 
546
  int dimy=cc->dimy;
 
547
  int dimz=cc->dimz;
 
548
  cc++;
 
549
  while (cc!=cubelist.end()) {
 
550
    if (cc->dimx!=dimx || cc->dimy!=dimy || cc->dimz!=dimz)
 
551
      return 0;
 
552
    cc++;
 
553
  }
 
554
  return 1;
 
555
}
 
556
 
 
557
// initfns
 
558
 
 
559
vbreturn op_shortzeros(list<Cube> &cubelist,tokenlist &)
 
560
{
 
561
  if (cubelist.size()==0)
 
562
    return 101;
 
563
  mycube.SetVolume(cubelist.front().dimx,cubelist.front().dimy,
 
564
                   cubelist.front().dimz,vb_short);
 
565
  return 0;
 
566
}
 
567
 
 
568
vbreturn op_allones(list<Cube> &cubelist,tokenlist &)
 
569
{
 
570
  if (cubelist.size()==0)
 
571
    return 101;
 
572
  reallyload(cubelist.front());
 
573
  mycube=cubelist.front();
 
574
  mycube.quantize(1.0);
 
575
  return 0;
 
576
}
 
577
 
 
578
vbreturn op_allzeros(list<Cube> &cubelist,tokenlist &)
 
579
{
 
580
  if (cubelist.size()==0)
 
581
    return 101;
 
582
  reallyload(cubelist.front());
 
583
  mycube=cubelist.front();
 
584
  mycube.zero();
 
585
  return 0;
 
586
}
 
587
 
 
588
vbreturn op_copy(list<Cube> &cubelist,tokenlist &)
 
589
{
 
590
  if (cubelist.size()==0)
 
591
    return 101;
 
592
  reallyload(cubelist.front());
 
593
  mycube=cubelist.front();
 
594
  return 0;
 
595
}
 
596
 
 
597
// procfn
 
598
 
 
599
vbreturn op_smoothvox(Cube &cube,tokenlist &args)
 
600
{
 
601
  reallyload(cube);
 
602
  smoothCube(cube,strtod(args[1]),strtod(args[2]),strtod(args[3]));
 
603
  return 0;
 
604
}
 
605
 
 
606
vbreturn op_smoothmm(Cube &cube,tokenlist &)
 
607
{
 
608
  reallyload(cube);
 
609
  printf("smoothmm not implemented\n");
 
610
  exit(5);
 
611
  return 0;
 
612
}
 
613
 
 
614
vbreturn op_thresh(Cube &cube,tokenlist &args)
 
615
{
 
616
  reallyload(cube);
 
617
  cube.thresh(strtod(args[1]));
 
618
  return 0;
 
619
}
 
620
 
 
621
vbreturn op_rotate(Cube &cube,tokenlist &args)
 
622
{
 
623
  reallyload(cube);
 
624
  rotatecube(cube,strtod(args[1]),strtod(args[2]),strtod(args[3]));
 
625
  return 0;
 
626
}
 
627
 
 
628
vbreturn op_regionat(Cube &cube,tokenlist &args)
 
629
{
 
630
  reallyload(cube);
 
631
  Cube mask;
 
632
  mask.SetVolume(cube.dimx,cube.dimy,cube.dimz,vb_short);
 
633
  mask+=1;
 
634
  VBRegion rr;
 
635
  rr=growregion(strtol(args[1]),strtol(args[2]),strtol(args[3]),cube,mask,vb_agt,FLT_MIN);
 
636
  mask*=0;
 
637
  for (VI myvox=rr.begin(); myvox!=rr.end(); myvox++) {
 
638
    mask.setValue(myvox->second.x,myvox->second.y,myvox->second.z,(int32)1);
 
639
  }
 
640
  cube=mask;
 
641
  return 0;
 
642
}
 
643
 
 
644
vbreturn op_splitregions(Cube &cube,tokenlist &args)
 
645
{
 
646
  reallyload(cube);
 
647
  vector<VBRegion> regions;
 
648
  regions=findregions(cube,vb_gt,0.0);
 
649
  Cube mask;
 
650
  mask.SetVolume(cube.dimx,cube.dimy,cube.dimz,vb_byte);
 
651
  int index=0;
 
652
  string fname;
 
653
  vbforeach(VBRegion &rr,regions) {
 
654
    mask.zero();
 
655
    for (VI myvox=rr.begin(); myvox!=rr.end(); myvox++)
 
656
      mask.setValue<char>(myvox->second.x,myvox->second.y,myvox->second.z,1);
 
657
    fname=args[1];
 
658
    string num=(format("%05d")%index).str();
 
659
    replace_string(fname,"XXX",num);
 
660
    mask.WriteFile(fname);
 
661
    index++;
 
662
  }
 
663
  cube=mask;
 
664
  return 0;
 
665
}
 
666
 
 
667
vbreturn op_threshabs(Cube &cube,tokenlist &args)
 
668
{
 
669
  reallyload(cube);
 
670
  cube.threshabs(strtod(args[1]));
 
671
  return 0;
 
672
}
 
673
 
 
674
vbreturn op_cutoff(Cube &cube,tokenlist &args)
 
675
{
 
676
  reallyload(cube);
 
677
  cube.cutoff(strtod(args[1]));
 
678
  return 0;
 
679
}
 
680
 
 
681
vbreturn op_orient(Cube &cube,tokenlist &)
 
682
{
 
683
  reallyload(cube);
 
684
  printf("orient not implemented\n");
 
685
  exit(5);
 
686
  return 0;
 
687
}
 
688
 
 
689
vbreturn op_quantize(Cube &cube,tokenlist &args)
 
690
{
 
691
  reallyload(cube);
 
692
  cube.quantize(strtod(args[1]));
 
693
  return 0;
 
694
}
 
695
 
 
696
vbreturn op_invert(Cube &cube,tokenlist &)
 
697
{
 
698
  reallyload(cube);
 
699
  cube.invert();
 
700
  return 0;
 
701
}
 
702
 
 
703
vbreturn op_flipx(Cube &cube,tokenlist &)
 
704
{
 
705
  reallyload(cube);
 
706
  cube.flipx();
 
707
  return 0;
 
708
}
 
709
 
 
710
vbreturn op_flipy(Cube &cube,tokenlist &)
 
711
{
 
712
  reallyload(cube);
 
713
  cube.flipy();
 
714
  return 0;
 
715
}
 
716
 
 
717
vbreturn op_flipz(Cube &cube,tokenlist &)
 
718
{
 
719
  reallyload(cube);
 
720
  cube.flipz();
 
721
  return 0;
 
722
}
 
723
 
 
724
vbreturn op_zeroleft(Cube &cube,tokenlist &)
 
725
{
 
726
  reallyload(cube);
 
727
  cube.rightify();
 
728
  return 0;
 
729
}
 
730
 
 
731
vbreturn op_zeroright(Cube &cube,tokenlist &)
 
732
{
 
733
  reallyload(cube);
 
734
  cube.leftify();
 
735
  return 0;
 
736
}
 
737
 
 
738
vbreturn op_bigendian(Cube &cube,tokenlist &)
 
739
{
 
740
  reallyload(cube);
 
741
  cube.filebyteorder=ENDIAN_BIG;
 
742
  return 0;
 
743
}
 
744
 
 
745
vbreturn op_littleendian(Cube &cube,tokenlist &)
 
746
{
 
747
  reallyload(cube);
 
748
  cube.filebyteorder=ENDIAN_LITTLE;
 
749
  return 0;
 
750
}
 
751
 
 
752
vbreturn op_byteswap(Cube &cube,tokenlist &)
 
753
{
 
754
  reallyload(cube);
 
755
  cube.byteswap();
 
756
  return 0;
 
757
}
 
758
 
 
759
vbreturn op_combine(Cube &cube,tokenlist &args)
 
760
{
 
761
  reallyload(cube);
 
762
  if (cube_combine(cube,args)) {
 
763
    printf("[E] vbimagemunge: error downsampling\n");
 
764
    exit(212);
 
765
  }
 
766
  return 0;
 
767
}
 
768
 
 
769
vbreturn op_convert(Cube &cube,tokenlist &args)
 
770
{
 
771
  reallyload(cube);
 
772
  Cube tmpc;
 
773
  string newtype=vb_tolower(args[1]);
 
774
  int err=1;
 
775
  if (newtype=="byte")
 
776
    err=cube.convert_type(vb_byte,VBSETALT);
 
777
  if (newtype=="int16")
 
778
    err=cube.convert_type(vb_short,VBSETALT);
 
779
  if (newtype=="int32")
 
780
    err=cube.convert_type(vb_long,VBSETALT);
 
781
  if (newtype=="float")
 
782
    err=cube.convert_type(vb_float,VBSETALT|VBNOSCALE);
 
783
  if (newtype=="double")
 
784
    err=cube.convert_type(vb_double,VBSETALT|VBNOSCALE);
 
785
 
 
786
  if (err) {
 
787
    printf("[E] vbimagemunge: error converting datatype\n");
 
788
    exit(212);
 
789
  }
 
790
  cube=tmpc;
 
791
  return 0;
 
792
}
 
793
 
 
794
vbreturn op_removesmallregions(Cube &cube,tokenlist &)
 
795
{
 
796
  reallyload(cube);
 
797
  return 0;
 
798
}
 
799
 
 
800
vbreturn op_add(Cube &cube,tokenlist &args)
 
801
{
 
802
  reallyload(cube);
 
803
  if (vb_fileexists(args[1])) {
 
804
    Cube tmp;
 
805
    if (tmp.ReadFile(args[1])) {
 
806
      printf("[E] vbimagemunge: couldn't read file %s\n",args(1));
 
807
      exit(202);
 
808
    }
 
809
    cube+=tmp;
 
810
  }
 
811
  else
 
812
    cube+=strtod(args[1]);
 
813
  return 0;
 
814
}
 
815
 
 
816
vbreturn op_sub(Cube &cube,tokenlist &args)
 
817
{
 
818
  reallyload(cube);
 
819
  if (vb_fileexists(args[1])) {
 
820
    Cube tmp;
 
821
    if (tmp.ReadFile(args[1])) {
 
822
      printf("[E] vbimagemunge: couldn't read file %s\n",args(1));
 
823
      exit(202);
 
824
    }
 
825
    cube-=tmp;
 
826
  }
 
827
  else
 
828
    cube-=strtod(args[1]);
 
829
  return 0;
 
830
}
 
831
 
 
832
vbreturn op_mult(Cube &cube,tokenlist &args)
 
833
{
 
834
  reallyload(cube);
 
835
  if (vb_fileexists(args[1])) {
 
836
    Cube tmp;
 
837
    if (tmp.ReadFile(args[1])) {
 
838
      printf("[E] vbimagemunge: couldn't read file %s\n",args(1));
 
839
      exit(202);
 
840
    }
 
841
    cube*=tmp;
 
842
  }
 
843
  else
 
844
    cube*=strtod(args[1]);
 
845
  return 0;
 
846
}
 
847
 
 
848
vbreturn op_div(Cube &cube,tokenlist &args)
 
849
{
 
850
  reallyload(cube);
 
851
  if (vb_fileexists(args[1])) {
 
852
    Cube tmp;
 
853
    if (tmp.ReadFile(args[1])) {
 
854
      printf("[E] vbimagemunge: couldn't read file %s\n",args(1));
 
855
      exit(202);
 
856
    }
 
857
    cube/=tmp;
 
858
  }
 
859
  else
 
860
    cube/=strtod(args[1]);
 
861
  return 0;
 
862
}
 
863
 
 
864
vbreturn op_nminus(Cube &cube,tokenlist &args)
 
865
{
 
866
  reallyload(cube);
 
867
  if (vb_fileexists(args[1])) {
 
868
    Cube tmp;
 
869
    if (tmp.ReadFile(args[1])) {
 
870
      printf("[E] vbimagemunge: couldn't read file %s\n",args(1));
 
871
      exit(202);
 
872
    }
 
873
    cube=tmp-cube;
 
874
  }
 
875
  else {
 
876
    cube-=strtod(args[1]);
 
877
    cube*=-1.0;
 
878
  }
 
879
  return 0;
 
880
}
 
881
 
 
882
vbreturn op_random01(Cube &cube,tokenlist &)
 
883
{
 
884
  reallyload(cube);
 
885
  cube_random01(cube);
 
886
  return 0;
 
887
}
 
888
 
 
889
vbreturn op_write(Cube &cube,tokenlist &args)
 
890
{
 
891
  reallyload(cube);
 
892
  if (cube.WriteFile(args[1]))
 
893
    printf("[E] vbimagemunge: error writing file %s\n",args(1));
 
894
  else
 
895
    printf("[I] vbimagemunge: wrote file %s\n",args(1));
 
896
  return 0;
 
897
}
 
898
 
 
899
vbreturn op_writeprefixed(Cube &cube,tokenlist &args)
 
900
{
 
901
  reallyload(cube);
 
902
  string xd=xdirname(cube.GetFileName());
 
903
  string xf=xfilename(cube.GetFileName()); 
 
904
  string fname=xd+"/"+args[1]+xf;
 
905
  if (xd==".")
 
906
    fname=args[1]+xf;
 
907
  if (cube.WriteFile(fname))
 
908
    printf("[E] vbimagemunge: error writing file %s\n",fname.c_str());
 
909
  else
 
910
    printf("[I] vbimagemunge: wrote file %s\n",fname.c_str());
 
911
  return 0;
 
912
}
 
913
 
 
914
vbreturn op_regioninfo(Cube &cube,tokenlist &)
 
915
{
 
916
  reallyload(cube);
 
917
  cube_printregioninfo(cube);
 
918
  return 0;
 
919
}
 
920
 
 
921
vbreturn op_remap(Cube &cube,tokenlist &args)
 
922
{
 
923
  reallyload(cube);
 
924
  ifstream fs;
 
925
  fs.open(args(1),ios::in);
 
926
  if (!fs) {
 
927
    printf("[E] vbimagemunge: couldn't read map file %s\n",args(1));
 
928
    return 1;
 
929
  }
 
930
  const short bufsz=1024;
 
931
  char buf[bufsz];
 
932
  map<int,float> mymap;
 
933
  int fromval;
 
934
  float toval;
 
935
  tokenlist mapargs;
 
936
  while (fs.good()) {
 
937
    fs.getline(buf,bufsz);
 
938
    if (fs.gcount()==0) break;
 
939
    mapargs.ParseLine(buf);
 
940
    if (mapargs.size()!=2) continue;
 
941
    fromval=strtol(mapargs[0]);
 
942
    toval=strtod(mapargs[1]);
 
943
    if (mymap.count(fromval)) {
 
944
      printf("[E] vbimagemunge: duplicate value %d found in map %s\n",
 
945
             fromval,args(1));
 
946
      fs.close();
 
947
      return 1;
 
948
    }
 
949
    mymap[fromval]=toval;
 
950
  }
 
951
  fs.close();
 
952
 
 
953
  // now build the new cube
 
954
  Cube newcube;
 
955
  newcube.SetVolume(cube.dimx,cube.dimy,cube.dimz,vb_float);
 
956
  int mappedcnt=0,unmappedcnt=0;
 
957
  for (int i=0; i<cube.dimx; i++) {
 
958
    for (int j=0; j<cube.dimy; j++) {
 
959
      for (int k=0; k<cube.dimz; k++) {
 
960
        int key=cube.getValue<int32>(i,j,k);
 
961
        if (mymap.count(key)) {
 
962
          newcube.SetValue(i,j,k,mymap[key]);
 
963
          mappedcnt++;
 
964
        }
 
965
        else {
 
966
          newcube.SetValue(i,j,k,key);
 
967
          unmappedcnt++;
 
968
        }
 
969
      }
 
970
    }
 
971
  }
 
972
  cube=newcube;
 
973
  printf("[I] vbimagemunge: applied map file %s (%d voxels mapped, %d unmapped)\n",
 
974
         args(1),mappedcnt,unmappedcnt);
 
975
  return 0;
 
976
}
 
977
 
 
978
vbreturn op_info(Cube &cube,tokenlist &)
 
979
{
 
980
  reallyload(cube);
 
981
  printf("[I] cube %s: min=%f max=%f infs/nans=%d\n",
 
982
         cube.GetFileName().c_str(),
 
983
         cube.get_minimum(),cube.get_maximum(),(int)cube.get_nonfinites());
 
984
  return 0;
 
985
}
 
986
 
 
987
// procfns for combining operations
 
988
 
 
989
vbreturn op_sum(Cube &cube,tokenlist &)
 
990
{
 
991
  reallyload(cube);
 
992
  mycube+=cube;
 
993
  return 0;
 
994
}
 
995
 
 
996
vbreturn op_product(Cube &cube,tokenlist &)
 
997
{
 
998
  reallyload(cube);
 
999
  mycube*=cube;
 
1000
  return 0;
 
1001
}
 
1002
 
 
1003
vbreturn op_multi(Cube &cube,tokenlist &)
 
1004
{
 
1005
  reallyload(cube);
 
1006
  cube.quantize(1.0);
 
1007
  cube.invert();
 
1008
  mycube*=cube;
 
1009
  cube.invert();
 
1010
  cube*=multi_index;
 
1011
  mycube+=cube;
 
1012
  multi_index++;
 
1013
  return 0;
 
1014
}
 
1015
 
 
1016
vbreturn op_union(Cube &cube,tokenlist &)
 
1017
{
 
1018
  reallyload(cube);
 
1019
  cube.quantize(1.0);
 
1020
  mycube.unionmask(cube);
 
1021
  return 0;
 
1022
}
 
1023
 
 
1024
vbreturn op_intersect(Cube &cube,tokenlist &)
 
1025
{
 
1026
  reallyload(cube);
 
1027
  mycube.intersect(cube);
 
1028
  return 0;
 
1029
}
 
1030
 
 
1031
vbreturn op_count(Cube &cube,tokenlist &)
 
1032
{
 
1033
  reallyload(cube);
 
1034
  cube.quantize(1.0);
 
1035
  mycube+=cube;
 
1036
  return 0;
 
1037
}
 
1038
 
 
1039
// finishfns for combining operations
 
1040
 
 
1041
vbreturn op_overlapmask(list <Cube> &cubelist,tokenlist &)
 
1042
{
 
1043
  if (cubelist.size()!=2 && cubelist.size()!=3) {
 
1044
    printf("[E] vbimagemunge: overlapmask can only work with 2 or 3 input masks\n");
 
1045
    return 101;
 
1046
  }
 
1047
  CUBI cc1,cc2,cc3;
 
1048
  cc1=cubelist.begin();
 
1049
  cc2=cc1; cc2++;
 
1050
  cc3=cc2; cc3++;
 
1051
  reallyload(*cc1);
 
1052
  reallyload(*cc2);
 
1053
  cc1->filename=xsetextension(xfilename(cc1->filename),"");
 
1054
  cc2->filename=xsetextension(xfilename(cc2->filename),"");
 
1055
  if (cubelist.size()==3) {
 
1056
    reallyload(*cc3);
 
1057
    cc3->filename=xsetextension(xfilename(cc3->filename),"");
 
1058
  }
 
1059
 
 
1060
  // check dims
 
1061
  int dimx=cc1->dimx;
 
1062
  int dimy=cc1->dimy;
 
1063
  int dimz=cc1->dimz;
 
1064
  for (CUBI cc=++(cubelist.begin()); cc!=cubelist.end(); cc++) {
 
1065
    //  for (size_t j=1; j<cubelist.size(); j++) {
 
1066
    if (cc->dimx!=dimx || cc->dimy!=dimy || cc->dimz!=dimz) {
 
1067
      printf("[E] vbimagemunge: overlapmask requires masks of matching dimensions\n");
 
1068
      return 101;
 
1069
    }
 
1070
  }
 
1071
 
 
1072
  Cube outmask;
 
1073
  outmask.SetVolume(dimx,dimy,dimz,vb_byte);
 
1074
  uint8 val=0;
 
1075
  for (int i=0; i<dimx*dimy*dimz; i++) {
 
1076
    val=0;
 
1077
    if (cc1->testValue(i)) val+=1;
 
1078
    if (cc2->testValue(i)) val+=2;
 
1079
    if (cubelist.size()==3)
 
1080
      if (cc3->testValue(i)) val+=4;
 
1081
    outmask.setValue<char>(i,val);
 
1082
  }
 
1083
  outmask.AddHeader((string)"vb_maskspec: 1 255 0 0 "+cc1->filename);
 
1084
  outmask.AddHeader((string)"vb_maskspec: 2 0 0 255 "+cc2->filename);
 
1085
  outmask.AddHeader((string)"vb_maskspec: 3 170 0 170 \""+cc1->filename+"/"+cc2->filename+"\"");
 
1086
  if (cubelist.size()==3) {
 
1087
    outmask.AddHeader((string)"vb_maskspec: 4 220 220 0 "+cc3->filename);
 
1088
    outmask.AddHeader((string)"vb_maskspec: 5 255 140 0 \""+cc1->filename+"/"+cc3->filename+"\"");
 
1089
    outmask.AddHeader((string)"vb_maskspec: 6 0 255 0 \""+cc2->filename+"/"+cc3->filename+"\"");
 
1090
    outmask.AddHeader((string)"vb_maskspec: 7 140 140 80\""+cc1->filename+"/"+cc2->filename+"/"+cc3->filename+"\"");
 
1091
  }
 
1092
  cubelist.clear();
 
1093
  cubelist.push_back(outmask);
 
1094
  return 0;
 
1095
}
 
1096
 
 
1097
vbreturn op_null(list <Cube> &cubelist,tokenlist &)
 
1098
{
 
1099
  cubelist.clear();
 
1100
  cubelist.push_back(mycube);
 
1101
  return 0;
 
1102
}
 
1103
 
 
1104
vbreturn op_scale(list<Cube> &cubelist,tokenlist &)
 
1105
{
 
1106
  if (cubelist.size()==0)
 
1107
    return 101;
 
1108
  mycube/=cubelist.size();
 
1109
  cubelist.clear();
 
1110
  cubelist.push_back(mycube);
 
1111
  return 0;
 
1112
}
 
1113
 
 
1114
vbreturn op_select(list<Cube> &cubelist,tokenlist &args)
 
1115
{
 
1116
  selectcubes(cubelist,strtol(args[1]));
 
1117
  return 0;
 
1118
}
 
1119
 
 
1120
vbreturn op_include(list<Cube> &cubelist,tokenlist &args)
 
1121
{
 
1122
  vector<int> nlist=numberlist(args[1]);
 
1123
  bool mymap[cubelist.size()];
 
1124
  memset(mymap,0,cubelist.size()*sizeof(bool));
 
1125
  for (size_t i=0; i<nlist.size(); i++)
 
1126
    if (nlist[i]<(int)cubelist.size())
 
1127
      mymap[nlist[i]]=1;
 
1128
  list<Cube>::iterator cc=cubelist.begin();
 
1129
  int osize=cubelist.size();
 
1130
  for (int i=0; i<osize; i++) {
 
1131
    if (!mymap[i])
 
1132
      cc=cubelist.erase(cc);
 
1133
    else
 
1134
      cc++;
 
1135
  }
 
1136
  return 0;
 
1137
}
 
1138
 
 
1139
vbreturn op_exclude(list<Cube> &cubelist,tokenlist &args)
 
1140
{
 
1141
  vector<int> nlist=numberlist(args[1]);
 
1142
  bool mymap[cubelist.size()];
 
1143
  memset(mymap,0,cubelist.size()*sizeof(bool));
 
1144
  for (size_t i=0; i<nlist.size(); i++)
 
1145
    if (nlist[i]<(int)cubelist.size())
 
1146
      mymap[nlist[i]]=1;
 
1147
  list<Cube>::iterator cc=cubelist.begin();
 
1148
  int osize=cubelist.size();
 
1149
  for (int i=0; i<osize; i++) {
 
1150
    if (mymap[i])
 
1151
      cc=cubelist.erase(cc);
 
1152
    else
 
1153
      cc++;    
 
1154
  }
 
1155
  return 0;
 
1156
}
 
1157
 
 
1158
vbreturn op_write4d(list<Cube> &cubelist,tokenlist &args)
 
1159
{
 
1160
  if (cubelist.size()==0)
 
1161
    return 101;
 
1162
  Tes ts;
 
1163
  Cube *cb=&(cubelist.front());
 
1164
  ts.SetVolume(cb->dimx,cb->dimy,cb->dimz,cubelist.size(),cb->datatype);
 
1165
  int ind=0;
 
1166
  for (CUBI cc=cubelist.begin(); cc!=cubelist.end(); cc++) {
 
1167
    reallyload(*cc);
 
1168
    ts.SetCube(ind++,*cc);
 
1169
    cc->invalidate();
 
1170
  }
 
1171
  if (ts.WriteFile(args[1]))
 
1172
    printf("[E] vbimagemunge: couldn't write file %s\n",args[1].c_str());
 
1173
  else
 
1174
    printf("[I] vbimagemunge: wrote file %s\n",args[1].c_str());
 
1175
  return 0;
 
1176
}
 
1177
 
 
1178
void
 
1179
build_oplist()
 
1180
{
 
1181
  // non-combining ones are easy  
 
1182
  oplist["smoothvox"]=imageop("smoothvox",3,op_smoothvox);
 
1183
  oplist["smoothmm"]=imageop("smoothmm",3,op_smoothmm);
 
1184
  oplist["thresh"]=imageop("thresh",1,op_thresh);
 
1185
  oplist["rotate"]=imageop("rotate",3,op_rotate);
 
1186
  oplist["regionat"]=imageop("regionat",3,op_regionat);
 
1187
  oplist["threshabs"]=imageop("threshabs",1,op_threshabs);
 
1188
  oplist["cutoff"]=imageop("cutoff",1,op_cutoff);
 
1189
  oplist["invert"]=imageop("invert",0,op_invert);
 
1190
  oplist["splitregions"]=imageop("splitregions",1,op_splitregions);
 
1191
 
 
1192
  oplist["flipx"]=imageop("flipx",0,op_flipx);
 
1193
  oplist["flipy"]=imageop("flipy",0,op_flipy);
 
1194
  oplist["flipz"]=imageop("flipz",0,op_flipz);
 
1195
  oplist["zeroleft"]=imageop("zeroleft",0,op_zeroleft);
 
1196
  oplist["zeroright"]=imageop("zeroright",0,op_zeroright);
 
1197
  oplist["bigendian"]=imageop("bigendian",0,op_bigendian);
 
1198
  oplist["littleendian"]=imageop("zeroright",0,op_littleendian);
 
1199
  oplist["byteswap"]=imageop("zeroright",0,op_byteswap);
 
1200
 
 
1201
  oplist["orient"]=imageop("",99,op_orient);
 
1202
  oplist["quantize"]=imageop("quantize",1,op_quantize);
 
1203
  oplist["remap"]=imageop("remap",1,op_remap);
 
1204
  oplist["combine"]=imageop("combine",6,op_combine);
 
1205
  oplist["convert"]=imageop("convert",1,op_convert);
 
1206
  oplist["removesmallregions"]=imageop("removesmallregions",1,op_removesmallregions);
 
1207
  oplist["add"]=imageop("add",1,op_add);
 
1208
  oplist["sub"]=imageop("sub",1,op_sub);
 
1209
  oplist["mult"]=imageop("mult",1,op_mult);
 
1210
  oplist["div"]=imageop("div",1,op_div);
 
1211
  oplist["nminus"]=imageop("nminus",1,op_nminus);
 
1212
  oplist["random01"]=imageop("random01",0,op_random01);
 
1213
  oplist["write"]=imageop("write",1,op_write);
 
1214
  oplist["writeprefixed"]=imageop("writeprefixed",1,op_writeprefixed);
 
1215
  oplist["regioninfo"]=imageop("regioninfo",0,op_regioninfo);
 
1216
  oplist["info"]=imageop("info",0,op_info);
 
1217
 
 
1218
  // combining operators take a little more work
 
1219
  imageop tmp;
 
1220
 
 
1221
  // by doing it this way, all the cubes are in core before we select.
 
1222
  // if we create an initfn that stores the inclusion map and a procfn
 
1223
  // that invalidates the data for the ones we don't want, we can be
 
1224
  // tighter about memory
 
1225
 
 
1226
  tmp=imageop("overlapmask",0,NULL);
 
1227
  tmp.finishfn=op_overlapmask;  
 
1228
  oplist[tmp.name]=tmp;
 
1229
 
 
1230
  tmp=imageop("select",1,NULL);
 
1231
  tmp.finishfn=op_select;  
 
1232
  oplist[tmp.name]=tmp;
 
1233
 
 
1234
  tmp=imageop("include",1,NULL);
 
1235
  tmp.finishfn=op_include;  
 
1236
  oplist[tmp.name]=tmp;
 
1237
 
 
1238
  tmp=imageop("exclude",1,NULL);
 
1239
  tmp.finishfn=op_exclude;  
 
1240
  oplist[tmp.name]=tmp;
 
1241
 
 
1242
  tmp=imageop("sum",0,NULL);
 
1243
  tmp.initfn=op_allzeros;
 
1244
  tmp.procfn=op_sum;
 
1245
  tmp.finishfn=op_null;
 
1246
  oplist[tmp.name]=tmp;
 
1247
 
 
1248
  tmp=imageop("average",0,NULL);
 
1249
  tmp.initfn=op_allzeros;
 
1250
  tmp.procfn=op_sum;
 
1251
  tmp.finishfn=op_scale;
 
1252
  oplist[tmp.name]=tmp;
 
1253
 
 
1254
  tmp=imageop("multi",0,NULL);
 
1255
  tmp.initfn=op_shortzeros;
 
1256
  tmp.procfn=op_multi;
 
1257
  tmp.finishfn=op_null;
 
1258
  oplist[tmp.name]=tmp;
 
1259
 
 
1260
  tmp=imageop("union",0,NULL);
 
1261
  tmp.initfn=op_shortzeros;
 
1262
  tmp.procfn=op_union;
 
1263
  tmp.finishfn=op_null;
 
1264
  oplist[tmp.name]=tmp;
 
1265
 
 
1266
  tmp=imageop("intersect",0,NULL);
 
1267
  tmp.initfn=op_allones;
 
1268
  tmp.procfn=op_intersect;
 
1269
  tmp.finishfn=op_null;
 
1270
  oplist[tmp.name]=tmp;
 
1271
 
 
1272
  tmp=imageop("count",0,NULL);
 
1273
  tmp.initfn=op_shortzeros;
 
1274
  tmp.procfn=op_count;
 
1275
  tmp.finishfn=op_null;
 
1276
  oplist[tmp.name]=tmp;
 
1277
 
 
1278
  tmp=imageop("product",0,NULL);
 
1279
  tmp.initfn=op_allones;
 
1280
  tmp.procfn=op_product;
 
1281
  tmp.finishfn=op_null;
 
1282
  oplist[tmp.name]=tmp;
 
1283
 
 
1284
  tmp=imageop("write4d",1,NULL);
 
1285
  tmp.finishfn=op_write4d;
 
1286
  oplist[tmp.name]=tmp;
 
1287
}
 
1288
 
 
1289
 
 
1290
void reallyload(Cube &cube)
 
1291
{
 
1292
  if (cube.data) return;
 
1293
  if (cube.ReadData(cube.GetFileName()))
 
1294
    printf("[E] vbimagemunge: error reading file %s\n",cube.GetFileName().c_str());
 
1295
}
 
1296
 
 
1297
void
 
1298
vbimagemunge_help()
 
1299
{
 
1300
  cout << boost::format(myhelp) % vbversion;
 
1301
}
 
1302
 
 
1303
void
 
1304
vbimagemunge_version()
 
1305
{
 
1306
  printf("VoxBo vbimagemunge (v%s)\n",vbversion.c_str());
 
1307
}