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
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
|
#include "fileio.h"
#ifdef VTKSUPPORT
#include <vtkStructuredPoints.h>
#include <vtkStructuredPointsReader.h>
#include <vtkStructuredPointsWriter.h>
//////////////////////////////////////////////////////////////
struct VtkFormat : public FileFormat {
STD_string description() const {return "Visualization Toolkit, vtkStructuredPoints";}
svector suffix() const {
svector result; result.resize(1);
result[0]="vtk";
return result;
}
svector dialects() const {return svector();}
int read(Data<float,4>& data, const STD_string& filename, const FileReadOpts& opts, Protocol& prot) {
Log<FileIO> odinlog("VtkFormat","read");
vtkStructuredPointsReader *in_read=vtkStructuredPointsReader::New();
vtkStructuredPoints *in=vtkStructuredPoints::New();
in_read->SetOutput(in);
in_read->SetFileName(filename.c_str());
if(!in_read->IsFileStructuredPoints()) {
ODINLOG(odinlog,errorLog) << "Not a valid vtkStructuredPoints file" << STD_endl;
return -1;
}
in_read->Update(); // read data
int dims[3];
in->GetDimensions(dims);
int nread=dims[0];
int nphase=dims[1];
int nslice=dims[2];
ODINLOG(odinlog,normalDebug) << "nslice/nphase/nread " << nslice << "/" << nphase << "/" << nread << STD_endl;
ODINLOG(odinlog,normalDebug) << "type=" << in->GetScalarTypeAsString() << STD_endl;
data.resize(1, nslice, nphase, nread);
for(int i=0;i<nread;i++) {
for(int j=0;j<nphase;j++) {
for(int k=0;k<nslice;k++) {
data(0,k,j,i)=in->GetScalarComponentAsFloat(i,j,k,0);
}
}
}
Geometry& geo=prot.geometry;
double spacing[3];
in->GetSpacing(spacing);
geo.set_FOV(readDirection, nread*spacing[0]);
geo.set_FOV(phaseDirection, nphase*spacing[1]);
geo.set_FOV(sliceDirection, nslice*spacing[2]);
geo.set_sliceThickness(spacing[2]);
geo.set_sliceDistance(spacing[2]); //no interslice distance in vtk - so use thickness
in->Delete();
in_read->Delete();
return nslice;
}
int write(const Data<float,4>& data, const STD_string& filename, const FileWriteOpts& opts, const Protocol& prot) {
Log<FileIO> odinlog("VtkFormat","write");
vtkStructuredPointsWriter *out_write=vtkStructuredPointsWriter::New();
vtkStructuredPoints *out=vtkStructuredPoints::New();
#if VTK_MAJOR_VERSION > 5
out_write->SetInputData(out);
#else
out_write->SetInput(out);
#endif
int nrep=data.extent(timeDim);
int nslice=data.extent(sliceDim);
int nphase=data.extent(phaseDim);
int nread=data.extent(readDim);
ODINLOG(odinlog,normalDebug) << "nrep/nslice/nphase/nread "<< nrep << "/" << nslice << "/" << nphase << "/" << nread << STD_endl;
const Geometry& geo=prot.geometry;
out->SetDimensions(nread,nphase,nslice);
out->SetSpacing(voxel_extent(geo, readDirection, nread), voxel_extent(geo, phaseDirection, nphase),voxel_extent(geo, sliceDirection, nslice));
out->SetOrigin(0,0,0);
// setting format (and allocate data on vtk6)
STD_string datatype=select_write_datatype(prot,opts);
#if VTK_MAJOR_VERSION > 5
if(datatype==TypeTraits::type2label((float)0)) out->AllocateScalars(VTK_FLOAT,1);
if(datatype==TypeTraits::type2label((double)0)) out->AllocateScalars(VTK_DOUBLE,1);
if(datatype==TypeTraits::type2label((s32bit)0)) out->AllocateScalars(VTK_INT,1);
if(datatype==TypeTraits::type2label((u32bit)0)) out->AllocateScalars(VTK_UNSIGNED_INT,1);
if(datatype==TypeTraits::type2label((s16bit)0)) out->AllocateScalars(VTK_SHORT,1);
if(datatype==TypeTraits::type2label((u16bit)0)) out->AllocateScalars(VTK_UNSIGNED_SHORT,1);
if(datatype==TypeTraits::type2label((s8bit)0)) out->AllocateScalars(VTK_SIGNED_CHAR,1);
if(datatype==TypeTraits::type2label((u8bit)0)) out->AllocateScalars(VTK_UNSIGNED_CHAR,1);
#else
if(datatype==TypeTraits::type2label((float)0)) out->SetScalarTypeToFloat();
if(datatype==TypeTraits::type2label((double)0)) out->SetScalarTypeToDouble();
if(datatype==TypeTraits::type2label((s32bit)0)) out->SetScalarTypeToInt();
if(datatype==TypeTraits::type2label((u32bit)0)) out->SetScalarTypeToUnsignedInt();
if(datatype==TypeTraits::type2label((s16bit)0)) out->SetScalarTypeToShort();
if(datatype==TypeTraits::type2label((u16bit)0)) out->SetScalarTypeToUnsignedShort();
if(datatype==TypeTraits::type2label((s8bit)0)) out->SetScalarTypeToSignedChar();
if(datatype==TypeTraits::type2label((u8bit)0)) out->SetScalarTypeToUnsignedChar();
out->SetNumberOfScalarComponents(1);
#endif
LDRfileName fname(filename);
for(int irep=0; irep<nrep; irep++) {
// assign data to vtk array
for(int i=0;i<nread;i++) {
for(int j=0;j<nphase;j++) {
for(int k=0;k<nslice;k++) {
out->SetScalarComponentFromFloat(i,j,k,0,data(irep,k,j,i));
}
}
}
STD_string onefilename=fname.get_dirname()+SEPARATOR_STR+fname.get_basename_nosuffix();
if(nrep>1) onefilename+="_time"+itos(irep,nrep-1);
onefilename+="."+fname.get_suffix();
ODINLOG(odinlog,normalDebug) << "writing file "<< onefilename << STD_endl;
out_write->SetFileName(onefilename.c_str());
out_write->SetScalarsName(onefilename.c_str());
// if(ascii) out_write->SetFileTypeToASCII();
// else out_write->SetFileTypeToBinary();
ostream *fp=out_write->OpenVTKFile();
out_write->Write();
out_write->CloseVTKFile(fp);
}
out->Delete();
out_write->Delete();
return nslice;
}
};
#endif
//////////////////////////////////////////////////////////////
void register_vtk_format() {
#ifdef VTKSUPPORT
static VtkFormat vf;
vf.register_format();
#endif
}
|
