← Back to branch summary

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

« back to all changes in this revision

Viewing changes to crunch/resample.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:
crunch/resample.cpp
 
1
 
 
2
// resample.cpp
 
3
// sinc resampling of volumes, adapted from spm
 
4
// Copyright (c) 1998-2004 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
// note that the guts of the underlying sinc resample code were ported
 
26
// from SPM96b source code written for MATLAB by John Ashburner
 
27
 
 
28
using namespace std;
 
29
 
 
30
#include "vbcrunch.h"
 
31
 
 
32
class Resample {
 
33
private:
 
34
  string imagename;     // name of the image to be normalized
 
35
  string refname;       // name of the reference image
 
36
  string outname;       // name of the output image
 
37
 
 
38
  vector<string> newheaders;  // headers to be added to output file
 
39
 
 
40
  enum {vb_3d,vb_4d} mode;
 
41
  enum {vb_nn,vb_sinc} method;
 
42
  tokenlist args;                // structure to be used for argument parsing
 
43
  CrunchCube *mycube,*newcube;
 
44
  Tes *mytes,*newtes;
 
45
  double zsize;             // size of slices in mm
 
46
  double x1,y1,z1;          // start voxel (0 indexed) in x,y,z
 
47
  int nx,ny,nz;             // number of voxels in resampled image
 
48
  int dimx,dimy,dimz;       // original dimensions
 
49
  double voxsize[3];        // original voxsize
 
50
  double xstep,ystep,zstep; // resample interval in voxels
 
51
public:
 
52
  Resample(int argc,char **argv);
 
53
  ~Resample();
 
54
  void init();
 
55
  int Setup();           // set up parameters
 
56
  int ResampleFile();    // detect filetype
 
57
  int Resample3D();      // set up for 3d
 
58
  int Resample4D();      // set up for 4d
 
59
 
 
60
  int ResampleCube();       // dispatch to sinc or nn resample
 
61
  int SincResampleCube();   // do it for each cube
 
62
  int NNResampleCube();     // do it for each cube
 
63
  void AdjustCornerAndOrigin(VBImage &im);
 
64
 
 
65
  // methods for working out the resample parameters
 
66
  int UseZ(const string &refname,double zsize);
 
67
  int UseDims(const string &refname);
 
68
  int UseTLHC(const string &refname);
 
69
  int UseCorner(const string &refname);
 
70
  int SetAlign(const string mode,const string ref);
 
71
 
 
72
};
 
73
 
 
74
void resample_help();
 
75
 
 
76
int
 
77
main(int argc,char *argv[])
 
78
{
 
79
  tzset();                     // make sure all times are timezone corrected
 
80
  if (argc < 2) {              // not enough args, display autodocs
 
81
    resample_help();
 
82
    exit(0);
 
83
  }
 
84
 
 
85
  Resample *r=new Resample(argc,argv);    // init resample object, parse args
 
86
  if (!r) {
 
87
    printf("*** resample: couldn't allocate a tiny structure\n");
 
88
    exit(5);
 
89
  }
 
90
  int err;
 
91
  err=r->Setup();              // set the start, stop, nslices, etc.
 
92
  if (!err)
 
93
    err = r->ResampleFile();   // do it
 
94
  delete r;                    // clean up
 
95
  if (err)
 
96
    printf("*** resample: error resampling (%d).\n",err);
 
97
  else
 
98
    printf("resample: done.\n");
 
99
  exit(err);
 
100
}
 
101
 
 
102
void
 
103
Resample::init()
 
104
{
 
105
  imagename=refname="";
 
106
  zsize=0.0;
 
107
  x1=y1=z1=0;
 
108
  xstep=ystep=zstep=1.0;
 
109
  nx=41;
 
110
  ny=51;
 
111
  nz=27;
 
112
  method=vb_sinc;
 
113
}
 
114
 
 
115
Resample::Resample(int argc,char *argv[])
 
116
{
 
117
  init();
 
118
  args.Transfer(argc-1,argv+1);
 
119
}
 
120
 
 
121
int
 
122
Resample::Setup()
 
123
{
 
124
  if (args.size() < 2) {
 
125
    printf("*** resample: requires at least an input and an output filename\n");
 
126
    return 110;
 
127
  }
 
128
  switch(EligibleFileTypes(args[0])[0].getDimensions()) {
 
129
  case 4:
 
130
    mytes = new Tes;
 
131
    mytes->ReadFile(args[0]);
 
132
    if (!mytes->data_valid) {
 
133
      printf("*** resample: invalid 4D file: %s\n",args[0].c_str());
 
134
      delete mytes;
 
135
      mytes=NULL;
 
136
      return 100;
 
137
    }
 
138
    dimx=mytes->dimx;
 
139
    dimy=mytes->dimy;
 
140
    dimz=mytes->dimz;
 
141
    voxsize[0]=mytes->voxsize[0];
 
142
    voxsize[1]=mytes->voxsize[1];
 
143
    voxsize[2]=mytes->voxsize[2];
 
144
    mode=vb_4d;
 
145
    break;
 
146
  case 3:
 
147
    mycube = new CrunchCube;
 
148
    mycube->ReadFile(args[0]);
 
149
    if (!mycube->data_valid) {
 
150
      printf("*** resample: invalid 3D file: %s\n",args[0].c_str());
 
151
      delete mycube;
 
152
      mycube=NULL;
 
153
      return 101;
 
154
    }
 
155
    dimx=mycube->dimx;
 
156
    dimy=mycube->dimy;
 
157
    dimz=mycube->dimz;
 
158
    voxsize[0]=mycube->voxsize[0];
 
159
    voxsize[1]=mycube->voxsize[1];
 
160
    voxsize[2]=mycube->voxsize[2];
 
161
    mode=vb_3d;
 
162
    break;
 
163
  default:
 
164
    printf("*** resample: can only resample 3D and 4D data\n");
 
165
    return 102;
 
166
    break;
 
167
  }
 
168
 
 
169
  x1=y1=z1=0;
 
170
  xstep=ystep=zstep=1.0;
 
171
  nx=dimx;
 
172
  ny=dimy;
 
173
  nz=dimz;
 
174
  
 
175
  for (int i=0; i<args.size(); i++) {
 
176
    if (args[i]=="-rr") {
 
177
      if (i+1 < args.size()) {
 
178
        UseZ(args[i+1],0);
 
179
        i++;
 
180
      }
 
181
    }
 
182
    else if (args[i]=="-ra") {
 
183
      if (i+1 < args.size()) {
 
184
        UseCorner(args[i+1]);
 
185
        i++;
 
186
      }
 
187
    }
 
188
    else if (args[i]=="-rz") {
 
189
      if (i+2 < args.size()) {
 
190
        zsize=strtod(args[i+2].c_str(),NULL);
 
191
        UseZ(args[i+1],zsize);
 
192
        i++;
 
193
      }
 
194
    }
 
195
    else if (args[i]=="-rd") {
 
196
      if (i+1 < args.size()) {
 
197
        UseDims(args[i+1]);
 
198
        i++;
 
199
      }
 
200
    }
 
201
    else if (args[i]=="-rc") {
 
202
      if (i+1 < args.size()) {
 
203
        UseTLHC(args[i+1]);
 
204
        i++;
 
205
      }
 
206
    }
 
207
    else if (args[i]=="-r") {
 
208
      if (i+1 < args.size()) {
 
209
        refname=args[i+1];
 
210
      }
 
211
    }
 
212
    else if (args[i]=="-aa" && i+1<args.size()) {
 
213
      SetAlign(args[i],args[i+1]);
 
214
    }
 
215
    else if (args[i]=="-d") {
 
216
      if (i+3 < args.size()) {
 
217
        nx=strtol(args[i+1]);
 
218
        ny=strtol(args[i+2]);
 
219
        nz=strtol(args[i+3]);
 
220
        // FIXME is the following correct?
 
221
        xstep=(double)(dimx)/nx;
 
222
        ystep=(double)(dimy)/ny;
 
223
        zstep=(double)(dimz)/nz;
 
224
        i+=3;
 
225
      }
 
226
    }
 
227
    else if (args[i]=="-nn") {
 
228
      method=vb_nn;
 
229
    }
 
230
    else if (args[i]=="-xx") {
 
231
      if (i+3 < args.size()) {
 
232
        x1=strtod(args[i+1]);
 
233
        xstep=strtod(args[i+2]);
 
234
        nx=strtol(args[i+3]);
 
235
        i+=3;
 
236
      }
 
237
    }
 
238
    else if (args[i]=="-yy") {
 
239
      if (i+3 < args.size()) {
 
240
        y1=strtod(args[i+1]);
 
241
        ystep=strtod(args[i+2]);
 
242
        ny=strtol(args[i+3]);
 
243
        i+=3;
 
244
      }
 
245
    }
 
246
    else if (args[i]=="-zz") {
 
247
      if (i+3 < args.size()) {
 
248
        z1=strtod(args[i+1]);
 
249
        zstep=strtod(args[i+2]);
 
250
        nz=strtol(args[i+3]);
 
251
        i+=3;
 
252
      }
 
253
    }
 
254
    // new flags: padleft, padright, padtop, padbottom, padback, padfront
 
255
    //            croplr
 
256
    //            flipx, flipy, flipz
 
257
    else if (args[i]=="--flipx") {
 
258
      x1=mycube->dimx-1;
 
259
      xstep=-1;
 
260
      nx=mycube->dimx;
 
261
    }
 
262
    else if (args[i]=="--flipy") {
 
263
      y1=mycube->dimy-1;
 
264
      ystep=-1;
 
265
      ny=mycube->dimy;
 
266
    }
 
267
    else if (args[i]=="--flipz") {
 
268
      z1=mycube->dimz-1;
 
269
      zstep=-1;
 
270
      nz=mycube->dimz;
 
271
    }
 
272
  }
 
273
 
 
274
  char tmps[STRINGLEN];
 
275
  sprintf(tmps,"resample_x: start %.6f step %.2f count %d",x1,xstep,nx);
 
276
  newheaders.push_back(tmps);
 
277
  printf("[I] resample: %s\n",tmps);
 
278
 
 
279
  sprintf(tmps,"resample_y: start %.6f step %.2f count %d",y1,ystep,ny);
 
280
  newheaders.push_back(tmps);
 
281
  printf("[I] resample: %s\n",tmps);
 
282
 
 
283
  sprintf(tmps,"resample_z: start %.6f step %.2f count %d",z1,zstep,nz);
 
284
  newheaders.push_back(tmps);
 
285
  printf("[I] resample: %s\n",tmps);
 
286
 
 
287
  return 0;
 
288
}
 
289
 
 
290
Resample::~Resample()
 
291
{
 
292
}
 
293
 
 
294
int
 
295
Resample::SetAlign(const string mode,const string ref)
 
296
{
 
297
  Cube refcube;
 
298
  if (refcube.ReadFile(refname))
 
299
    return 100;
 
300
  if (mode=="-aa" || mode=="-ax") {
 
301
    if (ref=="origin")
 
302
      x1=(float)mycube->origin[0]-(((float)refcube.origin[0]*refcube.voxsize[0])/mycube->voxsize[0]);
 
303
  }
 
304
  if (mode=="-aa" || mode=="-ay") {
 
305
    if (ref=="origin")
 
306
      y1=(float)mycube->origin[1]-(((float)refcube.origin[1]*refcube.voxsize[1])/mycube->voxsize[1]);
 
307
  }
 
308
  if (mode=="-aa" || mode=="-az") {
 
309
    if (ref=="origin")
 
310
      z1=(float)mycube->origin[2]-(((float)refcube.origin[2]*refcube.voxsize[2])/mycube->voxsize[2]);
 
311
  }
 
312
  return 0;
 
313
}
 
314
 
 
315
int
 
316
Resample::UseCorner(const string &refname)
 
317
{
 
318
  Cube refcube;
 
319
  stringstream tmps;
 
320
  if (refcube.ReadFile(refname)) {
 
321
    Tes reftes;
 
322
    if (!(reftes.ReadFile(refname))) {
 
323
      refcube=reftes[0];
 
324
    }
 
325
    else {
 
326
      printf("*** resample: invalid reference file: %s\n",refname.c_str());
 
327
      return 501;       // panic and exit
 
328
    }
 
329
  }
 
330
  // first get the corner coordinates of both images
 
331
  tokenlist ourline,refline;
 
332
  ourline.ParseLine(mycube->GetHeader("AbsoluteCornerPosition:"));
 
333
  refline.ParseLine(refcube.GetHeader("AbsoluteCornerPosition:"));
 
334
  if (ourline.size()!=3) {
 
335
    tmps.str("");
 
336
    tmps << "resample: input image file doesn't have an absolute corner position tag";
 
337
    printErrorMsg(VB_ERROR,tmps.str());
 
338
    return 101;
 
339
  }
 
340
  if (refline.size()!=3) {
 
341
    tmps.str("");
 
342
    tmps << "resample: reference image file doesn't have an absolute corner position tag";
 
343
    printErrorMsg(VB_ERROR,tmps.str());
 
344
    return 101;
 
345
  }
 
346
  double ourpos[3],refpos[3];
 
347
  ourpos[0]=strtod(ourline[0]);
 
348
  ourpos[1]=strtod(ourline[1]);
 
349
  ourpos[2]=strtod(ourline[2]);
 
350
  refpos[0]=strtod(refline[0]);
 
351
  refpos[1]=strtod(refline[1]);
 
352
  refpos[2]=strtod(refline[2]);
 
353
  
 
354
  // find the start point
 
355
  x1=(refpos[0]-ourpos[0])/mycube->voxsize[0];
 
356
  y1=(refpos[1]-ourpos[1])/mycube->voxsize[1];
 
357
  z1=(refpos[2]-ourpos[2])/mycube->voxsize[2];
 
358
  // the sampling interval should provide for 4x resolution in-plane
 
359
  xstep=(refcube.voxsize[0]/4.0)/mycube->voxsize[0];
 
360
  ystep=(refcube.voxsize[1]/4.0)/mycube->voxsize[1];
 
361
  zstep=refcube.voxsize[2]/mycube->voxsize[2];
 
362
  // the number of voxels should be the width of the target image
 
363
  nx=refcube.dimx*4;
 
364
  ny=refcube.dimy*4;
 
365
  nz=refcube.dimz;
 
366
  return 0;
 
367
}
 
368
 
 
369
int
 
370
Resample::UseZ(const string &refname,double zsize)
 
371
{
 
372
  CrunchCube *refcube;
 
373
 
 
374
  refcube = new CrunchCube;
 
375
  refcube->ReadFile(refname);
 
376
  if (!refcube->data_valid) {
 
377
    printf("*** resample: invalid 3D file: %s\n",refname.c_str());
 
378
    delete refcube;
 
379
    return 501;       // panic and exit
 
380
  }
 
381
  double ourstart,ourend,refstart,refend;
 
382
  // old style: startloc/endloc
 
383
  ourstart=strtod(mycube->GetHeader("StartLoc:"));
 
384
  ourend=strtod(mycube->GetHeader("EndLoc:"));
 
385
  refstart=strtod(refcube->GetHeader("StartLoc:"));
 
386
  refend=strtod(refcube->GetHeader("EndLoc:"));
 
387
  // new style: zrange
 
388
  string refzrange=refcube->GetHeader("ZRange:");
 
389
  string ourzrange=mycube->GetHeader("ZRange:");
 
390
  if (refzrange.size()) {
 
391
    tokenlist range(refzrange);
 
392
    refstart=strtod(range[0]);
 
393
    refend=strtod(range[1]);
 
394
  }
 
395
  if (ourzrange.size()) {
 
396
    tokenlist range(ourzrange);
 
397
    ourstart=strtod(range[0]);
 
398
    ourend=strtod(range[1]);
 
399
  }
 
400
 
 
401
  if (zsize < .001)
 
402
    zsize=refcube->voxsize[2];
 
403
  nx=dimx;
 
404
  ny=dimy;
 
405
  z1=(refstart-ourstart)/voxsize[2];
 
406
  nz=(int)((abs(refend-refstart)/zsize)+0.5)+1;
 
407
  zstep=zsize/voxsize[2];
 
408
  delete refcube;
 
409
  return 0;
 
410
}
 
411
 
 
412
int
 
413
Resample::UseDims(const string &refname)
 
414
{
 
415
  CrunchCube *refcube;
 
416
  refcube = new CrunchCube;
 
417
  refcube->ReadFile(refname);
 
418
  if (!refcube->data_valid) {
 
419
    printf("*** resample: invalid 3D file: %s\n",args[1].c_str());
 
420
    delete refcube;
 
421
    return 501;       // panic and exit
 
422
  }
 
423
  nx=refcube->dimx;
 
424
  ny=refcube->dimy;
 
425
  nz=refcube->dimz;
 
426
  // FIXME is the following correct?
 
427
  xstep=(double)(dimx)/nx;
 
428
  ystep=(double)(dimy)/ny;
 
429
  zstep=(double)(dimz)/nz;
 
430
  delete refcube;
 
431
  return 0;
 
432
}
 
433
 
 
434
int
 
435
Resample::UseTLHC(const string &refname)
 
436
{
 
437
  CrunchCube *refcube;
 
438
 
 
439
  refcube = new CrunchCube;
 
440
  refcube->ReadFile(refname);
 
441
  if (!refcube->data_valid) {
 
442
    printf("*** resample: invalid 3D file: %s\n",refname.c_str());
 
443
    delete refcube;
 
444
    return 501;       // panic and exit
 
445
  }
 
446
  double ourLR,refLR,ourAP,refAP;
 
447
  ourLR=refLR=ourAP=refAP=0.0;
 
448
 
 
449
  // new style: zrange
 
450
  string reftlhc=refcube->GetHeader("im_tlhc:");
 
451
  string ourtlhc=mycube->GetHeader("im_tlhc:");
 
452
 
 
453
  if (reftlhc.size()) {
 
454
    tokenlist range(reftlhc);
 
455
    refLR=strtod(range[0]);
 
456
    refAP=strtod(range[1]);
 
457
  }
 
458
  if (ourtlhc.size()) {
 
459
    tokenlist range(ourtlhc);
 
460
    ourLR=strtod(range[0]);
 
461
    ourAP=strtod(range[1]);
 
462
  }
 
463
 
 
464
  nx=dimx;
 
465
  ny=dimy;
 
466
  nz=dimz;
 
467
  x1=y1=z1=0;
 
468
  xstep=ystep=zstep=1.0;
 
469
 
 
470
  if (abs(ourLR-refLR)>0.001) {
 
471
    x1=(ourLR-refLR)/mycube->voxsize[0];
 
472
  }
 
473
  if (abs(ourAP-refAP)>0.001) {
 
474
    y1=(refAP-ourAP)/mycube->voxsize[1];
 
475
  }
 
476
  if (x1==0 && y1==0) {
 
477
    printf("resample: no fov adjustment neeeded\n");
 
478
  }
 
479
 
 
480
  delete refcube;
 
481
  return 0;
 
482
}
 
483
 
 
484
int
 
485
Resample::ResampleFile()
 
486
{
 
487
  printf("resample: resampling %s\n",args[0].c_str());
 
488
  if (mode == vb_3d)
 
489
    return Resample3D();
 
490
  else if (mode == vb_4d)
 
491
    return Resample4D();
 
492
  else {
 
493
    printf("*** resample: error\n");
 
494
    return 200;
 
495
  }
 
496
}
 
497
 
 
498
int
 
499
Resample::Resample3D()
 
500
{
 
501
  stringstream tmps;
 
502
 
 
503
  if (ResampleCube())
 
504
    return 101;
 
505
  newcube->SetFileName(args[1]);
 
506
  newcube->SetFileFormat(mycube->GetFileFormat());
 
507
  newcube->header=mycube->header;
 
508
 
 
509
  if (mycube->origin[0]) {
 
510
    newcube->origin[0]=lround((mycube->origin[0]-x1)/xstep);
 
511
    newcube->origin[1]=lround((mycube->origin[1]-y1)/ystep);
 
512
    newcube->origin[2]=lround((mycube->origin[2]-z1)/zstep);
 
513
  }
 
514
 
 
515
  // add new headers
 
516
  for (int i=0; i<(int)newheaders.size(); i++)
 
517
    newcube->AddHeader(newheaders[i]);
 
518
  newcube->AddHeader((string)"resample_date: "+timedate());
 
519
 
 
520
  AdjustCornerAndOrigin(*newcube);
 
521
  if (newcube->WriteFile())
 
522
    printf("*** resample: error writing %s\n",args[1].c_str());
 
523
  else
 
524
    printf("resample: wrote %s.\n",args[1].c_str());
 
525
  delete newcube;
 
526
  return 0;
 
527
}
 
528
 
 
529
int
 
530
Resample::Resample4D()
 
531
{
 
532
  newtes=new Tes();
 
533
  newtes->SetVolume(nx,ny,nz,mytes->dimt,mytes->datatype);
 
534
  newtes->voxsize[0]=fabs(xstep*mytes->voxsize[0]);
 
535
  newtes->voxsize[1]=fabs(ystep*mytes->voxsize[1]);
 
536
  newtes->voxsize[2]=fabs(zstep*mytes->voxsize[2]);
 
537
 
 
538
  if (mytes->origin[0]) {
 
539
    newtes->origin[0]=lround((mytes->origin[0]-x1)/xstep);
 
540
    newtes->origin[1]=lround((mytes->origin[1]-y1)/ystep);
 
541
    newtes->origin[2]=lround((mytes->origin[2]-z1)/zstep);
 
542
  }
 
543
 
 
544
  for (int i=0; i<mytes->dimt; i++) {
 
545
    mycube=new CrunchCube((*mytes)[i]);
 
546
    if (ResampleCube())
 
547
      return 102;
 
548
    newtes->SetCube(i,newcube);
 
549
    delete newcube;
 
550
    delete mycube;
 
551
  }
 
552
 
 
553
  // add new headers
 
554
  for (int i=0; i<(int)newheaders.size(); i++)
 
555
    newtes->AddHeader(newheaders[i]);
 
556
  newtes->AddHeader((string)"resample_date: "+timedate());
 
557
 
 
558
  newtes->SetFileName(args[1]);
 
559
  newtes->SetFileFormat(mytes->GetFileFormat());
 
560
  newtes->header=mytes->header;
 
561
  AdjustCornerAndOrigin(*newtes);
 
562
  if (int err=newtes->WriteFile())
 
563
    printf("*** resample: error %d writing %s\n",err,args[1].c_str());
 
564
  else
 
565
    printf("resample: wrote %s.\n",args[1].c_str());
 
566
  return 0;
 
567
}
 
568
 
 
569
void
 
570
Resample::AdjustCornerAndOrigin(VBImage &im)
 
571
{
 
572
  vector<string> newheader;
 
573
  tokenlist args;
 
574
  for (int i=0; i<(int)im.header.size(); i++) {
 
575
    args.ParseLine(im.header[i]);
 
576
    if (args[0]!="AbsoluteCornerPosition:")
 
577
      newheader.push_back(im.header[i]);
 
578
  }
 
579
  double x,y,z;
 
580
  im.GetCorner(x,y,z);
 
581
  x+=x1*im.voxsize[0];
 
582
  y+=y1*im.voxsize[1];
 
583
  z+=z1*im.voxsize[2];
 
584
  stringstream tmps;
 
585
  tmps << "AbsoluteCornerPosition: " << x << " " << y << " " << z;
 
586
  newheader.push_back(tmps.str());
 
587
  im.header=newheader;
 
588
  // the origin is easier...
 
589
//   im.origin[0]-=lround(x1);
 
590
//   im.origin[1]-=lround(y1);
 
591
//   im.origin[2]-=lround(z1);
 
592
}
 
593
 
 
594
int
 
595
Resample::ResampleCube()
 
596
{
 
597
  if (method==vb_nn)
 
598
    return NNResampleCube();
 
599
  else
 
600
    return SincResampleCube();
 
601
}
 
602
 
 
603
int
 
604
Resample::NNResampleCube()
 
605
{
 
606
  int i,j,k;
 
607
 
 
608
  if (!mycube)
 
609
    return 5;
 
610
 
 
611
  newcube=new CrunchCube();
 
612
  newcube->SetVolume(nx,ny,nz,mycube->datatype);
 
613
  newcube->voxsize[0]=fabs(xstep*mycube->voxsize[0]);
 
614
  newcube->voxsize[1]=fabs(ystep*mycube->voxsize[1]);
 
615
  newcube->voxsize[2]=fabs(zstep*mycube->voxsize[2]);
 
616
 
 
617
  for (k=0; k<nz; k++) {
 
618
    for (i=0; i<nx; i++) {
 
619
      for (j=0; j<ny; j++) {
 
620
        int c1=lround(x1+(xstep * i));
 
621
        int c2=lround(y1+(ystep * j));
 
622
        int c3=lround(z1+(zstep * k));
 
623
        newcube->SetValue(i,j,k,mycube->GetValue(c1,c2,c3));
 
624
      }
 
625
    }
 
626
  }
 
627
 
 
628
  return 0;   // no error
 
629
}
 
630
 
 
631
int
 
632
Resample::SincResampleCube()
 
633
{
 
634
  int i,j,k;
 
635
 
 
636
  if (!mycube)
 
637
    return 5;
 
638
 
 
639
  newcube=new CrunchCube();
 
640
  newcube->SetVolume(nx,ny,nz,mycube->datatype);
 
641
  newcube->voxsize[0]=fabs(xstep*mycube->voxsize[0]);
 
642
  newcube->voxsize[1]=fabs(ystep*mycube->voxsize[1]);
 
643
  newcube->voxsize[2]=fabs(zstep*mycube->voxsize[2]);
 
644
 
 
645
  RowVector c1(1),c2(1),c3(1),out(1);
 
646
 
 
647
  for (k=0; k<nz; k++) {
 
648
    for (i=0; i<nx; i++) {
 
649
      for (j=0; j<ny; j++) {
 
650
        c1(0)=x1+(xstep * i)+1;   // +1 because the algorithm 1-indexes
 
651
        c2(0)=y1+(ystep * j)+1;
 
652
        c3(0)=z1+(zstep * k)+1;
 
653
        switch (mycube->datatype) {
 
654
        case vb_byte:
 
655
          resample_sinc(1,(unsigned char *)mycube->data,out,c1,c2,c3,
 
656
                        mycube->dimx,mycube->dimy,mycube->dimz,5,
 
657
                        0.0,1.0);
 
658
          break;
 
659
        case vb_short:
 
660
          resample_sinc(1,(short *)mycube->data,out,c1,c2,c3,
 
661
                        mycube->dimx,mycube->dimy,mycube->dimz,5,
 
662
                        0.0,1.0);
 
663
          break;
 
664
        case vb_long:
 
665
          resample_sinc(1,(int *)mycube->data,out,c1,c2,c3,
 
666
                        mycube->dimx,mycube->dimy,mycube->dimz,5,
 
667
                        0.0,1.0);
 
668
          break;
 
669
        case vb_float:
 
670
          resample_sinc(1,(float *)mycube->data,out,c1,c2,c3,
 
671
                        mycube->dimx,mycube->dimy,mycube->dimz,5,
 
672
                        0.0,1.0);
 
673
          break;
 
674
        case vb_double:
 
675
          resample_sinc(1,(double *)mycube->data,out,c1,c2,c3,
 
676
                        mycube->dimx,mycube->dimy,mycube->dimz,5,
 
677
                        0.0,1.0);
 
678
          break;
 
679
        }
 
680
        newcube->SetValue(i,j,k,out(0));
 
681
      }
 
682
    }
 
683
  }
 
684
 
 
685
  return 0;   // no error
 
686
}
 
687
 
 
688
void
 
689
resample_help()
 
690
{
 
691
  printf("\nVoxBo resample (v%s)\n",vbversion.c_str());
 
692
  printf("usage: resample <in> <out> [flags]\n");
 
693
  printf("flags:\n");
 
694
  printf("   -rr <ref>      ref=image in desired space\n");
 
695
  printf("   -rc <ref>      ref=image with desired field of view\n");
 
696
  printf("   -rz <ref> [z]  same but specify size in z\n");
 
697
  printf("   -rd <ref>      ref=image of desired dimensions\n");
 
698
  printf("   -ra <ref>      aligns to ref on all dims, using corner position\n");
 
699
  printf("   -d  <x y z>    desired dimensions\n");
 
700
  printf("   -nn            use nearest neighbor interpolation (default is sinc)\n");
 
701
  printf("   -xx <a b c>    a=start voxel (from 0); b=interval in voxels; c=# of voxels\n");
 
702
  printf("   -yy <a b c>    same for y dimension\n");
 
703
  printf("   -zz <a b c>    same for z dimension\n");
 
704
  printf("\n");
 
705
  printf("   -r     generic ref cube\n");
 
706
  printf("   -aa, -ax, -ay, -az  [origin/corner]\n");
 
707
  printf("   -sa, -sx, -sy, -sz  [n]\n");
 
708
  printf("\n");
 
709
  printf("Some examples for an axial image of dimensions 64x64x21:\n");
 
710
  printf("  invert the order of the slices: resample in.cub out.cub -zz 20 -1 21\n");
 
711
  printf("  swap left-right: resample in.cub out.cub -xx 63 -1 64\n");
 
712
  printf("  interpolate in-plane: resample in.cub out.cub -xx 0 .25 256 -yy 0 .25 256\n");
 
713
  printf("  crop to 40x64x21: resample in.cub out.cub -xx 12 1 40\n");
 
714
  printf("\n");
 
715
}