vtkMNITransformReader.cxx 25.6 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkMNITransformReader.cxx

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/
/*=========================================================================

Copyright (c) 2006 Atamai, Inc.

Use, modification and redistribution of the software, in source or
binary forms, are permitted provided that the following terms and
conditions are met:

1) Redistribution of the source code, in verbatim or modified
   form, must retain the above copyright notice, this license,
   the following disclaimer, and any notices that refer to this
   license and/or the following disclaimer.

2) Redistribution in binary form must include the above copyright
   notice, a copy of this license and the following disclaimer
   in the documentation or with other materials provided with the
   distribution.

3) Modified copies of the source code must be clearly marked as such,
   and must not be misrepresented as verbatim copies of the source code.

THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES PROVIDE THE SOFTWARE "AS IS"
WITHOUT EXPRESSED OR IMPLIED WARRANTY INCLUDING, BUT NOT LIMITED TO,
THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
PURPOSE.  IN NO EVENT SHALL ANY COPYRIGHT HOLDER OR OTHER PARTY WHO MAY
MODIFY AND/OR REDISTRIBUTE THE SOFTWARE UNDER THE TERMS OF THIS LICENSE
BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, LOSS OF DATA OR DATA BECOMING INACCURATE
OR LOSS OF PROFIT OR BUSINESS INTERRUPTION) ARISING IN ANY WAY OUT OF
THE USE OR INABILITY TO USE THE SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGES.

=========================================================================*/

#include "vtkMNITransformReader.h"

#include "vtkObjectFactory.h"

#include "vtkImageData.h"
#include "vtkMINCImageReader.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkCollection.h"
#include "vtkTransform.h"
#include "vtkGeneralTransform.h"
#include "vtkThinPlateSplineTransform.h"
#include "vtkGridTransform.h"
#include "vtkDoubleArray.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkPoints.h"

66
#include <cctype>
67

68 69 70
#include <string>
#include <vector>
#include <vtksys/SystemTools.hxx>
71 72 73 74 75 76 77

//--------------------------------------------------------------------------
vtkStandardNewMacro(vtkMNITransformReader);

//-------------------------------------------------------------------------
vtkMNITransformReader::vtkMNITransformReader()
{
78 79
  this->FileName = nullptr;
  this->Transform = nullptr;
80 81
  this->Transforms = vtkCollection::New();
  this->LineNumber = 0;
82
  this->Comments = nullptr;
83 84 85 86 87 88
}

//-------------------------------------------------------------------------
vtkMNITransformReader::~vtkMNITransformReader()
{
  if (this->Transforms)
89
  {
90
    this->Transforms->Delete();
91
  }
92
  if (this->Transform)
93
  {
94
    this->Transform->Delete();
95
  }
96 97
  delete [] this->FileName;
  delete [] this->Comments;
98 99 100 101 102 103 104 105 106 107 108
}

//-------------------------------------------------------------------------
void vtkMNITransformReader::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os,indent);

  os << indent << "FileName: "
     << (this->FileName ? this->FileName : "none") << "\n";
  os << indent << "Transform: " << this->Transform << "\n";
  if (this->Transform)
109
  {
110
    this->Transform->PrintSelf(os, indent.GetNextIndent());
111
  }
112 113 114 115 116 117 118 119 120 121 122
  os << indent << "NumberOfTransforms: "
     << this->Transforms->GetNumberOfItems() << "\n";
  os << indent << "Comments: "
     << (this->Comments ? this->Comments : "none") << "\n";
}

//-------------------------------------------------------------------------
int vtkMNITransformReader::CanReadFile(const char* fname)
{
  // First make sure the file exists.  This prevents an empty file
  // from being created on older compilers.
123 124
  vtksys::SystemTools::Stat_t fs;
  if (vtksys::SystemTools::Stat(fname, &fs) != 0)
125
  {
126
    return 0;
127
  }
128 129 130 131 132 133 134

  // Try to read the first line of the file.
  int status = 0;

  ifstream infile(fname);

  if (infile.good())
135
  {
136 137 138 139
    status = 1;
    char linetext[256];
    infile.getline(linetext, 256);
    if (strncmp(linetext, "MNI Transform File", 18) != 0)
140
    {
141
      status = 0;
142
    }
143 144

    infile.close();
145
  }
146 147 148 149 150 151 152 153 154 155 156 157 158 159

  return status;
}

//-------------------------------------------------------------------------
// Internal function to read in a line up to 256 characters and then
// skip to the next line in the file.
int vtkMNITransformReader::ReadLine(
  istream &infile, char result[256])
{
  this->LineNumber++;

  infile.getline(result,256);
  if (infile.fail())
160
  {
161
    if (infile.eof())
162
    {
163
      return 0;
164
    }
165
    if (infile.gcount() == 255)
166
    {
167 168 169 170 171 172
      // Read 256 chars; ignoring the rest of the line.
      infile.clear();
      infile.ignore(VTK_INT_MAX, '\n');
      vtkWarningMacro("Overlength line (limit is 255) in "
                      << this->FileName << ":" << this->LineNumber);
    }
173
  }
174 175 176 177 178 179 180 181 182 183 184

  return 1;
}

//-------------------------------------------------------------------------
// Skip all blank lines or comment lines and return the first useful line
int vtkMNITransformReader::ReadLineAfterComments(
  istream &infile, char result[256])
{
  // Skip over any comment lines or blank lines.
  // Comment lines start with '%'
185
  std::string comments;
186
  do
187
  {
188 189 190
    this->ReadLine(infile, result);
    const char *cp = result;
    while (*cp && isspace(*cp))
191
    {
192
      cp++;
193
    }
194
    if (result[0] == '%')
195
    {
196
      if (comments.length() > 0)
197
      {
198 199
        comments.append("\n");
      }
200 201
      comments.append(result);
    }
202
    else if (*cp != '\0')
203
    {
204
      delete [] this->Comments;
205
      this->Comments = new char[comments.length() + 1];
206 207
      strncpy(this->Comments, comments.c_str(), comments.length());
      this->Comments[comments.length()] = '\0';
208 209
      return 1;
    }
210
  }
211 212 213 214 215 216 217 218 219 220 221 222 223
  while (infile.good());

  return 0;
}

//-------------------------------------------------------------------------
// Skip all whitespace, reading additional lines if necessary
int vtkMNITransformReader::SkipWhitespace(
  istream &infile, char linetext[256], char **cpp)
{
  char *cp = *cpp;

  while (infile.good())
224
  {
225 226
    // Skip leading whitespace
    while (isspace(*cp))
227
    {
228
      cp++;
229
    }
230 231

    if (*cp != '\0')
232
    {
233 234
      *cpp = cp;
      return 1;
235
    }
236 237 238

    this->ReadLine(infile, linetext);
    cp = linetext;
239
  }
240 241 242 243 244 245 246 247 248 249 250 251 252 253 254

  return 0;
}

//-------------------------------------------------------------------------
// Read the left hand side of a statement, including the equals sign
// and any whitespace following the equals.
int vtkMNITransformReader::ParseLeftHandSide(
  istream &infile, char linetext[256], char **cpp, char identifier[256])
{
  int i = 0;
  char *cp = *cpp;

  // Read alphanumeric plus underscore
  if (!isdigit(*cp))
255
  {
256
    while ((isalnum(*cp) || *cp == '_') && i < 255)
257
    {
258 259
      identifier[i++] = *cp++;
    }
260
  }
261 262 263 264
  identifier[i] = '\0';

  // Skip trailing whitespace
  while (isspace(*cp))
265
  {
266
    cp++;
267
  }
268 269 270 271

  // Check for equals
  this->SkipWhitespace(infile, linetext, &cp);
  if (*cp != '=')
272
  {
273 274 275
    vtkErrorMacro("Missing \'=\' " << this->FileName
                  << ":" << this->LineNumber);
    return 0;
276
  }
277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301
  cp++;

  // Skip ahead to the value part of the statement
  this->SkipWhitespace(infile, linetext, &cp);

  *cpp = cp;

  return 1;
}

//-------------------------------------------------------------------------
// Read a string value.  The terminating semicolon will be read, but
// won't be included in the output string.  Neither will any
// whitespace occurring before the semicolon. The string may not be
// split across multiple lines.
int vtkMNITransformReader::ParseStringValue(
  istream &infile, char linetext[256], char **cpp, char data[256])
{
  int i = 0;
  char *cp = *cpp;

  this->SkipWhitespace(infile, linetext, &cp);

  // Read until end of the line or semicolon
  while (*cp && *cp != ';' && i < 255)
302
  {
303
    data[i++] = *cp++;
304
  }
305 306 307

  // Remove trailing whitespace
  while (i > 0 && isspace(data[i-1]))
308
  {
309
    i--;
310
  }
311 312 313 314 315

  data[i] = '\0';

  this->SkipWhitespace(infile, linetext, &cp);
  if (*cp != ';')
316
  {
317 318 319
    vtkErrorMacro("Missing semicolon " << this->FileName
                  << ":" << this->LineNumber);
    return 0;
320
  }
321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337
  cp++;

  *cpp = cp;

  return 1;
}

//-------------------------------------------------------------------------
// Read floating-point values into a vtkDoubleArray until a semicolon
// is reached.  The semicolon is also read.
int vtkMNITransformReader::ParseFloatValues(
  istream &infile, char linetext[256], char **cpp, vtkDoubleArray *array)
{
  char *cp = *cpp;

  this->SkipWhitespace(infile, linetext, &cp);
  while (infile.good() && *cp != ';')
338
  {
339 340 341
    char *tmp = cp;
    double val = strtod(cp, &cp);
    if (cp == tmp)
342
    {
343 344 345
      vtkErrorMacro("Syntax error " << this->FileName
                    << ":" << this->LineNumber);
      return 0;
346
    }
347 348
    array->InsertNextValue(val);
    this->SkipWhitespace(infile, linetext, &cp);
349
  }
350 351

  if (*cp != ';')
352
  {
353 354 355
    vtkErrorMacro("Missing semicolon " << this->FileName
                  << ":" << this->LineNumber);
    return 0;
356
  }
357 358 359 360 361 362 363 364 365 366 367 368 369 370
  cp++;

  *cpp = cp;

  return 1;
}

//-------------------------------------------------------------------------
int vtkMNITransformReader::ParseInvertFlagValue(
  istream &infile, char linetext[256], char **cpp, int *invertFlag)
{
  char data[256];

  if (!this->ParseStringValue(infile, linetext, cpp, data))
371
  {
372
    return 0;
373
  }
374
  if (strcmp(data, "False") == 0)
375
  {
376
    *invertFlag = 0;
377
  }
378
  else if (strcmp(data, "True") == 0)
379
  {
380
    *invertFlag = 1;
381
  }
382
  else
383
  {
384 385 386
    vtkErrorMacro("Invert_Flag must be \'True\' or \'False\' "
                  << this->FileName << ":" << this->LineNumber);
    return 0;
387
  }
388 389 390 391 392 393 394 395 396 397 398 399 400

  return 1;
}


//-------------------------------------------------------------------------
int vtkMNITransformReader::ReadLinearTransform(
  istream &infile, char linetext[256], char **cpp)
{
  // Read the first variable
  this->SkipWhitespace(infile, linetext, cpp);
  char identifier[256];
  if (!this->ParseLeftHandSide(infile, linetext, cpp, identifier))
401
  {
402
    return 0;
403
  }
404 405 406 407

  // Check for Invert_Flag
  int invertFlag = 0;
  if (strcmp(identifier, "Invert_Flag") == 0)
408
  {
409
    if (!this->ParseInvertFlagValue(infile, linetext, cpp, &invertFlag))
410
    {
411
      return 0;
412
    }
413 414 415

    this->SkipWhitespace(infile, linetext, cpp);
    if (!this->ParseLeftHandSide(infile, linetext, cpp, identifier))
416
    {
417 418
      return 0;
    }
419
  }
420 421 422

  // Check for Linear_Transform
  if (strcmp(identifier, "Linear_Transform") != 0)
423
  {
424 425 426
    vtkErrorMacro("Expected \'Linear_Transform\' in "
                  << this->FileName << ":" << this->LineNumber);
    return 0;
427
  }
428 429 430 431

  // Read twelve array elements from the file
  vtkDoubleArray *array = vtkDoubleArray::New();
  if (!this->ParseFloatValues(infile, linetext, cpp, array))
432
  {
433
    return 0;
434
  }
435 436

  if (array->GetNumberOfTuples() != 12)
437
  {
438 439 440 441
    vtkErrorMacro("Linear transform must have exactly 12 elements "
                  << this->FileName << ":" << this->LineNumber);
    array->Delete();
    return 0;
442
  }
443 444 445 446 447 448 449 450 451 452 453 454

  // Fill in the last row of the 4x4 matrix
  array->InsertNextValue(0.0);
  array->InsertNextValue(0.0);
  array->InsertNextValue(0.0);
  array->InsertNextValue(1.0);

  // Create the transform
  vtkTransform *transform = vtkTransform::New();
  transform->Concatenate(array->GetPointer(0));
  array->Delete();
  if (invertFlag)
455
  {
456
    transform->Inverse();
457
  }
458 459 460 461 462 463 464 465 466 467 468 469 470 471 472

  this->Transforms->AddItem(transform);
  transform->Delete();

  return 1;
}

//-------------------------------------------------------------------------
int vtkMNITransformReader::ReadThinPlateSplineTransform(
  istream &infile, char linetext[256], char **cpp)
{
  // Read the first variable
  this->SkipWhitespace(infile, linetext, cpp);
  char identifier[256];
  if (!this->ParseLeftHandSide(infile, linetext, cpp, identifier))
473
  {
474
    return 0;
475
  }
476 477 478 479

  // Check for Invert_Flag
  int invertFlag = 0;
  if (strcmp(identifier, "Invert_Flag") == 0)
480
  {
481
    if (!this->ParseInvertFlagValue(infile, linetext, cpp, &invertFlag))
482
    {
483
      return 0;
484
    }
485 486 487

    this->SkipWhitespace(infile, linetext, cpp);
    if (!this->ParseLeftHandSide(infile, linetext, cpp, identifier))
488
    {
489 490
      return 0;
    }
491
  }
492 493 494

  // Number_Dimensions: vtkThinPlateSplineTransform supports 2 and 3
  if (strcmp(identifier, "Number_Dimensions") != 0)
495
  {
496 497 498
    vtkErrorMacro("Expected \'Number_Dimensions\' in "
                  << this->FileName << ":" << this->LineNumber);
    return 0;
499
  }
500 501 502

  char data[256];
  if (!this->ParseStringValue(infile, linetext, cpp, data))
503
  {
504
    return 0;
505
  }
506 507 508

  int numDimensions = data[0] - '0';
  if (data[1] != '\0' || numDimensions < 2 || numDimensions > 3)
509
  {
510 511 512
    vtkErrorMacro("Number_Dimensions must be 2 or 3 in "
                  << this->FileName << ":" << this->LineNumber);
    return 0;
513
  }
514 515 516 517

  // Read the points
  this->SkipWhitespace(infile, linetext, cpp);
  if (!this->ParseLeftHandSide(infile, linetext, cpp, identifier))
518
  {
519
    return 0;
520
  }
521 522

  if (strcmp(identifier, "Points") != 0)
523
  {
524 525 526
    vtkErrorMacro("Expected \'Points\' in "
                  << this->FileName << ":" << this->LineNumber);
    return 0;
527
  }
528 529 530

  vtkDoubleArray *points = vtkDoubleArray::New();
  if (!this->ParseFloatValues(infile, linetext, cpp, points))
531
  {
532 533
    points->Delete();
    return 0;
534
  }
535 536

  if (points->GetNumberOfTuples() % numDimensions != 0)
537
  {
538 539 540 541
    points->Delete();
    vtkErrorMacro("Points list not divisible by Number_Dimensions in "
                  << this->FileName << ":" << this->LineNumber);
    return 0;
542
  }
543 544 545 546

  // Read the displacements
  this->SkipWhitespace(infile, linetext, cpp);
  if (!this->ParseLeftHandSide(infile, linetext, cpp, identifier))
547
  {
548 549
    points->Delete();
    return 0;
550
  }
551 552

  if (strcmp(identifier, "Displacements") != 0)
553
  {
554
    points->Delete();
555 556 557
    vtkErrorMacro("Expected \'Displacements\' in "
                  << this->FileName << ":" << this->LineNumber);
    return 0;
558
  }
559 560 561

  vtkDoubleArray *displacements = vtkDoubleArray::New();
  if (!this->ParseFloatValues(infile, linetext, cpp, displacements))
562
  {
563 564 565
    displacements->Delete();
    points->Delete();
    return 0;
566
  }
567 568 569

  if (displacements->GetNumberOfTuples() !=
      points->GetNumberOfTuples() + numDimensions*(numDimensions + 1))
570
  {
571 572
    displacements->Delete();
    points->Delete();
Kunda's avatar
Kunda committed
573
    vtkErrorMacro("Incorrect number of Displacements in "
574 575
                  << this->FileName << ":" << this->LineNumber);
    return 0;
576
  }
577 578 579 580 581 582 583 584 585 586 587 588 589 590

  // The vtkThinPlateSplineTransform expects two sets of points,
  // not a set of points and a set of displacements.  We need
  // to perform a thin-plate spline transform to get the points
  // that we need.

  int numPoints = points->GetNumberOfTuples()/numDimensions;
  int i = 0;
  int j = 0;

  // Convert points and displacements to 3D
  double (*q)[3] = new double[numPoints][3];
  double (*W)[3] = new double[numPoints][3];
  for (i = 0; i < numPoints; i++)
591
  {
592 593 594 595 596
    double *point = q[i];
    double *displacement = W[i];
    point[0] = point[1] = point[2] = 0.0;
    displacement[0] = displacement[1] = displacement[2] = 0.0;
    for (j = 0; j < numDimensions; j++)
597
    {
598 599 600
      point[j] = points->GetValue(i*numDimensions + j);
      displacement[j] = displacements->GetValue(i*numDimensions + j);
    }
601
  }
602 603 604 605 606

  // Get the translation from the TPS matrix
  double C[3];
  C[0] = C[1] = C[2] = 0.0;
  for (j = 0; j < numDimensions; j++)
607
  {
608
    C[j] = displacements->GetValue(numPoints*numDimensions + j);
609
  }
610 611 612 613 614 615

  // Get the square matrix portion of the TPS matrix
  double A[3][3];
  A[0][1] = A[0][2] = A[1][0] = A[1][2] = A[2][0] = A[2][1] = 0.0;
  A[0][0] = A[1][1] = A[2][2] = 1.0;
  for (i = 0; i < numDimensions; i++)
616
  {
617
    for (j = 0; j < numDimensions; j++)
618
    {
619 620 621
      A[i][j] = displacements->GetValue(
        (numPoints + 1 + i)*numDimensions + j);
    }
622
  }
623

624 625 626
  displacements->Delete();
  points->Delete();

627 628 629 630
  // Create the source and target point lists
  vtkPoints *source = vtkPoints::New();
  vtkPoints *target = vtkPoints::New();
  for (i = 0; i < numPoints; i++)
631
  {
632 633 634 635 636 637
    double *p = q[i];
    double x = 0.0;
    double y = 0.0;
    double z = 0.0;

    for (j = 0; j < numPoints; j++)
638
    {
639 640 641 642
      double dx = p[0] - q[j][0];
      double dy = p[1] - q[j][1];
      double dz = p[2] - q[j][2];
      double r = sqrt(dx*dx + dy*dy + dz*dz);
643
      double U = ((numDimensions == 2 && r != 0) ? r*r*log(r) : r);
644 645 646 647 648 649 650 651 652 653 654
      x += U*W[j][0];
      y += U*W[j][1];
      z += U*W[j][2];
    }

    x += C[0] + p[0]*A[0][0] + p[1]*A[1][0] + p[2]*A[2][0];
    y += C[1] + p[0]*A[0][1] + p[1]*A[1][1] + p[2]*A[2][1];
    z += C[2] + p[0]*A[0][2] + p[1]*A[1][2] + p[2]*A[2][2];

    source->InsertNextPoint(p);
    target->InsertNextPoint(x, y, z);
655
  }
656 657 658 659 660 661 662 663 664

  delete [] q;
  delete [] W;

  // Create the thin plate spline transform
  vtkThinPlateSplineTransform *transform = vtkThinPlateSplineTransform::New();
  transform->SetSourceLandmarks(source);
  transform->SetTargetLandmarks(target);
  if (numDimensions == 2)
665
  {
666
    transform->SetBasisToR2LogR();
667
  }
668
  else
669
  {
670
    transform->SetBasisToR();
671
  }
672 673

  if (invertFlag)
674
  {
675
    transform->Inverse();
676
  }
677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694

  source->Delete();
  target->Delete();

  this->Transforms->AddItem(transform);
  transform->Delete();

  return 1;
}

//-------------------------------------------------------------------------
int vtkMNITransformReader::ReadGridTransform(
  istream &infile, char linetext[256], char **cpp)
{
  // Read the first variable
  this->SkipWhitespace(infile, linetext, cpp);
  char identifier[256];
  if (!this->ParseLeftHandSide(infile, linetext, cpp, identifier))
695
  {
696
    return 0;
697
  }
698 699 700 701

  // Check for Invert_Flag
  int invertFlag = 0;
  if (strcmp(identifier, "Invert_Flag") == 0)
702
  {
703
    if (!this->ParseInvertFlagValue(infile, linetext, cpp, &invertFlag))
704
    {
705
      return 0;
706
    }
707 708 709

    this->SkipWhitespace(infile, linetext, cpp);
    if (!this->ParseLeftHandSide(infile, linetext, cpp, identifier))
710
    {
711 712
      return 0;
    }
713
  }
714 715 716

  // Displacement_Volume must be a minc file
  if (strcmp(identifier, "Displacement_Volume") != 0)
717
  {
718 719 720
    vtkErrorMacro("Expected \'Displacement_Volume\' in "
                  << this->FileName << ":" << this->LineNumber);
    return 0;
721
  }
722

Bill Lorensen's avatar
Bill Lorensen committed
723
  char filename[VTK_MAXPATH];
724
  if (!this->ParseStringValue(infile, linetext, cpp, filename))
725
  {
726
    return 0;
727
  }
728 729 730 731

  // Create the minc reader
  vtkMINCImageReader *reader = vtkMINCImageReader::New();

732 733
  std::vector<std::string> xfmpath;
  std::vector<std::string> mincpath;
734 735 736 737 738

  vtksys::SystemTools::SplitPath(this->FileName, xfmpath);
  vtksys::SystemTools::SplitPath(filename, mincpath);

  // Join minc filename to this->FileName if filename is relative
739
  if (mincpath[0].empty())
740
  {
741 742 743
    xfmpath.pop_back();
    xfmpath.insert(xfmpath.end(), mincpath.begin()+1, mincpath.end());
    reader->SetFileName(vtksys::SystemTools::JoinPath(xfmpath).c_str());
744
  }
745
  else
746
  {
747
    reader->SetFileName(filename);
748
  }
749 750

  // Read the minc file now, rather than later
751
  reader->Update();
752 753 754

  // Create the transform
  vtkGridTransform *transform = vtkGridTransform::New();
755
  transform->SetDisplacementGridConnection(reader->GetOutputPort());
756 757 758 759 760 761
  transform->SetDisplacementShift(reader->GetRescaleIntercept());
  transform->SetDisplacementScale(reader->GetRescaleSlope());
  transform->SetInverseTolerance(0.05);
  transform->SetInterpolationModeToCubic();

  if (invertFlag)
762
  {
763
    transform->Inverse();
764
  }
765 766 767 768 769 770 771 772 773 774 775 776 777 778

  reader->Delete();

  this->Transforms->AddItem(transform);
  transform->Delete();

  return 1;
}

//-------------------------------------------------------------------------
int vtkMNITransformReader::ReadNextTransform(istream &infile, char linetext[256])
{
  // Check for errors
  if (infile.eof())
779
  {
780
    return 1;
781
  }
782
  else if (infile.fail())
783
  {
784 785
    vtkErrorMacro("IO error while reading " << this->FileName);
    return 0;
786
  }
787 788 789 790 791 792 793 794

  // Parse the file
  char *cp = linetext;

  // Check for Transform_Type
  char identifier[256];
  this->SkipWhitespace(infile, linetext, &cp);
  if (!this->ParseLeftHandSide(infile, linetext, &cp, identifier))
795
  {
796
    return 0;
797
  }
798 799

  if (strcmp(identifier, "Transform_Type") != 0)
800
  {
801 802 803
    vtkErrorMacro("Expected Transform_Type in "
                  << this->FileName << ":" << this->LineNumber);
    return 0;
804
  }
805 806 807 808

  // Read the transform type
  char transformType[256];
  if (!this->ParseStringValue(infile, linetext, &cp, transformType))
809
  {
810
    return 0;
811
  }
812 813 814

  // Check the transform type
  if (strcmp(transformType, "Linear") == 0)
815
  {
816
    return this->ReadLinearTransform(infile, linetext, &cp);
817
  }
818
  else if (strcmp(transformType, "Thin_Plate_Spline_Transform") == 0)
819
  {
820
    return this->ReadThinPlateSplineTransform(infile, linetext, &cp);
821
  }
822
  else if (strcmp(transformType, "Grid_Transform") == 0)
823
  {
824
    return this->ReadGridTransform(infile, linetext, &cp);
825
  }
826 827 828 829 830 831 832 833 834 835

  vtkErrorMacro("Unrecognized type " << transformType << " in "
                << this->FileName << ":" << this->LineNumber);
  return 0;
}

//-------------------------------------------------------------------------
int vtkMNITransformReader::ReadFile()
{
  this->Transforms->RemoveAllItems();
836
  this->SetTransform(nullptr);
837 838 839

  // Check that the file name has been set.
  if (!this->FileName)
840
  {
841 842
    vtkErrorMacro("ReadFile: No file name has been set");
    return 0;
843
  }
844 845

  // Make sure that the file exists.
846 847
  vtksys::SystemTools::Stat_t fs;
  if (vtksys::SystemTools::Stat(this->FileName, &fs) != 0)
848
  {
849 850
    vtkErrorMacro("ReadFile: Can't open file " << this->FileName);
    return 0;
851
  }
852 853 854 855 856

  // Make sure that the file is readable.
  ifstream infile(this->FileName);

  if (infile.fail())
857
  {
858 859
    vtkErrorMacro("ReadFile: Can't read the file " << this->FileName);
    return 0;
860
  }
861 862 863 864 865 866

  // Read the first line
  char linetext[256];
  this->LineNumber = 0;
  this->ReadLine(infile, linetext);
  if (strncmp(linetext, "MNI Transform File", 18) != 0)
867
  {
868 869 870
    vtkErrorMacro("ReadFile: File is not a MNI xfm file: " << this->FileName);
    infile.close();
    return 0;
871
  }
872 873 874 875 876 877

  // Read the comments
  this->ReadLineAfterComments(infile, linetext);

  // Read the transforms
  while (infile.good())
878
  {
879
    if (this->ReadNextTransform(infile, linetext) == 0)
880
    {
881 882 883 884
      this->Transforms->RemoveAllItems();
      infile.close();
      return 0;
    }
885 886
    this->ReadLine(infile, linetext);
  }
887 888 889 890 891 892 893

  // Close the file
  infile.close();

  // Create the output transform.
  int n = this->Transforms->GetNumberOfItems();
  if (n == 1)
894
  {
895 896
    this->SetTransform(
      (vtkAbstractTransform *)this->Transforms->GetItemAsObject(0));
897
  }
898
  else
899
  {
900 901 902 903
    // Determine whether the full transform is linear
    int linear = 1;
    int i = 0;
    for (i = 0; i < n; i++)
904
    {
905
      if (!this->Transforms->GetItemAsObject(i)->IsA("vtkLinearTransform"))
906
      {
907 908 909
        linear = 0;
        break;
      }
910
    }
911 912 913 914

    // If linear, use vtkTransform to concatenate,
    // else use vtkGeneralTransform.
    if (linear)
915
    {
916 917 918
      vtkTransform *transform = vtkTransform::New();
      transform->PostMultiply();
      for (i = 0; i < n; i++)
919
      {
920 921 922
        vtkLinearTransform *linearTransform =
          (vtkLinearTransform *)this->Transforms->GetItemAsObject(i);
        transform->Concatenate(linearTransform->GetMatrix());
923
      }
924 925
      this->SetTransform(transform);
      transform->Delete();
926
    }
927
    else
928
    {
929 930 931
      vtkGeneralTransform *transform = vtkGeneralTransform::New();
      transform->PostMultiply();
      for (i = 0; i < n; i++)
932
      {
933 934 935
        vtkAbstractTransform *abstractTransform =
          (vtkAbstractTransform *)this->Transforms->GetItemAsObject(i);
        if (abstractTransform->IsA("vtkLinearTransform"))
936
        {
937 938
          transform->Concatenate(
            ((vtkLinearTransform *)abstractTransform)->GetMatrix());
939
        }
940
        else
941
        {
942 943
          transform->Concatenate(abstractTransform);
        }
944
      }
945 946 947
      this->SetTransform(transform);
      transform->Delete();
    }
948
  }
949 950 951 952 953 954 955 956 957 958

  return 1;
}

//-------------------------------------------------------------------------
int vtkMNITransformReader::ProcessRequest(vtkInformation *request,
                                    vtkInformationVector **inputVector,
                                    vtkInformationVector *outputVector)
{
  if(request->Has(vtkDemandDrivenPipeline::REQUEST_DATA()))
959
  {
960
    return this->ReadFile();
961
  }
962 963 964 965 966 967 968 969

  return this->Superclass::ProcessRequest(request, inputVector, outputVector);
}

//-------------------------------------------------------------------------
void vtkMNITransformReader::SetTransform(vtkAbstractTransform *transform)
{
  if (this->Transform != transform)
970
  {
971
    if (this->Transform)
972
    {
973
      this->Transform->Delete();
974
    }
975
    if (transform)
976
    {
977 978
      transform->Register(this);
    }
979 980
    this->Transform = transform;
  }
981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004
}

//-------------------------------------------------------------------------
vtkAbstractTransform *vtkMNITransformReader::GetTransform()
{
  this->Update();

  return this->Transform;
}

//-------------------------------------------------------------------------
int vtkMNITransformReader::GetNumberOfTransforms()
{
  this->Update();

  return this->Transforms->GetNumberOfItems();
}

//-------------------------------------------------------------------------
vtkAbstractTransform *vtkMNITransformReader::GetNthTransform(int i)
{
  this->Update();

  if (i < 0 || i >= this->Transforms->GetNumberOfItems())
1005
  {
1006
    return nullptr;
1007
  }
1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018

  return (vtkAbstractTransform *)this->Transforms->GetItemAsObject(i);
}

//-------------------------------------------------------------------------
const char *vtkMNITransformReader::GetComments()
{
  this->Update();

  return this->Comments;
}