vtkMINCImageWriter.cxx 63.2 KB
Newer Older
1 2 3 4 5
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkMINCImageWriter.cxx

6 7 8 9 10 11 12 13 14 15 16
  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.

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

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 66 67 68 69 70 71
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 "vtkMINCImageWriter.h"

#include "vtkObjectFactory.h"

#include "vtkImageData.h"
#include "vtkStringArray.h"
#include "vtkCharArray.h"
#include "vtkSignedCharArray.h"
#include "vtkUnsignedCharArray.h"
#include "vtkShortArray.h"
#include "vtkIntArray.h"
#include "vtkFloatArray.h"
#include "vtkDoubleArray.h"
#include "vtkIdTypeArray.h"
#include "vtkMatrix4x4.h"
#include "vtkSmartPointer.h"
#include "vtkMath.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkStreamingDemandDrivenPipeline.h"

#include "vtkType.h"

72
#include "vtkMINCImageAttributes.h"
73
#include "vtkMINC.h"
74
#include "vtk_netcdf.h"
75

76 77 78 79 80 81 82
#ifdef _WIN32
#include "vtkWindows.h"
#else
#include <sys/types.h>
#include <unistd.h>
#endif

83
#include <cstdlib>
84
#include <cfloat>
85
#include <ctime>
86 87 88
#include <string>
#include <vector>
#include <map>
89 90 91 92 93 94

#define VTK_MINC_MAX_DIMS 8

//--------------------------------------------------------------------------
vtkStandardNewMacro(vtkMINCImageWriter);

95
vtkCxxSetObjectMacro(vtkMINCImageWriter,DirectionCosines,vtkMatrix4x4);
96 97
vtkCxxSetObjectMacro(vtkMINCImageWriter,ImageAttributes,
                     vtkMINCImageAttributes);
98 99 100 101

//-------------------------------------------------------------------------
vtkMINCImageWriter::vtkMINCImageWriter()
{
102
  this->DirectionCosines = nullptr;
103 104
  this->RescaleIntercept = 0.0;
  this->RescaleSlope = 0.0;
105 106
  this->InternalRescaleIntercept = 0.0;
  this->InternalRescaleSlope = 0.0;
107 108 109 110 111

  this->MINCImageType = 0;
  this->MINCImageTypeSigned = 1;
  this->MINCImageMinMaxDims = 0;

112 113 114 115
  this->FileDataType = 0;
  this->FileValidRange[0] = 0.0;
  this->FileValidRange[1] = 1.0;
  this->ComputeValidRangeFromScalarRange = 0;
116

117 118 119 120 121 122 123
  this->DataUpdateExtent[0] = 0;
  this->DataUpdateExtent[1] = 0;
  this->DataUpdateExtent[2] = 0;
  this->DataUpdateExtent[3] = 0;
  this->DataUpdateExtent[4] = 0;
  this->DataUpdateExtent[5] = 0;

124
  this->FileDimensionNames = vtkStringArray::New();
125

126
  this->ImageAttributes = nullptr;
127 128

  this->StrictValidation = 1;
129

130
  this->MismatchedInputs = 0;
131

132
  this->HistoryAddition = nullptr;
133 134 135 136 137
}

//-------------------------------------------------------------------------
vtkMINCImageWriter::~vtkMINCImageWriter()
{
138
  if (this->DirectionCosines)
139
  {
140
    this->DirectionCosines->Delete();
141
    this->DirectionCosines = nullptr;
142
  }
143
  if (this->FileDimensionNames)
144
  {
145
    this->FileDimensionNames->Delete();
146
    this->FileDimensionNames = nullptr;
147
  }
148
  if (this->ImageAttributes)
149
  {
150
    this->ImageAttributes->Delete();
151
    this->ImageAttributes = nullptr;
152
  }
153
  this->SetHistoryAddition(nullptr);
154 155 156 157 158 159 160
}

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

161 162
  os << indent << "DirectionCosines: " << this->DirectionCosines << "\n";
  if (this->DirectionCosines)
163
  {
164
    this->DirectionCosines->PrintSelf(os, indent.GetNextIndent());
165
  }
166 167 168 169
  os << indent << "RescaleSlope: " << this->RescaleSlope << "\n";
  os << indent << "RescaleIntercept: " << this->RescaleIntercept << "\n";
  os << indent << "StrictValidation: " <<
    (this->StrictValidation ? "On\n" : "Off\n");
170
  os << indent << "HistoryAddition: " <<
171
    (this->HistoryAddition ? this->HistoryAddition : "(None)") << "\n";
172 173 174 175 176 177 178 179 180 181 182 183 184
}

//-------------------------------------------------------------------------
void vtkMINCImageWriter::SetFileName(const char *name)
{
  this->Superclass::SetFileName(name);
}

//-------------------------------------------------------------------------
int vtkMINCImageWriter::OpenNetCDFFile(const char *filename, int& ncid)
{
  int status = 0;

185
  if (filename == nullptr)
186
  {
187 188
    vtkErrorMacro("No filename was set");
    return 0;
189
  }
190 191 192

  status = nc_create(filename, 0, &ncid);
  if (status != NC_NOERR)
193
  {
194 195 196
    vtkErrorMacro("Could not open the MINC file:\n"
                  << nc_strerror(status));
    return 0;
197
  }
198 199 200 201 202 203 204 205 206 207

  return 1;
}

//-------------------------------------------------------------------------
int vtkMINCImageWriter::CloseNetCDFFile(int ncid)
{
  int status = 0;
  status = nc_close(ncid);
  if (status != NC_NOERR)
208
  {
209 210 211
    vtkErrorMacro("Could not close the MINC file:\n"
                  << nc_strerror(status));
    return 0;
212
  }
213 214 215 216 217 218 219 220

  return 1;
}

//-------------------------------------------------------------------------
// this is a macro so the vtkErrorMacro will report a useful line number
#define vtkMINCImageWriterFailAndClose(ncid, status) \
{ \
221
  if ((status) != NC_NOERR) \
222
  { \
223 224 225
    vtkErrorMacro("There was an error with the MINC file \"" \
                  << this->GetFileName() << "\":\n" \
                  << nc_strerror(status)); \
226
  } \
227 228 229 230 231 232 233 234
  nc_close(ncid); \
}

//-------------------------------------------------------------------------
// Function for getting VTK dimension index from file name.
int vtkMINCImageWriter::IndexFromDimensionName(const char *dimName)
{
  switch(dimName[0])
235
  {
236 237 238 239 240 241 242
    case 'x':
      return this->Permutation[0];
    case 'y':
      return this->Permutation[1];
    case 'z':
      return this->Permutation[2];
    default:
243
      if (strcmp(dimName, MIvector_dimension) == 0)
244
      {
245
        return -1;
246
      }
247
      break;
248
  }
249 250 251 252 253 254 255 256 257

  // Any unrecognized dimensions are returned as index 3
  return 3;
}

//-------------------------------------------------------------------------
// Compute the default dimension order from the direction cosines,
// and look for flips.
// The way the permutation should be used is as follows:
258 259
// If permutation[0] == 0 then MIxspace is VTK's X dimension.
// If permutation[0] == 2 then MIxspace is VTK's Z dimension.
260 261 262 263 264 265 266 267 268
// If the "flip" is set for a VTK, then that VTK dimension
// and its dircos will have to be flipped before the MINC
// file is written.
// For example, if flip[2] == 1, then the MINC dimension that
// maps to the VTK Z dimension will to be flipped along with
// its dircos.
void vtkMINCImageWriter::ComputePermutationFromOrientation(
  int permutation[3], int flip[3])
{
269
  vtkMatrix4x4 *matrix = this->DirectionCosines;
270
  if (matrix == nullptr)
271
  {
272 273 274 275 276 277 278 279
    permutation[0] = 0;
    permutation[1] = 1;
    permutation[2] = 2;
    flip[0] = 0;
    flip[1] = 0;
    flip[2] = 0;

    return;
280
  }
281 282 283 284 285 286 287 288 289 290 291 292 293

  // There are 6 permutations for 3 dimensions.  In addition,
  // if each of those dimensions can be flipped, then there are
  // 8 (two to the power of three) possible flips.  That would
  // give 48 different possibilities, but since we don't consider
  // any combinations that result in left-handed rotations, the
  // total number of combinations that we test is 24.

  // Convert the matrix into three column vectors
  double vectors[3][4];
  int i = 0;
  int j = 0;
  for (i = 0; i < 3; i++)
294
  {
295 296
    double *v = vectors[i];
    for (j = 0; j < 4; j++)
297
    {
298
      v[j] = 0.0;
299
    }
300 301
    v[i] = 1.0;
    matrix->MultiplyPoint(v, v);
302
  }
303 304 305 306

  // Here's how the algorithm works.  We want to find a matrix
  // composed only of permutations and flips that has the closest
  // possible orientation (in terms of absolute orientation angle)
307
  // to our DirectionCosines.
308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334

  // The orientation angle for any matrix A is given by:
  //
  //   cos(angle/2) = sqrt(1 + trace(A))/2
  //
  // Therefore, the minimum angle occurs when the trace is
  // at its maximum.

  // So our method is to calculate the traces of all the various
  // permutations and flips, and just use the one with the largest
  // trace.

  // First check if the matrix includes an odd number of flips,
  // since if it does, it specifies a left-handed rotation.
  double d = vtkMath::Determinant3x3(vectors[0], vectors[1], vectors[2]);
  int oddPermutation = (d < 0);

  // Calculate all the traces, including any combination of
  // permutations and flips that represent right-handed
  // orientations.
  int imax = 0;
  int jmax = 0;
  int kmax = 0;
  int lmax = 0;
  double maxtrace = -1e30;

  for (i = 0; i < 3; i++)
335
  {
336
    for (j = 0; j < 2; j++)
337
    {
338 339 340 341
      double xval = vectors[i][0];
      double yval = vectors[(i + 1 + j) % 3][1];
      double zval = vectors[(i + 2 - j) % 3][2];
      for (int k = 0; k < 2; k++)
342
      {
343
        for (int l = 0; l < 2; l++)
344
        {
345 346 347 348 349 350 351 352 353 354 355
          // The (1 - 2*k) gives a sign from a boolean.
          // For z, we want to set the sign that will
          // not change the handedness ("^" is XOR).
          double xtmp = xval * (1 - 2*k);
          double ytmp = yval * (1 - 2*l);
          double ztmp = zval * (1 - 2*(j ^ k ^ l ^ oddPermutation));

          double trace = xtmp + ytmp + ztmp;

          // Find maximum trace
          if (trace > maxtrace)
356
          {
357 358 359 360 361 362 363 364 365
            maxtrace = trace;
            imax = i;
            jmax = j;
            kmax = k;
            lmax = l;
          }
        }
      }
    }
366
  }
367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384

  // Find the permutation to map each column of the orientation
  // matrix to a spatial dimension x, y, or z.
  int xidx = imax;
  int yidx = (imax + 1 + jmax) % 3;
  int zidx = (imax + 2 - jmax) % 3;

  permutation[0] = xidx;
  permutation[1] = yidx;
  permutation[2] = zidx;

  flip[xidx] = kmax;
  flip[yidx] = lmax;
  flip[zidx] = (jmax ^ kmax ^ lmax ^ oddPermutation);
}

//-------------------------------------------------------------------------
// Create an identity string for a file.
385
std::string vtkMINCImageWriterCreateIdentString()
386 387
{
  // A static counter for this process.
388
  static int identx = 1;
389 390

  // The separator between element.
391
  static const char *itemsep = ":";
392 393

  // Get username and hostname
394
#ifdef _WIN32
395 396
  const char *username = nullptr;
  const char *hostname = nullptr;
397 398 399
  char usernametext[100];
  DWORD numchars = sizeof(usernametext);
  if (GetUserName(usernametext, &numchars))
400
  {
401
    username = usernametext;
402
  }
403 404 405
  char hostnametext[100];
  numchars = sizeof(hostnametext);
  if (GetComputerName(hostnametext, &numchars))
406
  {
407
    hostname = hostnametext;
408
  }
409
#else
410 411
  const char *username = getenv("LOGNAME");
  const char *hostname = getenv("HOSTNAME");
412
#endif
413
  if (username == nullptr)
414
  {
415
    username = "nobody";
416
  }
417
  if (hostname == nullptr)
418
  {
419
    hostname = "unknown";
420
  }
421
  std::string ident = username;
422 423 424
  ident.append(itemsep);
  ident.append(hostname);
  ident.append(itemsep);
425

426 427 428 429 430 431 432 433 434
  // Get the local time
  char buf[1024];
  time_t t;
  time(&t);
  strftime(buf, sizeof(buf), "%Y.%m.%d.%H.%M.%S", localtime(&t));
  ident.append(buf);
  ident.append(itemsep);

  // Get the process ID and the counter for this process.
435
#ifdef _WIN32
436
  int processId = GetCurrentProcessId();
437 438 439
#else
  int processId = getpid();
#endif
440
  snprintf(buf, 1024, "%i%s%i", processId, itemsep, identx++);
441
  ident.append(buf);
442

443 444 445 446 447 448 449 450 451 452 453 454
  return ident;
}

//-------------------------------------------------------------------------
nc_type vtkMINCImageWriterConvertVTKTypeToMINCType(
  int dataType,
  int &mincsigned)
{
  nc_type minctype = NC_BYTE;

  // Get the vtk type of the data.
  switch (dataType)
455
  {
456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490
    case VTK_CHAR:
    case VTK_SIGNED_CHAR:
      minctype = NC_BYTE;
      mincsigned = 1;
      break;
    case VTK_UNSIGNED_CHAR:
      minctype = NC_BYTE;
      mincsigned = 0;
      break;
    case VTK_SHORT:
      minctype = NC_SHORT;
      mincsigned = 1;
      break;
    case VTK_UNSIGNED_SHORT:
      minctype = NC_SHORT;
      mincsigned = 0;
      break;
    case VTK_INT:
      minctype = NC_INT;
      mincsigned = 1;
      break;
    case VTK_UNSIGNED_INT:
      minctype = NC_INT;
      mincsigned = 0;
      break;
    case VTK_FLOAT:
      minctype = NC_FLOAT;
      mincsigned = 1;
      break;
    case VTK_DOUBLE:
      minctype = NC_DOUBLE;
      mincsigned = 1;
      break;
    default:
      break;
491
  }
492 493 494 495 496

  return minctype;
}

//-------------------------------------------------------------------------
497 498 499 500 501
// These macros are only for use in WriteMINCFileAttributes.

// Note: Until VTK 7.0, this macro added a terminating null byte to all
// text attributes.  As of VTK 7.1, it does not.  The attribute length
// should be the string length, not the string length "plus one".
502 503
#define vtkMINCImageWriterPutAttributeTextMacro(name, text) \
  if (status == NC_NOERR) \
504
  { \
505
    status = nc_put_att_text(ncid, varid, name, strlen(text), text); \
506
  }
507 508 509

#define vtkMINCImageWriterPutAttributeDoubleMacro(name, count, ptr) \
  if (status == NC_NOERR) \
510
  { \
511
    status = nc_put_att_double(ncid, varid, name, NC_DOUBLE, count, ptr); \
512
  }
513 514 515 516

//-------------------------------------------------------------------------
// Allowed dimension variable names
static const char *vtkMINCDimVarNames[] = {
517 518
  MIxspace, MIyspace, MIzspace, MItime,
  MIxfrequency, MIyfrequency, MIzfrequency, MItfrequency,
519
  nullptr
520 521 522
};

//-------------------------------------------------------------------------
523
int vtkMINCImageWriter::CreateMINCDimensions(
524
  vtkImageData *input, int numTimeSteps, int *dimids)
525
{
526
  int wholeExtent[6];
527 528
  vtkStreamingDemandDrivenPipeline::GetWholeExtent(
    this->GetInputInformation(0, 0), wholeExtent);
529 530
  int numComponents = input->GetNumberOfScalarComponents();

531 532 533
  // Create a default dimension order using the direction cosines.
  this->ComputePermutationFromOrientation(this->Permutation, this->Flip);
  const char *defaultdims[3];
534 535 536
  defaultdims[this->Permutation[0]] = MIxspace;
  defaultdims[this->Permutation[1]] = MIyspace;
  defaultdims[this->Permutation[2]] = MIzspace;
537 538

  int hasTimeDim = 0;
539
  std::vector<std::string> dimensions;
540
  int nuserdims = 0;
541
  vtkStringArray *dimensionNames = nullptr;
542
  if (this->ImageAttributes)
543
  {
544 545
    dimensionNames = this->ImageAttributes->GetDimensionNames();
    nuserdims = dimensionNames->GetNumberOfValues();
546
  }
547
  for (int iuserdims = 0; iuserdims < nuserdims; iuserdims++)
548
  {
549
    const char *dimname = dimensionNames->GetValue(iuserdims);
550
    // Remove vector_dimension, we'll add it back if it is needed
551
    if (strcmp(dimname, MIvector_dimension) == 0)
552
    {
553
      continue;
554
    }
555 556
    // Check for time or tfrequency
    if (dimname[0] == 't')
557
    {
558
      hasTimeDim = 1;
559
    }
560
    // Ensure the dimension name is valid
561 562
    const char **tryname = nullptr;
    for (tryname = vtkMINCDimVarNames; *tryname !=nullptr; tryname++)
563
    {
564
      if (strcmp(dimname, *tryname) == 0)
565
      {
566 567
        break;
      }
568
    }
569
    if (*tryname == nullptr)
570
    {
571 572 573
      vtkErrorMacro("The dimension name " << dimname <<
                    " is not recognized.");
      return 0;
574
    }
575
    // Check for duplicated dimensions
576
    int ndim = static_cast<int>(dimensions.size());
577
    for (int idim = 0; idim < ndim; idim++)
578
    {
579
      if (dimensions[idim][0] == dimname[0])
580
      {
581 582 583 584
        vtkErrorMacro("Tried to create dimension " << dimname <<
                      " but " << dimensions[idim] << " already exists");
        return 0;
      }
585
    }
586 587 588

    // Add the dimension
    dimensions.push_back(dimname);
589
  }
590 591

  // Make sure number of dimensions matches the dimensionality
592
  int timeDimensions = ( numTimeSteps > 1);
593 594 595 596
  int spatialDimensions = ((wholeExtent[0] < wholeExtent[1]) +
                           (wholeExtent[2] < wholeExtent[3]) +
                           (wholeExtent[4] < wholeExtent[5]));
  if (spatialDimensions < 2)
597
  {
598
    spatialDimensions = 2;
599
  }
600 601
  // Insert dimension names until we have all spatial dimensions
  while (static_cast<int>(dimensions.size()) < spatialDimensions+hasTimeDim)
602
  {
603 604
    // Make sure we don't insert a dimension that is already there
    for (int i = 0; i < 3; i++)
605
    {
606
      int idim = 0;
607
      int ndims = static_cast<int>(dimensions.size());
608
      for (idim = 0; idim < ndims; idim++)
609
      {
610
        if (defaultdims[i][0] == dimensions[idim][0])
611
        {
612 613
          break;
        }
614
      }
615
      if (idim == ndims)
616
      {
617 618 619
        dimensions.insert(dimensions.begin(), defaultdims[i]);
      }
    }
620
  }
621 622
  // Make sure we have a time dimension if we need one
  if (timeDimensions == 1 && hasTimeDim == 0)
623
  {
624
    dimensions.insert(dimensions.begin(), MItime);
625
  }
626 627
  // Check for vector_dimension
  if (numComponents > 1)
628
  {
629
    dimensions.push_back(MIvector_dimension);
630
  }
631 632 633 634 635 636 637

  // ------------------------
  // Create the NetCDF dimensions

  int ncid = this->MINCFileId;
  int status = NC_NOERR;

638
  int ndim = static_cast<int>(dimensions.size());
639
  this->FileDimensionNames->SetNumberOfValues(ndim);
640
  for (int idim = 0; idim < ndim; idim++)
641
  {
642
    const char *dimname = dimensions[idim].c_str();
643
    this->FileDimensionNames->SetValue(idim, dimname);
644
    int dimIndex = this->IndexFromDimensionName(dimname);
645
    size_t length = numTimeSteps;
646
    if (dimIndex >= 0 && dimIndex < 3)
647
    {
648
      length = wholeExtent[2*dimIndex+1] - wholeExtent[2*dimIndex] + 1;
649
    }
650
    else if (strcmp(dimname, MIvector_dimension) == 0)
651
    {
652
      length = numComponents;
653
    }
654 655
    status = nc_def_dim(ncid, dimname, length, &dimids[idim]);
    if (status != NC_NOERR)
656
    {
657 658 659 660
      vtkMINCImageWriterFailAndClose(ncid, status);
      this->MINCFileId = 0;
      return 0;
    }
661
  }
662 663 664 665 666

  return 1;
}

//-------------------------------------------------------------------------
667
int vtkMINCImageWriter::CreateMINCVariables(
668
  vtkImageData *input, int vtkNotUsed(numTimeSteps), int *dimids)
669 670 671
{
  // Allowed standard variable names
  static const char *stdVarNames[] = {
672 673
    MIrootvariable, MIimage, MIimagemin, MIimagemax,
    MIpatient, MIstudy, MIacquisition,
674
    nullptr
675
  };
676

677
  std::vector<std::string> variables;
678

679 680 681 682 683 684 685 686
  // Get the information from the input
  double spacing[3];
  double origin[3];
  int wholeExtent[6];
  int numComponents = input->GetNumberOfScalarComponents();
  int imageDataType = input->GetScalarType();
  input->GetSpacing(spacing);
  input->GetOrigin(origin);
687 688
  vtkStreamingDemandDrivenPipeline::GetWholeExtent(
    this->GetInputInformation(0, 0), wholeExtent);
689

690
  // Add all dimensions onto the list of variables
691
  int ndim = this->FileDimensionNames->GetNumberOfValues();
692
  for (int dimidx = 0; dimidx < ndim; dimidx++)
693
  {
694
    const char *dimname = this->FileDimensionNames->GetValue(dimidx);
695
    // vector_dimension isn't ever included as a variable
696
    if (strcmp(dimname, MIvector_dimension) != 0)
697
    {
698
      variables.push_back(this->FileDimensionNames->GetValue(dimidx));
699
    }
700
  }
701
  // Reset ndim so that it only includes dimensions with variables
702
  ndim = static_cast<int>(variables.size());
703

704 705
  variables.push_back(MIimage);
  variables.push_back(MIrootvariable);
706

707
  // Not all MINC images need image-min and image-max.
708
  this->MINCImageMinMaxDims = 0;
709
  if (this->InternalRescaleSlope != 0)
710
  {
711 712 713 714 715
    // Check whether slice-by-slice rescaling is needed
    if ((imageDataType == VTK_FLOAT ||
         imageDataType == VTK_DOUBLE) &&
        (this->MINCImageType != NC_FLOAT &&
         this->MINCImageType != NC_DOUBLE))
716
    {
717
      this->MINCImageMinMaxDims = ndim - 2;
718
    }
719 720
    variables.push_back(MIimagemin);
    variables.push_back(MIimagemax);
721
  }
722 723

  // Add user-defined variables
724
  int nuservars = 0;
725
  vtkStringArray *variableNames = nullptr;
726
  if (this->ImageAttributes)
727
  {
728 729
    variableNames = this->ImageAttributes->GetVariableNames();
    nuservars = variableNames->GetNumberOfValues();
730
  }
731
  for (int iuservars = 0; iuservars < nuservars; iuservars++)
732
  {
733
    const char *varname = variableNames->GetValue(iuservars);
734
    int ivar;
735
    int nvars = static_cast<int>(variables.size());
736
    for (ivar = 0; ivar < nvars; ivar++)
737
    {
738
      if (strcmp(variables[ivar].c_str(), varname) == 0)
739
      {
740 741
        break;
      }
742
    }
743
    if (ivar == nvars) // wasn't already in the list
744
    {
745 746
      // Check if the variable name is a dimension that isn't one
      // of the selected dimensions for this image
747
      for (const char **tryname = vtkMINCDimVarNames; *tryname !=nullptr; tryname++)
748
      {
749
        if (strcmp(varname, *tryname) == 0)
750
        {
751 752 753 754 755
          vtkErrorMacro("The variable " << varname
                        << " is not a dimension of this image");
          return 0;
        }
      }
756
        variables.push_back(varname);
757
    }
758
  }
759 760 761

  // ------------------------
  // Find the children of the root variable
762
  std::string rootChildren;
763

764
  int nvars = static_cast<int>(variables.size());
765 766
  int ivar = 0;
  for (ivar = 0; ivar < nvars; ivar++)
767
  {
768
    const char *varname = variables[ivar].c_str();
769 770 771
    if (strcmp(varname, MIrootvariable) == 0 ||
        strcmp(varname, MIimagemin) == 0 ||
        strcmp(varname, MIimagemax) == 0)
772
    {
773
      continue;
774
    }
775
    for (const char **tryname = stdVarNames; *tryname !=nullptr; tryname++)
776
    {
777
      if (strcmp(varname, *tryname) == 0)
778
      {
779
        if (!rootChildren.empty())
780
        {
781
          rootChildren.append(MI_CHILD_SEPARATOR);
782
        }
783 784 785 786
        rootChildren.append(varname);
        break;
      }
    }
787
  }
788 789 790 791 792 793

  // ------------------------
  // Create the variables and write the attributes.
  // Start at -1, which stands for global attributes.
  int ncid = this->MINCFileId;
  int status = NC_NOERR;
794
  nvars = static_cast<int>(variables.size());
795
  for (ivar = -1; ivar < nvars; ivar++)
796
  {
797 798
    const char *varname = MI_EMPTY_STRING;
    const char *vartype = MI_EMPTY_STRING;
799 800 801
    int varid = -1;

    if (ivar >= 0)
802
    {
803 804
      nc_type cdftype = NC_INT;
      varname = variables[ivar].c_str();
805
      const char *parent = MIrootvariable;
806
      const char *children = nullptr;
807 808 809 810 811
      int vardims = 0;

      // The dimensions are the first variables (note that ndim
      // does not include the vector_dimension)
      if (ivar < ndim)
812
      {
813
        vartype = MI_DIMENSION;
814
      }
815
      else
816
      {
817
        for (const char **tryname = stdVarNames; *tryname != nullptr; tryname++)
818
        {
819
          if (strcmp(varname, *tryname) == 0)
820
          {
821
            vartype = MI_GROUP;
822 823
          }
        }
824
      }
825 826

      // Check if this is an image-related variable
827
      if (strcmp(varname, MIimage) == 0)
828
      {
829 830
        cdftype = (nc_type)this->MINCImageType;
        vardims = ndim + (numComponents > 1);
831
      }
832 833
      else if (strcmp(varname, MIimagemin) == 0 ||
               strcmp(varname, MIimagemax) == 0)
834
      {
835 836
        parent = MIimage;
        vartype = MI_VARATT;
837 838
        cdftype = NC_DOUBLE;
        vardims = this->MINCImageMinMaxDims;
839
      }
840 841

      // Check if this is the rootvariable
842
      if (strcmp(varname, MIrootvariable) == 0)
843
      {
844
        parent = MI_EMPTY_STRING;
845
        children = rootChildren.c_str();
846
      }
847 848 849 850 851 852

      // Create the NetCDF variable
      status = nc_def_var(ncid, varname, cdftype, vardims, dimids,
                          &varid);

      if (status != NC_NOERR)
853
      {
854 855 856
        vtkMINCImageWriterFailAndClose(ncid, status);
        this->MINCFileId = 0;
        return 0;
857
      }
858

859
      // Variables of known type get standard MINC attributes
860
      if (strcmp(vartype, MI_EMPTY_STRING) != 0)
861
      {
862 863 864
        vtkMINCImageWriterPutAttributeTextMacro(MIvarid,    MI_STDVAR);
        vtkMINCImageWriterPutAttributeTextMacro(MIversion,  MI_VERSION_1_0);
        vtkMINCImageWriterPutAttributeTextMacro(MIvartype,  vartype);
865
      }
866 867

      int dimIndex = 0;
868
      if (strcmp(vartype, MI_DIMENSION) == 0)
869
      {
870 871 872 873
        static const char *dimensionComments[] = {
          "X increases from patient left to right",
          "Y increases from patient posterior to anterior",
          "Z increases from patient inferior to superior",
874
          nullptr
875 876 877 878 879 880
        };

        dimIndex = this->IndexFromDimensionName(varname);
        double start = 0.0;
        double step = 1.0;
        if (dimIndex >= 0 && dimIndex < 3)
881
        {
882
          vtkMINCImageWriterPutAttributeTextMacro(
883
            MIcomments, dimensionComments[dimIndex]);
884 885 886
          start = origin[dimIndex];
          step = spacing[dimIndex];
          if (this->Flip[dimIndex])
887
          {
888
            // Switch the MIstart to the other end and change sign
889 890 891 892
            double length = (wholeExtent[2*dimIndex+1] -
                             wholeExtent[2*dimIndex] + 1);
            start = -(start + step*(length-1));
          }
893
        }
894

895 896 897 898 899
        vtkMINCImageWriterPutAttributeDoubleMacro(MIstart, 1, &start);
        vtkMINCImageWriterPutAttributeDoubleMacro(MIstep,  1, &step);
        vtkMINCImageWriterPutAttributeTextMacro(MIspacing,   MI_REGULAR);
        vtkMINCImageWriterPutAttributeTextMacro(MIspacetype, MI_NATIVE);
        vtkMINCImageWriterPutAttributeTextMacro(MIalignment, MI_CENTRE);
900 901 902

        // Extra attributes for spatial dimensions
        if (dimIndex >= 0 && dimIndex < 3)
903
        {
904
          vtkMatrix4x4 *matrix = this->GetDirectionCosines();
905
          if (matrix)
906
          {
907 908 909 910 911 912
            double dircos[3];
            // need to take permutation into account here
            dircos[0] = matrix->GetElement(0, dimIndex);
            dircos[1] = matrix->GetElement(1, dimIndex);
            dircos[2] = matrix->GetElement(2, dimIndex);
            if (this->Flip[dimIndex])
913
            {
914 915
              // Flip the dimension direction
              for (int idx = 0; idx < 3; idx++)
916
              {
917
                if (dircos[idx] != 0)
918
                {
919 920 921
                  dircos[idx] = -dircos[idx];
                }
              }
922
            }
923 924
            vtkMINCImageWriterPutAttributeDoubleMacro(MIdirection_cosines,
                                                      3, dircos);
925 926
          }
        }
927
      }
928
      else if (strcmp(vartype, MI_VARATT) == 0)
929
      {
930
        vtkMINCImageWriterPutAttributeTextMacro(MIparent,  parent);
931
        if (children)
932
        {
933
          vtkMINCImageWriterPutAttributeTextMacro(MIchildren, children);
934
        }
935
        if (strcmp(varname, MIimagemin) == 0)
936
        {
937
          double val = 0.0;
938
          vtkMINCImageWriterPutAttributeDoubleMacro(MI_FillValue, 1, &val);
939
        }
940
        else if (strcmp(varname, MIimagemax) == 0)
941
        {
942
          double val = 1.0;
943
          vtkMINCImageWriterPutAttributeDoubleMacro(MI_FillValue, 1, &val);
944
        }
945
      }
946
      else if (strcmp(vartype, MI_GROUP) == 0)
947
      {
948
        vtkMINCImageWriterPutAttributeTextMacro(MIparent,   parent);
949
        if (children)
950
        {
951
          vtkMINCImageWriterPutAttributeTextMacro(MIchildren, children);
952
        }
953

954
        if (strcmp(varname, MIimage) == 0)
955
        {
956
          const char *signType = MI_SIGNED;
957
          if (this->MINCImageTypeSigned == 0)
958
          {
959
            signType = MI_UNSIGNED;
960
          }
961
          double *validRange = this->FileValidRange;
962

963
          vtkMINCImageWriterPutAttributeTextMacro(MIcomplete,  MI_TRUE);
964 965 966 967

          // Only produce signtype and valid_range for integer data
          if (this->MINCImageType != NC_FLOAT &&
              this->MINCImageType != NC_DOUBLE)
968
          {
969 970 971 972 973
            vtkMINCImageWriterPutAttributeTextMacro(MIsigntype,  signType);

            // Don't set valid_range if the default is suitable
            if (this->ComputeValidRangeFromScalarRange ||
                (this->ImageAttributes &&
974
                 vtkArrayDownCast<vtkDoubleArray>(
975 976
                   this->ImageAttributes->GetAttributeValueAsArray(
                     MIimage, MIvalid_range))))
977
            {
978 979
              vtkMINCImageWriterPutAttributeDoubleMacro(MIvalid_range,2,
                                                   validRange);
980
            }
981
          }
982

983
          // The image-min, image-max will not always be present
984
          if (this->InternalRescaleSlope != 0)
985
          {
986 987 988 989
            vtkMINCImageWriterPutAttributeTextMacro(
              MIimagemin, MI_VARATT_POINTER_PREFIX MIimagemin);
            vtkMINCImageWriterPutAttributeTextMacro(
              MIimagemax, MI_VARATT_POINTER_PREFIX MIimagemax);
990 991 992
          }
        }
      }
993
    }
994
    else
995
    {
996 997
      // Set the varid for global variables
      varid = -1;
998

999
      // Global attributes: ident and history
1000
      std::string ident = vtkMINCImageWriterCreateIdentString();
1001
      vtkMINCImageWriterPutAttributeTextMacro(MIident, ident.c_str());
1002

1003
      // For history, include any previous history
1004
      std::string history;
1005
      const char *previousHistory = nullptr;
1006
      if (this->ImageAttributes)
1007
      {
1008 1009 1010
        previousHistory = this->ImageAttributes->GetAttributeValueAsString(
          MI_EMPTY_STRING, MIhistory);
        if (previousHistory)
1011
        {
1012
          history.append(previousHistory);
1013
        }
1014
      }
1015

1016
      if (history.size() > 1 && history[history.size()-1] != '\n')
1017
      {
1018
        history.append("\n");
1019
      }
1020

1021 1022
      time_t t;
      time(&t);
1023
      std::string timestamp = ctime(&t);
1024 1025
      history.append(timestamp.substr(0, timestamp.size()-1) + ">>>");
      if (this->HistoryAddition)
1026
      {
1027
        history = history + this->HistoryAddition + "\n";
1028
      }
1029
      else
1030
      {
1031
        history = history + "Created by " + this->GetClassName() + "\n";
1032
      }
1033 1034
      vtkMINCImageWriterPutAttributeTextMacro(MIhistory, history.c_str());
    }
1035 1036

    // Write out user-defined attributes for this variable
1037
    vtkStringArray *attArray = nullptr;
1038
    if (this->ImageAttributes)
1039
    {
1040
      attArray = this->ImageAttributes->GetAttributeNames(varname);
1041
    }
1042
    if (attArray)
1043
    {
1044
      std::string varpath = MI_GRPNAME MI_GRP_SEP;
1045 1046
      int natts = attArray->GetNumberOfValues();
      for (int iatt = 0; iatt < natts; iatt++)
1047
      {
1048
        const char *attname = attArray->GetValue(iatt);
1049
        vtkDataArray *array =
1050 1051
          this->ImageAttributes->GetAttributeValueAsArray(
            varname, attname);
1052

1053 1054
        int result = this->ImageAttributes->ValidateAttribute(
          varname, attname, array);
1055 1056

        if (result == 0)
1057
        {
1058 1059 1060
          // A result of zero means that this attribute has already
          // been automatically generated, or is mis-formatted
          continue;
1061
        }
1062
        else if (result > 1 && this->StrictValidation)
1063
        {
1064 1065
          vtkWarningMacro("Attribute " << varname << ":" << attname
                          << " is not recognized");
1066
        }
1067
        else if (strcmp(attname, MIdirection_cosines) == 0 &&
1068
                 this->DirectionCosines)
1069
        {
1070
          // Let DirectionCosines override the attributes setting
1071
          continue;
1072
        }
1073
        else
1074
        {
1075 1076 1077 1078
          // Write out the attribute
          int dataType = array->GetDataType();
          size_t size = array->GetNumberOfTuples();
          switch (dataType)
1079
          {
1080 1081 1082
            case VTK_CHAR:
              status = nc_put_att_text(ncid, varid, attname, size,
                ((vtkCharArray *)array)->GetPointer(0));
1083
              break;
1084 1085 1086
            case VTK_INT:
              status = nc_put_att_int(ncid, varid, attname, NC_INT, size,
                ((vtkIntArray *)array)->GetPointer(0));
1087
              break;
1088 1089 1090 1091 1092
            case VTK_DOUBLE:
              status = nc_put_att_double(ncid, varid, attname, NC_DOUBLE,
                size, ((vtkDoubleArray *)array)->GetPointer(0));
              break;
            default:
1093
            {
1094 1095 1096
              vtkWarningMacro("Attribute " << varname << ":" << attname
                              << " has bad data type " << dataType << ".");
            }