1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkMINCImageReader.cxx
5
6 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 All rights reserved.
8 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9
10 This software is distributed WITHOUT ANY WARRANTY; without even
11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12 PURPOSE. See the above copyright notice for more information.
13
14 =========================================================================*/
15 /*=========================================================================
16
17 Copyright (c) 2006 Atamai, Inc.
18
19 Use, modification and redistribution of the software, in source or
20 binary forms, are permitted provided that the following terms and
21 conditions are met:
22
23 1) Redistribution of the source code, in verbatim or modified
24 form, must retain the above copyright notice, this license,
25 the following disclaimer, and any notices that refer to this
26 license and/or the following disclaimer.
27
28 2) Redistribution in binary form must include the above copyright
29 notice, a copy of this license and the following disclaimer
30 in the documentation or with other materials provided with the
31 distribution.
32
33 3) Modified copies of the source code must be clearly marked as such,
34 and must not be misrepresented as verbatim copies of the source code.
35
36 THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES PROVIDE THE SOFTWARE "AS IS"
37 WITHOUT EXPRESSED OR IMPLIED WARRANTY INCLUDING, BUT NOT LIMITED TO,
38 THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
39 PURPOSE. IN NO EVENT SHALL ANY COPYRIGHT HOLDER OR OTHER PARTY WHO MAY
40 MODIFY AND/OR REDISTRIBUTE THE SOFTWARE UNDER THE TERMS OF THIS LICENSE
41 BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL OR CONSEQUENTIAL DAMAGES
42 (INCLUDING, BUT NOT LIMITED TO, LOSS OF DATA OR DATA BECOMING INACCURATE
43 OR LOSS OF PROFIT OR BUSINESS INTERRUPTION) ARISING IN ANY WAY OUT OF
44 THE USE OR INABILITY TO USE THE SOFTWARE, EVEN IF ADVISED OF THE
45 POSSIBILITY OF SUCH DAMAGES.
46
47 =========================================================================*/
48
49 #include "vtkMINCImageReader.h"
50
51 #include "vtkObjectFactory.h"
52
53 #include "vtkImageData.h"
54 #include "vtkStringArray.h"
55 #include "vtkCharArray.h"
56 #include "vtkUnsignedCharArray.h"
57 #include "vtkShortArray.h"
58 #include "vtkIntArray.h"
59 #include "vtkFloatArray.h"
60 #include "vtkDoubleArray.h"
61 #include "vtkIdTypeArray.h"
62 #include "vtkMatrix4x4.h"
63 #include "vtkSmartPointer.h"
64 #include "vtkMath.h"
65 #include "vtkStreamingDemandDrivenPipeline.h"
66
67 #include "vtkType.h"
68
69 #include "vtkMINCImageAttributes.h"
70 #include "vtkMINC.h"
71 #include "vtk_netcdf.h"
72
73 #include <cstdlib>
74 #include <ctype.h>
75 #include <float.h>
76 #include <string>
77 #include <map>
78
79 #define VTK_MINC_MAX_DIMS 8
80
81 //--------------------------------------------------------------------------
82 vtkStandardNewMacro(vtkMINCImageReader);
83
84 //-------------------------------------------------------------------------
vtkMINCImageReader()85 vtkMINCImageReader::vtkMINCImageReader()
86 {
87 this->NumberOfTimeSteps = 1;
88 this->TimeStep = 0;
89 this->DirectionCosines = vtkMatrix4x4::New();
90 this->RescaleIntercept = 0.0;
91 this->RescaleSlope = 1.0;
92 this->RescaleRealValues = 0;
93
94 this->MINCImageType = 0;
95 this->MINCImageTypeSigned = 1;
96
97 this->ValidRange[0] = 0.0;
98 this->ValidRange[1] = 1.0;
99
100 this->ImageRange[0] = 0.0;
101 this->ImageRange[1] = 1.0;
102
103 this->DataRange[0] = 0.0;
104 this->DataRange[1] = 1.0;
105
106 this->ImageAttributes = vtkMINCImageAttributes::New();
107 this->ImageAttributes->ValidateAttributesOff();
108
109 this->FileNameHasChanged = 0;
110 }
111
112 //-------------------------------------------------------------------------
~vtkMINCImageReader()113 vtkMINCImageReader::~vtkMINCImageReader()
114 {
115 if (this->DirectionCosines)
116 {
117 this->DirectionCosines->Delete();
118 this->DirectionCosines = 0;
119 }
120 if (this->ImageAttributes)
121 {
122 this->ImageAttributes->Delete();
123 this->ImageAttributes = 0;
124 }
125 }
126
127 //-------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)128 void vtkMINCImageReader::PrintSelf(ostream& os, vtkIndent indent)
129 {
130 this->Superclass::PrintSelf(os,indent);
131
132 os << indent << "ImageAttributes: " << this->ImageAttributes << "\n";
133 if (this->ImageAttributes)
134 {
135 this->ImageAttributes->PrintSelf(os, indent.GetNextIndent());
136 }
137 os << indent << "DirectionCosines: " << this->DirectionCosines << "\n";
138 if (this->DirectionCosines)
139 {
140 this->DirectionCosines->PrintSelf(os, indent.GetNextIndent());
141 }
142 os << indent << "RescaleSlope: " << this->RescaleSlope << "\n";
143 os << indent << "RescaleIntercept: " << this->RescaleIntercept << "\n";
144 os << indent << "RescaleRealValues: "
145 << (this->RescaleRealValues ? "On" : "Off") << "\n";
146 os << indent << "DataRange: (" << this->DataRange[0]
147 << ", " << this->DataRange[1] << ")\n";
148
149 os << indent << "NumberOfTimeSteps: " << this->NumberOfTimeSteps << "\n";
150 os << indent << "TimeStep: " << this->TimeStep << "\n";
151 }
152
153 //-------------------------------------------------------------------------
SetFileName(const char * name)154 void vtkMINCImageReader::SetFileName(const char *name)
155 {
156 // Set FileNameHasChanged even if the file name hasn't changed,
157 // because it is possible that the user is re-reading a file after
158 // changing it.
159 if (!(name == 0 && this->GetFileName() == 0))
160 {
161 this->FileNameHasChanged = 1;
162 }
163
164 this->Superclass::SetFileName(name);
165 }
166
167 //-------------------------------------------------------------------------
CanReadFile(const char * fname)168 int vtkMINCImageReader::CanReadFile(const char* fname)
169 {
170 // First do a very rapid check of the magic number
171 FILE *fp = fopen(fname, "rb");
172 if (!fp)
173 {
174 return 0;
175 }
176
177 char magic[4];
178 size_t count = fread(magic, 4, 1, fp);
179 fclose(fp);
180
181 if (count != 1 ||
182 magic[0] != 'C' ||
183 magic[1] != 'D' ||
184 magic[2] != 'F' ||
185 magic[3] != '\001')
186 {
187 return 0;
188 }
189
190 // Do a more thorough check of the image:version attribute, since
191 // there are lots of NetCDF files out there that aren't minc files.
192 int status = NC_NOERR;
193 int ncid = 0;
194 status = nc_open(fname, 0, &ncid);
195 if (status != NC_NOERR)
196 {
197 return 0;
198 }
199
200 int ndims = 0;
201 int nvars = 0;
202 int ngatts = 0;
203 int unlimdimid = 0;
204 status = nc_inq(ncid, &ndims, &nvars, &ngatts, &unlimdimid);
205 if (status != NC_NOERR)
206 {
207 return 0;
208 }
209
210 int varid = 0;
211 char varname[NC_MAX_NAME+1];
212 nc_type vartype = NC_INT;
213 int nvardims;
214 int dimids[VTK_MINC_MAX_DIMS];
215 int nvaratts = 0;
216 for (varid = 0; varid < nvars && status == NC_NOERR; varid++)
217 {
218 status = nc_inq_var(ncid, varid, varname, &vartype, &nvardims,
219 dimids, &nvaratts);
220 if (status == NC_NOERR && strcmp(varname, MIimage) == 0)
221 {
222 nc_type atttype = NC_INT;
223 size_t attlength = 0;
224 status = nc_inq_att(ncid, varid, MIversion, &atttype, &attlength);
225 if (status == NC_NOERR && atttype == NC_CHAR && attlength < 32)
226 {
227 char verstring[32];
228 status = nc_get_att_text(ncid, varid, MIversion, verstring);
229 if (status == NC_NOERR && strncmp(verstring, "MINC ", 5) == 0)
230 {
231 nc_close(ncid);
232 return 1;
233 }
234 }
235 break;
236 }
237 }
238
239 nc_close(ncid);
240
241 return 0;
242 }
243
244 //-------------------------------------------------------------------------
GetDirectionCosines()245 vtkMatrix4x4 *vtkMINCImageReader::GetDirectionCosines()
246 {
247 this->ReadMINCFileAttributes();
248 return this->DirectionCosines;
249 }
250
251 //-------------------------------------------------------------------------
GetRescaleSlope()252 double vtkMINCImageReader::GetRescaleSlope()
253 {
254 this->ReadMINCFileAttributes();
255 this->FindRangeAndRescaleValues();
256 return this->RescaleSlope;
257 }
258
259 //-------------------------------------------------------------------------
GetRescaleIntercept()260 double vtkMINCImageReader::GetRescaleIntercept()
261 {
262 this->ReadMINCFileAttributes();
263 this->FindRangeAndRescaleValues();
264 return this->RescaleIntercept;
265 }
266
267 //-------------------------------------------------------------------------
GetDataRange()268 double *vtkMINCImageReader::GetDataRange()
269 {
270 this->ReadMINCFileAttributes();
271 this->FindRangeAndRescaleValues();
272 return this->DataRange;
273 }
274
275 //-------------------------------------------------------------------------
GetNumberOfTimeSteps()276 int vtkMINCImageReader::GetNumberOfTimeSteps()
277 {
278 this->ReadMINCFileAttributes();
279 return this->NumberOfTimeSteps;
280 }
281
282 //-------------------------------------------------------------------------
GetImageAttributes()283 vtkMINCImageAttributes *vtkMINCImageReader::GetImageAttributes()
284 {
285 this->ReadMINCFileAttributes();
286 return this->ImageAttributes;
287 }
288
289 //-------------------------------------------------------------------------
OpenNetCDFFile(const char * filename,int & ncid)290 int vtkMINCImageReader::OpenNetCDFFile(const char *filename, int& ncid)
291 {
292 int status = 0;
293
294 if (filename == 0)
295 {
296 vtkErrorMacro("No filename was set");
297 return 0;
298 }
299
300 status = nc_open(filename, 0, &ncid);
301 if (status != NC_NOERR)
302 {
303 vtkErrorMacro("Could not open the MINC file:\n"
304 << nc_strerror(status));
305 return 0;
306 }
307
308 return 1;
309 }
310
311 //-------------------------------------------------------------------------
CloseNetCDFFile(int ncid)312 int vtkMINCImageReader::CloseNetCDFFile(int ncid)
313 {
314 int status = 0;
315 status = nc_close(ncid);
316 if (status != NC_NOERR)
317 {
318 vtkErrorMacro("Could not close the MINC file:\n"
319 << nc_strerror(status));
320 return 0;
321 }
322
323 return 1;
324 }
325
326 //-------------------------------------------------------------------------
327 // this is a macro so the vtkErrorMacro will report a useful line number
328 #define vtkMINCImageReaderFailAndClose(ncid, status) \
329 { \
330 if (status != NC_NOERR) \
331 { \
332 vtkErrorMacro("There was an error with the MINC file:\n" \
333 << this->GetFileName() << "\n" \
334 << nc_strerror(status)); \
335 } \
336 nc_close(ncid); \
337 }
338
339
340 //-------------------------------------------------------------------------
341 // Function for getting VTK dimension index from the dimension name.
IndexFromDimensionName(const char * dimName)342 int vtkMINCImageReader::IndexFromDimensionName(const char *dimName)
343 {
344 switch(dimName[0])
345 {
346 case 'x':
347 return 0;
348 case 'y':
349 return 1;
350 case 'z':
351 return 2;
352 default:
353 if (strcmp(dimName, MIvector_dimension) == 0)
354 {
355 return -1;
356 }
357 break;
358 }
359
360 // Any unrecognized dimensions are returned as index 3
361 return 3;
362 }
363
364 //-------------------------------------------------------------------------
ReadMINCFileAttributes()365 int vtkMINCImageReader::ReadMINCFileAttributes()
366 {
367 // If the filename hasn't changed since the last time the attributes
368 // were read, don't read them again.
369 if (!this->FileNameHasChanged)
370 {
371 return 1;
372 }
373
374 // Reset the MINC information for the file.
375 this->MINCImageType = 0;
376 this->MINCImageTypeSigned = 1;
377
378 this->NumberOfTimeSteps = 1;
379 this->DirectionCosines->Identity();
380
381 // Orientation set tells us which direction cosines were found
382 int orientationSet[3];
383 orientationSet[0] = 0;
384 orientationSet[1] = 0;
385 orientationSet[2] = 0;
386
387 this->ImageAttributes->Reset();
388
389 // Miscellaneous NetCDF variables
390 int status = 0;
391 int ncid = 0;
392 int dimid = 0;
393 int varid = 0;
394 int ndims = 0;
395 int nvars = 0;
396 int ngatts = 0;
397 int unlimdimid = 0;
398
399 if (this->OpenNetCDFFile(this->GetFileName(), ncid) == 0)
400 {
401 return 0;
402 }
403
404 // Get the basic information for the file. The ndims are
405 // ignored here, because we only want the dimensions that
406 // belong to the image variable.
407 status = nc_inq(ncid, &ndims, &nvars, &ngatts, &unlimdimid);
408 if (status != NC_NOERR)
409 {
410 vtkMINCImageReaderFailAndClose(ncid, status);
411 return 0;
412 }
413 if (ndims > VTK_MINC_MAX_DIMS)
414 {
415 vtkErrorMacro("MINC file has " << ndims << ", but this reader"
416 " only supports " << VTK_MINC_MAX_DIMS << ".");
417 return 0;
418 }
419
420 // Go through all the variables in the MINC file. A varid of -1
421 // is used to signal global attributes.
422 for (varid = -1; varid < nvars; varid++)
423 {
424 char varname[NC_MAX_NAME+1];
425 int dimids[VTK_MINC_MAX_DIMS];
426 nc_type vartype = NC_SHORT;
427 int nvardims = 0;
428 int nvaratts = 0;
429
430 if (varid == -1) // for global attributes
431 {
432 nvaratts = ngatts;
433 varname[0] = '\0';
434 }
435 else
436 {
437 status = nc_inq_var(ncid, varid, varname, &vartype, &nvardims,
438 dimids, &nvaratts);
439 if (status != NC_NOERR)
440 {
441 vtkMINCImageReaderFailAndClose(ncid, status);
442 return 0;
443 }
444 }
445
446 // Get all the variable attributes
447 for (int j = 0; j < nvaratts; j++)
448 {
449 char attname[NC_MAX_NAME+1];
450 nc_type atttype;
451 size_t attlength = 0;
452
453 status = nc_inq_attname(ncid, varid, j, attname);
454 if (status != NC_NOERR)
455 {
456 vtkMINCImageReaderFailAndClose(ncid, status);
457 return 0;
458 }
459 status = nc_inq_att(ncid, varid, attname, &atttype, &attlength);
460 if (status != NC_NOERR)
461 {
462 vtkMINCImageReaderFailAndClose(ncid, status);
463 return 0;
464 }
465
466 // Get the attribute values as a vtkDataArray.
467 vtkDataArray *dataArray = 0;
468 switch (atttype)
469 {
470 case NC_BYTE:
471 {
472 // NetCDF leaves it up to us to decide whether NC_BYTE
473 // should be signed.
474 vtkUnsignedCharArray *ucharArray = vtkUnsignedCharArray::New();
475 ucharArray->SetNumberOfValues(attlength);
476 nc_get_att_uchar(ncid, varid, attname,
477 ucharArray->GetPointer(0));
478 dataArray = ucharArray;
479 }
480 break;
481 case NC_CHAR:
482 {
483 // The NC_CHAR type is for text.
484 vtkCharArray *charArray = vtkCharArray::New();
485 // The netcdf standard doesn't enforce null-termination
486 // of string attributes, so we add a null here.
487 charArray->SetNumberOfValues(attlength+1);
488 charArray->SetValue(attlength, 0);
489 charArray->SetNumberOfValues(attlength);
490 nc_get_att_text(ncid, varid, attname,
491 charArray->GetPointer(0));
492 dataArray = charArray;
493 }
494 break;
495 case NC_SHORT:
496 {
497 vtkShortArray *shortArray = vtkShortArray::New();
498 shortArray->SetNumberOfValues(attlength);
499 nc_get_att_short(ncid, varid, attname,
500 shortArray->GetPointer(0));
501 dataArray = shortArray;
502 }
503 break;
504 case NC_INT:
505 {
506 vtkIntArray *intArray = vtkIntArray::New();
507 intArray->SetNumberOfValues(attlength);
508 nc_get_att_int(ncid, varid, attname,
509 intArray->GetPointer(0));
510 dataArray = intArray;
511 }
512 break;
513 case NC_FLOAT:
514 {
515 vtkFloatArray *floatArray = vtkFloatArray::New();
516 floatArray->SetNumberOfValues(attlength);
517 nc_get_att_float(ncid, varid, attname,
518 floatArray->GetPointer(0));
519 dataArray = floatArray;
520 }
521 break;
522 case NC_DOUBLE:
523 {
524 vtkDoubleArray *doubleArray = vtkDoubleArray::New();
525 doubleArray->SetNumberOfValues(attlength);
526 nc_get_att_double(ncid, varid, attname,
527 doubleArray->GetPointer(0));
528 dataArray = doubleArray;
529 }
530 break;
531 default:
532 break;
533 }
534 if (dataArray)
535 {
536 this->ImageAttributes->SetAttributeValueAsArray(
537 varname, attname, dataArray);
538 dataArray->Delete();
539 }
540 }
541
542 // Special treatment of image variable.
543 if (strcmp(varname, MIimage) == 0)
544 {
545 // Set the type of the data.
546 this->MINCImageType = vartype;
547
548 // Find the sign of the data, default to "signed"
549 int signedType = 1;
550 // Except for bytes, where default is "unsigned"
551 if (vartype == NC_BYTE)
552 {
553 signedType = 0;
554 }
555 const char *signtype =
556 this->ImageAttributes->GetAttributeValueAsString(
557 MIimage, MIsigntype);
558 if (signtype)
559 {
560 if (strcmp(signtype, MI_UNSIGNED) == 0)
561 {
562 signedType = 0;
563 }
564 }
565 this->MINCImageTypeSigned = signedType;
566
567 for (int i = 0; i < nvardims; i++)
568 {
569 char dimname[NC_MAX_NAME+1];
570 size_t dimlength = 0;
571
572 dimid = dimids[i];
573
574 status = nc_inq_dim(ncid, dimid, dimname, &dimlength);
575 if (status != NC_NOERR)
576 {
577 vtkMINCImageReaderFailAndClose(ncid, status);
578 return 0;
579 }
580
581 this->ImageAttributes->AddDimension(dimname, dimlength);
582
583 int dimIndex = this->IndexFromDimensionName(dimname);
584
585 if (dimIndex >= 0 && dimIndex < 3)
586 {
587 // Set the orientation matrix from the direction_cosines
588 vtkDoubleArray *doubleArray =
589 vtkDoubleArray::SafeDownCast(
590 this->ImageAttributes->GetAttributeValueAsArray(
591 dimname, MIdirection_cosines));
592 if (doubleArray && doubleArray->GetNumberOfTuples() == 3)
593 {
594 double *dimDirCos = doubleArray->GetPointer(0);
595 this->DirectionCosines->SetElement(0, dimIndex, dimDirCos[0]);
596 this->DirectionCosines->SetElement(1, dimIndex, dimDirCos[1]);
597 this->DirectionCosines->SetElement(2, dimIndex, dimDirCos[2]);
598 orientationSet[dimIndex] = 1;
599 }
600 }
601 else if (strcmp(dimname, MIvector_dimension) != 0)
602 {
603 // Set the NumberOfTimeSteps to the product of all dimensions
604 // that are neither spatial dimensions nor vector dimensions.
605 this->NumberOfTimeSteps *= static_cast<int>(dimlength);
606 }
607 }
608 }
609 else if (strcmp(varname, MIimagemin) == 0 ||
610 strcmp(varname, MIimagemax) == 0)
611 {
612 // Read the image-min and image-max.
613 this->ImageAttributes->SetNumberOfImageMinMaxDimensions(nvardims);
614
615 vtkDoubleArray *doubleArray = vtkDoubleArray::New();
616 if (strcmp(varname, MIimagemin) == 0)
617 {
618 this->ImageAttributes->SetImageMin(doubleArray);
619 }
620 else
621 {
622 this->ImageAttributes->SetImageMax(doubleArray);
623 }
624 doubleArray->Delete();
625
626 vtkIdType size = 1;
627 size_t start[VTK_MINC_MAX_DIMS];
628 size_t count[VTK_MINC_MAX_DIMS];
629
630 for (int i = 0; i < nvardims; i++)
631 {
632 char dimname[NC_MAX_NAME+1];
633 size_t dimlength = 0;
634
635 dimid = dimids[i];
636
637 status = nc_inq_dim(ncid, dimid, dimname, &dimlength);
638 if (status != NC_NOERR)
639 {
640 vtkMINCImageReaderFailAndClose(ncid, status);
641 return 0;
642 }
643
644 start[i] = 0;
645 count[i] = dimlength;
646
647 size *= dimlength;
648 }
649
650 doubleArray->SetNumberOfValues(size);
651 status = nc_get_vara_double(ncid, varid, start, count,
652 doubleArray->GetPointer(0));
653 if (status != NC_NOERR)
654 {
655 vtkMINCImageReaderFailAndClose(ncid, status);
656 return 0;
657 }
658 }
659 }
660
661 // Check to see if only 2 spatial dimensions were included,
662 // since we'll have to make up the third dircos if that is the case
663 int numDirCos = 0;
664 int notSetIndex = 0;
665 for (int dcount = 0; dcount < 3; dcount++)
666 {
667 if (orientationSet[dcount])
668 {
669 numDirCos++;
670 }
671 else
672 {
673 notSetIndex = dcount;
674 }
675 }
676 // If only two were set, use cross product to get the third
677 if (numDirCos == 2)
678 {
679 int idx1 = (notSetIndex + 1) % 3;
680 int idx2 = (notSetIndex + 2) % 3;
681 double v1[4];
682 double v2[4];
683 double v3[3];
684 for (int tmpi = 0; tmpi < 4; tmpi++)
685 {
686 v1[tmpi] = v2[tmpi] = 0.0;
687 }
688 v1[idx1] = 1.0;
689 v2[idx2] = 1.0;
690 this->DirectionCosines->MultiplyPoint(v1, v1);
691 this->DirectionCosines->MultiplyPoint(v2, v2);
692 vtkMath::Cross(v1, v2, v3);
693 this->DirectionCosines->SetElement(0, notSetIndex, v3[0]);
694 this->DirectionCosines->SetElement(1, notSetIndex, v3[1]);
695 this->DirectionCosines->SetElement(2, notSetIndex, v3[2]);
696 }
697
698 // Get the data type
699 int dataType = this->ConvertMINCTypeToVTKType(this->MINCImageType,
700 this->MINCImageTypeSigned);
701 this->ImageAttributes->SetDataType(dataType);
702
703 // Get the name from the file name by removing the path and
704 // the extension.
705 const char *fileName = this->FileName;
706 char name[128];
707 name[0] = '\0';
708 int startChar = 0;
709 int endChar = static_cast<int>(strlen(fileName));
710
711 for (startChar = endChar-1; startChar > 0; startChar--)
712 {
713 if (fileName[startChar] == '.')
714 {
715 endChar = startChar;
716 }
717 if (fileName[startChar-1] == '/'
718 #ifdef _WIN32
719 || fileName[startChar-1] == '\\'
720 #endif
721 )
722 {
723 break;
724 }
725 }
726 if (endChar - startChar > 127)
727 {
728 endChar = startChar + 128;
729 }
730 if (endChar > startChar)
731 {
732 strncpy(name, &fileName[startChar], endChar-startChar);
733 name[endChar - startChar] = '\0';
734 }
735
736 this->ImageAttributes->SetName(name);
737
738 // We're done reading the attributes, so close the file.
739 if (this->CloseNetCDFFile(ncid) == 0)
740 {
741 return 0;
742 }
743
744 // Get the ValidRange and ImageRange.
745 this->ImageAttributes->FindValidRange(this->ValidRange);
746 this->ImageAttributes->FindImageRange(this->ImageRange);
747
748 // Don't have to do this again until the file name changes.
749 this->FileNameHasChanged = 0;
750
751 return 1;
752 }
753
754 //-------------------------------------------------------------------------
ConvertMINCTypeToVTKType(int minctype,int mincsigned)755 int vtkMINCImageReader::ConvertMINCTypeToVTKType(
756 int minctype,
757 int mincsigned)
758 {
759 int dataType = 0;
760
761 // Get the vtk type of the data.
762 switch (minctype)
763 {
764 case NC_BYTE:
765 dataType = VTK_UNSIGNED_CHAR;
766 if (mincsigned)
767 {
768 dataType = VTK_SIGNED_CHAR;
769 }
770 break;
771 case NC_SHORT:
772 dataType = VTK_UNSIGNED_SHORT;
773 if (mincsigned)
774 {
775 dataType = VTK_SHORT;
776 }
777 break;
778 case NC_INT:
779 dataType = VTK_UNSIGNED_INT;
780 if (mincsigned)
781 {
782 dataType = VTK_INT;
783 }
784 break;
785 case NC_FLOAT:
786 dataType = VTK_FLOAT;
787 break;
788 case NC_DOUBLE:
789 dataType = VTK_DOUBLE;
790 break;
791 default:
792 break;
793 }
794
795 return dataType;
796 }
797
798 //-------------------------------------------------------------------------
FindRangeAndRescaleValues()799 void vtkMINCImageReader::FindRangeAndRescaleValues()
800 {
801 // Set DataRange and Rescale values according to whether
802 // RescaleRealValues is set
803 if (this->RescaleRealValues)
804 {
805 // Set DataRange to ImageRange
806 this->DataRange[0] = this->ImageRange[0];
807 this->DataRange[1] = this->ImageRange[1];
808
809 // The output data values will be the real data values.
810 this->RescaleSlope = 1.0;
811 this->RescaleIntercept = 0.0;
812 }
813 else
814 {
815 // Set DataRange to ValidRange
816 this->DataRange[0] = this->ValidRange[0];
817 this->DataRange[1] = this->ValidRange[1];
818
819 // Set rescale parameters
820 this->RescaleSlope = ((this->ImageRange[1] - this->ImageRange[0])/
821 (this->ValidRange[1] - this->ValidRange[0]));
822
823 this->RescaleIntercept = (this->ImageRange[0] -
824 this->RescaleSlope*this->ValidRange[0]);
825 }
826 }
827
828 //-------------------------------------------------------------------------
ExecuteInformation()829 void vtkMINCImageReader::ExecuteInformation()
830 {
831 // Read the MINC attributes from the file.
832 if (this->ReadMINCFileAttributes() == 0)
833 {
834 return;
835 }
836
837 // Set the VTK information from the MINC information.
838 int dataExtent[6];
839 dataExtent[0] = dataExtent[1] = 0;
840 dataExtent[2] = dataExtent[3] = 0;
841 dataExtent[4] = dataExtent[5] = 0;
842
843 double dataSpacing[3];
844 dataSpacing[0] = dataSpacing[1] = dataSpacing[2] = 1.0;
845
846 double dataOrigin[3];
847 dataOrigin[0] = dataOrigin[1] = dataOrigin[2] = 0.0;
848
849 int numberOfComponents = 1;
850
851 int fileType = this->ConvertMINCTypeToVTKType(this->MINCImageType,
852 this->MINCImageTypeSigned);
853
854 if (fileType == 0)
855 {
856 vtkErrorMacro("Couldn't convert NetCDF data type " << this->MINCImageType
857 << (this->MINCImageTypeSigned ? " signed" : " unsigned")
858 << " to a VTK data type.");
859 return;
860 }
861
862 // Compute the DataRange, RescaleSlope, and RescaleIntercept
863 this->FindRangeAndRescaleValues();
864
865 // If we are rescaling the data, find the appropriate
866 // output data type. The data is only rescaled if the
867 // data has an ImageMin and ImageMax.
868 int dataType = fileType;
869 if (this->RescaleRealValues &&
870 this->ImageAttributes->GetImageMin() &&
871 this->ImageAttributes->GetImageMax())
872 {
873 switch (fileType)
874 {
875 case VTK_SIGNED_CHAR:
876 case VTK_UNSIGNED_CHAR:
877 case VTK_CHAR:
878 case VTK_SHORT:
879 case VTK_UNSIGNED_SHORT:
880 dataType = VTK_FLOAT;
881 break;
882 case VTK_INT:
883 case VTK_UNSIGNED_INT:
884 dataType = VTK_DOUBLE;
885 break;
886 default:
887 dataType = fileType;
888 break;
889 }
890 }
891
892 // Go through the image dimensions to discover data information.
893 vtkStringArray *dimensionNames =
894 this->ImageAttributes->GetDimensionNames();
895 vtkIdTypeArray *dimensionLengths =
896 this->ImageAttributes->GetDimensionLengths();
897
898 int numberOfDimensions = dimensionNames->GetNumberOfValues();
899 for (int i = 0; i < numberOfDimensions; i++)
900 {
901 const char *dimName = dimensionNames->GetValue(i);
902 vtkIdType dimLength = dimensionLengths->GetValue(i);
903
904 // Set the VTK dimension index.
905 int dimIndex = this->IndexFromDimensionName(dimName);
906
907 // Do special things with the spatial dimensions.
908 if (dimIndex >= 0 && dimIndex < 3)
909 {
910 // Set the spacing from the 'step' attribute.
911 double step = this->ImageAttributes->GetAttributeValueAsDouble(
912 dimName, MIstep);
913 if (step)
914 {
915 dataSpacing[dimIndex] = step;
916 }
917
918 // Set the origin from the 'start' attribute.
919 double start = this->ImageAttributes->GetAttributeValueAsDouble(
920 dimName, MIstart);
921 if (start)
922 {
923 dataOrigin[dimIndex] = start;
924 }
925
926 // Set the extent from the dimension length.
927 dataExtent[2*dimIndex + 1] = static_cast<int>(dimLength - 1);
928 }
929
930 // Check for vector_dimension.
931 else if (strcmp(dimName, MIvector_dimension) == 0)
932 {
933 numberOfComponents = dimLength;
934 }
935 }
936
937 this->SetDataExtent(dataExtent);
938 this->SetDataSpacing(dataSpacing[0], dataSpacing[1], dataSpacing[2]);
939 this->SetDataOrigin(dataOrigin[0], dataOrigin[1], dataOrigin[2]);
940 this->SetDataScalarType(dataType);
941 this->SetNumberOfScalarComponents(numberOfComponents);
942 }
943
944 //-------------------------------------------------------------------------
945 // Data conversion functions. The rounding is done using the same
946 // method as in the MINC libraries.
947 #define vtkMINCImageReaderConvertMacro(F, T, MIN, MAX) \
948 inline void vtkMINCImageReaderConvert(const F& inVal, T& outVal) \
949 { \
950 double val = inVal; \
951 if (val >= static_cast<double>(MIN)) \
952 { \
953 if (val <= static_cast<double>(MAX)) \
954 { \
955 outVal = static_cast<T>((val < 0) ? (val - 0.5) : (val + 0.5)); \
956 return; \
957 } \
958 outVal = static_cast<T>(MAX); \
959 return; \
960 } \
961 outVal = static_cast<T>(MIN); \
962 }
963
964 #define vtkMINCImageReaderConvertMacroFloat(F, T) \
965 inline void vtkMINCImageReaderConvert(const F &inVal, T &outVal) \
966 { \
967 outVal = static_cast<T>(inVal); \
968 }
969
970 vtkMINCImageReaderConvertMacro(double, signed char,
971 VTK_SIGNED_CHAR_MIN, VTK_SIGNED_CHAR_MAX);
972 vtkMINCImageReaderConvertMacro(double, unsigned char,
973 0, VTK_UNSIGNED_CHAR_MAX);
974 vtkMINCImageReaderConvertMacro(double, short,
975 VTK_SHORT_MIN, VTK_SHORT_MAX);
976 vtkMINCImageReaderConvertMacro(double, unsigned short,
977 0, VTK_UNSIGNED_SHORT_MAX);
978 vtkMINCImageReaderConvertMacro(double, int,
979 VTK_INT_MIN, VTK_INT_MAX);
980 vtkMINCImageReaderConvertMacro(double, unsigned int,
981 0, VTK_UNSIGNED_INT_MAX);
982 vtkMINCImageReaderConvertMacroFloat(double, float);
983 vtkMINCImageReaderConvertMacroFloat(double, double);
984
985 //-------------------------------------------------------------------------
986 // Overloaded functions for reading various data types.
987
988 // Handle most with a macro.
989 #define vtkMINCImageReaderReadChunkMacro(ncFunction, T) \
990 inline int vtkMINCImageReaderReadChunk( \
991 int ncid, int varid, size_t *start, size_t *count, T *buffer) \
992 { \
993 return ncFunction(ncid, varid, start, count, buffer); \
994 }
995
996 #define vtkMINCImageReaderReadChunkMacro2(ncFunction, T1, T2) \
997 inline int vtkMINCImageReaderReadChunk( \
998 int ncid, int varid, size_t *start, size_t *count, T1 *buffer) \
999 { \
1000 return ncFunction(ncid, varid, start, count, (T2 *)buffer); \
1001 }
1002
1003 vtkMINCImageReaderReadChunkMacro(nc_get_vara_schar, signed char);
1004 vtkMINCImageReaderReadChunkMacro(nc_get_vara_uchar, unsigned char);
1005 vtkMINCImageReaderReadChunkMacro(nc_get_vara_short, short);
1006 vtkMINCImageReaderReadChunkMacro2(nc_get_vara_short, unsigned short, short);
1007 vtkMINCImageReaderReadChunkMacro(nc_get_vara_int, int);
1008 vtkMINCImageReaderReadChunkMacro2(nc_get_vara_int, unsigned int, int);
1009 vtkMINCImageReaderReadChunkMacro(nc_get_vara_float, float);
1010 vtkMINCImageReaderReadChunkMacro(nc_get_vara_double, double);
1011
1012 //-------------------------------------------------------------------------
1013 template<class T1, class T2>
vtkMINCImageReaderExecuteChunk(T1 * outPtr,T2 * buffer,double slope,double intercept,int ncid,int varid,int ndims,size_t * start,size_t * count,vtkIdType * permutedInc)1014 void vtkMINCImageReaderExecuteChunk(
1015 T1 *outPtr, T2 *buffer, double slope, double intercept,
1016 int ncid, int varid, int ndims, size_t *start, size_t *count,
1017 vtkIdType *permutedInc)
1018 {
1019 // Read the chunk of data from the MINC file.
1020 vtkMINCImageReaderReadChunk(ncid, varid, start, count, buffer);
1021
1022 // Create space to save values during the copy loop.
1023 T1 *saveOutPtr[VTK_MINC_MAX_DIMS];
1024 size_t index[VTK_MINC_MAX_DIMS];
1025 int idim = 0;
1026 for (idim = 0; idim < ndims; idim++)
1027 {
1028 index[idim] = 0;
1029 saveOutPtr[idim] = outPtr;
1030 }
1031
1032 // See if there is a range of dimensions over which the
1033 // the MINC data and VTK data will be contiguous. The
1034 // lastdim is the dimension after which all dimensions
1035 // are contiguous between the MINC file and the output.
1036 int lastdim = ndims-1;
1037 int ncontiguous = 1;
1038 vtkIdType dimprod = 1;
1039 for (idim = ndims; idim > 0; )
1040 {
1041 idim--;
1042
1043 lastdim = idim;
1044 ncontiguous = dimprod;
1045
1046 if (dimprod != permutedInc[idim])
1047 {
1048 break;
1049 }
1050
1051 dimprod *= count[idim];
1052 }
1053
1054 // Save the count and permuted increment of this dimension.
1055 size_t lastdimcount = count[lastdim];
1056 size_t lastdimindex = 0;
1057 vtkIdType lastdimInc = permutedInc[lastdim];
1058 T1 *lastdimOutPtr = saveOutPtr[lastdim];
1059
1060 // Loop over all contiguous sections of the image.
1061 for (;;)
1062 {
1063 // Loop through one contiguous section
1064 vtkIdType k = ncontiguous;
1065 do
1066 {
1067 // Use special function for type conversion.
1068 vtkMINCImageReaderConvert((*buffer++)*slope + intercept, *outPtr++);
1069 }
1070 while (--k);
1071
1072 lastdimindex++;
1073 lastdimOutPtr += lastdimInc;
1074 outPtr = lastdimOutPtr;
1075
1076 // Continue until done lastdim.
1077 if (lastdimindex < lastdimcount)
1078 {
1079 continue;
1080 }
1081
1082 // Handle all dimensions that are lower than lastdim. Go down
1083 // the dimensions one at a time until we find one for which
1084 // the index is still less than the count.
1085 idim = lastdim;
1086 do
1087 {
1088 // We're done if the lowest dim's index has reached its count.
1089 if (idim == 0)
1090 {
1091 return;
1092 }
1093 // Reset the index to zero if it previously reached its count.
1094 index[idim--] = 0;
1095
1096 // Now increase the index for the next lower dimension;
1097 index[idim]++;
1098 saveOutPtr[idim] += permutedInc[idim];
1099
1100 // Continue the loop if this dim's index has reached its count.
1101 }
1102 while (index[idim] >= count[idim]);
1103
1104 // Increment back up to the lastdim, resetting the pointers.
1105 outPtr = saveOutPtr[idim];
1106 do
1107 {
1108 saveOutPtr[++idim] = outPtr;
1109 }
1110 while (idim < lastdim);
1111
1112 lastdimOutPtr = outPtr;
1113 lastdimindex = 0;
1114 }
1115 }
1116
1117 //-------------------------------------------------------------------------
1118 // Our own template that only includes MINC data types.
1119
1120 #define vtkMINCImageReaderTemplateMacro(call) \
1121 case VTK_DOUBLE: { typedef double VTK_TT; call; }; break; \
1122 case VTK_FLOAT: { typedef float VTK_TT; call; }; break; \
1123 case VTK_INT: { typedef int VTK_TT; call; }; break; \
1124 case VTK_UNSIGNED_INT: { typedef unsigned int VTK_TT; call; }; break; \
1125 case VTK_SHORT: { typedef short VTK_TT; call; }; break; \
1126 case VTK_UNSIGNED_SHORT: { typedef unsigned short VTK_TT; call; }; break; \
1127 case VTK_SIGNED_CHAR: { typedef signed char VTK_TT; call; }; break; \
1128 case VTK_UNSIGNED_CHAR: { typedef unsigned char VTK_TT; call; }; break
1129
1130 //-------------------------------------------------------------------------
ExecuteDataWithInformation(vtkDataObject * output,vtkInformation * outInfo)1131 void vtkMINCImageReader::ExecuteDataWithInformation(vtkDataObject *output,
1132 vtkInformation *outInfo)
1133 {
1134 vtkImageData *data = this->AllocateOutputData(output, outInfo);
1135 int scalarType = data->GetScalarType();
1136 int scalarSize = data->GetScalarSize();
1137 int numComponents = data->GetNumberOfScalarComponents();
1138 int outExt[6];
1139 memcpy(outExt, vtkStreamingDemandDrivenPipeline::GetUpdateExtent(
1140 this->GetOutputInformation(0)), 6*sizeof(int));
1141 vtkIdType outInc[3];
1142 data->GetIncrements(outInc);
1143 int outSize[3];
1144 data->GetDimensions(outSize);
1145
1146 void *outPtr = data->GetScalarPointerForExtent(outExt);
1147
1148 int timeStep = this->TimeStep;
1149 if (timeStep < 0 || timeStep >= this->NumberOfTimeSteps)
1150 {
1151 vtkWarningMacro("TimeStep is set to " << this->TimeStep <<
1152 " but there are only " << this->NumberOfTimeSteps <<
1153 " time steps.");
1154 timeStep = timeStep % this->NumberOfTimeSteps;
1155 }
1156
1157 int status = 0;
1158 int ncid = 0;
1159 int varid = 0;
1160
1161 if (this->OpenNetCDFFile(this->GetFileName(), ncid) == 0)
1162 {
1163 return;
1164 }
1165
1166 // Get the image variable.
1167 status = nc_inq_varid(ncid, MIimage, &varid);
1168 if (status != NC_NOERR)
1169 {
1170 vtkMINCImageReaderFailAndClose(ncid, status);
1171 return;
1172 }
1173
1174 // Get the dimensions.
1175 vtkStringArray *dimensionNames =
1176 this->ImageAttributes->GetDimensionNames();
1177 vtkIdTypeArray *dimensionLengths =
1178 this->ImageAttributes->GetDimensionLengths();
1179 int ndims = dimensionNames->GetNumberOfValues();
1180 int idim = 0;
1181 int nminmaxdims = this->ImageAttributes->GetNumberOfImageMinMaxDimensions();
1182 vtkIdType minmaxSize = 0;
1183 if (this->ImageAttributes->GetImageMin())
1184 {
1185 minmaxSize = this->ImageAttributes->GetImageMin()->GetNumberOfTuples();
1186 }
1187
1188 // The default dimensionality of the chunks that are used.
1189 int nchunkdims = ndims - nminmaxdims;
1190
1191 // All of these values will be changed in the following loop
1192 vtkIdType nchunks = 1;
1193 vtkIdType numTimeSteps = 1;
1194 vtkIdType chunkSize = 1;
1195 int hitChunkSizeLimit = 0;
1196 int nchunkdimsIsSet = 0;
1197
1198 // These arrays will be filled in by the following loop
1199 vtkIdType permutedInc[VTK_MINC_MAX_DIMS];
1200 size_t start[VTK_MINC_MAX_DIMS];
1201 size_t count[VTK_MINC_MAX_DIMS];
1202 size_t length[VTK_MINC_MAX_DIMS];
1203
1204 // Loop over the dimensions starting with the fastest-varying.
1205 for (idim = ndims; idim > 0; )
1206 {
1207 idim--;
1208
1209 const char *dimName = dimensionNames->GetValue(idim);
1210 vtkIdType dimLength = dimensionLengths->GetValue(idim);
1211 length[idim] = dimLength;
1212
1213 // Find the VTK dimension index.
1214 int dimIndex = this->IndexFromDimensionName(dimName);
1215
1216 if (dimIndex >= 0 && dimIndex < 3)
1217 {
1218 // Set start and count according to the update extent.
1219 start[idim] = outExt[2*dimIndex];
1220 count[idim] = outExt[2*dimIndex+1] - outExt[2*dimIndex] + 1;
1221 permutedInc[idim] = outInc[dimIndex];
1222 }
1223 else if (strcmp(dimName, MIvector_dimension) == 0)
1224 {
1225 // Vector dimension size is also stored in numComponents.
1226 start[idim] = 0;
1227 count[idim] = numComponents;
1228 permutedInc[idim] = 1;
1229 }
1230 else
1231 {
1232 // Use TimeStep to compute the index into the remaining dimensions.
1233 start[idim] = (timeStep / numTimeSteps) % dimLength;
1234 count[idim] = 1;
1235 numTimeSteps *= dimLength;
1236 permutedInc[idim] = 0;
1237 }
1238
1239 // For scalar minmax, use chunk sizes of 65536 or less,
1240 // unless this would force the chunk size to be 1
1241 if (nminmaxdims == 0 && chunkSize != 1 &&
1242 chunkSize*count[idim] > 65536)
1243 {
1244 hitChunkSizeLimit = 1;
1245 }
1246
1247 // If idim is one of the image-min/image-max dimensions, or if
1248 // we have reached the maximum chunk size, then increase the
1249 // number of chunks instead of increasing the chunk size
1250 if (idim < nminmaxdims || hitChunkSizeLimit)
1251 {
1252 // Number of chunks is product of dimensions in minmax.
1253 nchunks *= count[idim];
1254
1255 // Only set nchunkdims once
1256 if (nchunkdimsIsSet == 0)
1257 {
1258 nchunkdims = ndims - idim - 1;
1259 nchunkdimsIsSet = 1;
1260 }
1261 }
1262 else
1263 {
1264 chunkSize *= count[idim];
1265 }
1266 }
1267
1268 // Create a buffer for intermediate results.
1269 int fileType = this->ImageAttributes->GetDataType();
1270 void *buffer = 0;
1271 switch (fileType)
1272 {
1273 vtkMINCImageReaderTemplateMacro(buffer=(void *)(new VTK_TT[chunkSize]));
1274 }
1275
1276 // Initialize the min and max to the global min max.
1277 double *minPtr = &this->ImageRange[0];
1278 double *maxPtr = &this->ImageRange[1];
1279
1280 // If min and max arrays are not empty, use them instead.
1281 if (minmaxSize > 0)
1282 {
1283 minPtr = this->ImageAttributes->GetImageMin()->GetPointer(0);
1284 maxPtr = this->ImageAttributes->GetImageMax()->GetPointer(0);
1285 }
1286
1287 // Initialize the start and count to use for each chunk.
1288 size_t start2[VTK_MINC_MAX_DIMS];
1289 size_t count2[VTK_MINC_MAX_DIMS];
1290 for (idim = 0; idim < ndims; idim++)
1291 {
1292 start2[idim] = start[idim];
1293 count2[idim] = count[idim];
1294 }
1295
1296 // Go through all the chunks
1297 for (vtkIdType ichunk = 0; ichunk < nchunks; ichunk++)
1298 {
1299 // Find the start and count to use for each chunk.
1300 vtkIdType minmaxIdx = 0;
1301 vtkIdType minmaxInc = 1;
1302 vtkIdType chunkProd = 1;
1303 vtkIdType chunkOffset = 0;
1304 for (idim = ndims - nchunkdims; idim > 0; )
1305 {
1306 idim--;
1307 start2[idim] = start[idim] + (ichunk / chunkProd) % count[idim];
1308 count2[idim] = 1;
1309 if (idim < nminmaxdims)
1310 {
1311 minmaxIdx += start2[idim]*minmaxInc;
1312 minmaxInc *= length[idim];
1313 }
1314 chunkOffset += (start2[idim] - start[idim])*permutedInc[idim];
1315 chunkProd *= count[idim];
1316 }
1317
1318 // Get the min and max values to apply to this chunk
1319 double chunkRange[2];
1320 if (fileType == VTK_FLOAT || fileType == VTK_DOUBLE)
1321 {
1322 // minc files that are float or double use global scaling
1323 chunkRange[0] = this->ImageRange[0];
1324 chunkRange[1] = this->ImageRange[1];
1325 }
1326 else
1327 {
1328 // minc files of other types use slice-by-slice scaling
1329 chunkRange[0] = minPtr[minmaxIdx];
1330 chunkRange[1] = maxPtr[minmaxIdx];
1331 }
1332
1333 // Use the range to calculate a linear transformation
1334 // to apply to the data values of this chunk.
1335 double slope = ((chunkRange[1] - chunkRange[0])/
1336 ((this->ValidRange[1] - this->ValidRange[0])
1337 *this->RescaleSlope));
1338 double intercept = ((chunkRange[0] - this->RescaleIntercept)/
1339 this->RescaleSlope) - slope*this->ValidRange[0];
1340
1341 // set the output pointer to use for this chunk
1342 void *outPtr1 = (void *)(((char *)outPtr) + chunkOffset*scalarSize);
1343
1344 // Read in the chunks and permute them.
1345 if (scalarType == fileType)
1346 {
1347 switch (scalarType)
1348 {
1349 vtkMINCImageReaderTemplateMacro(
1350 vtkMINCImageReaderExecuteChunk(
1351 (VTK_TT *)outPtr1, (VTK_TT *)buffer, slope, intercept,
1352 ncid, varid, ndims, start2, count2, permutedInc));
1353 }
1354 }
1355 else if (scalarType == VTK_FLOAT)
1356 {
1357 switch (fileType)
1358 {
1359 vtkMINCImageReaderTemplateMacro(
1360 vtkMINCImageReaderExecuteChunk(
1361 (float *)outPtr1, (VTK_TT *)buffer, slope, intercept,
1362 ncid, varid, ndims, start2, count2, permutedInc));
1363 }
1364 }
1365 else if (scalarType == VTK_DOUBLE)
1366 {
1367 switch (fileType)
1368 {
1369 vtkMINCImageReaderTemplateMacro(
1370 vtkMINCImageReaderExecuteChunk(
1371 (double *)outPtr1, (VTK_TT *)buffer, slope, intercept,
1372 ncid, varid, ndims, start2, count2, permutedInc));
1373 }
1374 }
1375 }
1376
1377 switch (fileType)
1378 {
1379 vtkMINCImageReaderTemplateMacro(delete [] ((VTK_TT *)buffer));
1380 }
1381
1382 this->CloseNetCDFFile(ncid);
1383 }
1384