1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkGaussianCubeReader.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 #include "vtkGaussianCubeReader.h"
16
17 #include "vtkImageData.h"
18 #include "vtkInformation.h"
19 #include "vtkInformationVector.h"
20 #include "vtkPointData.h"
21 #include "vtkPoints.h"
22 #include "vtkFloatArray.h"
23 #include "vtkIdTypeArray.h"
24 #include "vtkObjectFactory.h"
25 #include "vtkStreamingDemandDrivenPipeline.h"
26 #include "vtkStringArray.h"
27 #include "vtkTransform.h"
28
29 #include <ctype.h>
30
31 vtkStandardNewMacro(vtkGaussianCubeReader);
32
33 //----------------------------------------------------------------------------
34 // Construct object with merging set to true.
vtkGaussianCubeReader()35 vtkGaussianCubeReader::vtkGaussianCubeReader()
36 {
37 this->FileName = NULL;
38 this->Transform = vtkTransform::New();
39 // Add the second output for the grid data
40
41 this->SetNumberOfOutputPorts(2);
42 vtkImageData *grid;
43 grid = vtkImageData::New();
44 grid->ReleaseData();
45 this->GetExecutive()->SetOutputData(1, grid);
46 grid->Delete();
47 }
48
49 //----------------------------------------------------------------------------
~vtkGaussianCubeReader()50 vtkGaussianCubeReader::~vtkGaussianCubeReader()
51 {
52 delete [] this->FileName;
53
54 this->Transform->Delete();
55 // must delete the second output added
56 }
57
58 //----------------------------------------------------------------------------
RequestData(vtkInformation * vtkNotUsed (request),vtkInformationVector ** vtkNotUsed (inputVector),vtkInformationVector * outputVector)59 int vtkGaussianCubeReader::RequestData(
60 vtkInformation *vtkNotUsed(request),
61 vtkInformationVector **vtkNotUsed(inputVector),
62 vtkInformationVector *outputVector)
63 {
64 vtkInformation *outInfo = outputVector->GetInformationObject(0);
65 vtkPolyData *output = vtkPolyData::SafeDownCast(
66 outInfo->Get(vtkDataObject::DATA_OBJECT()));
67
68 FILE *fp;
69 char title[256];
70 char data_name[256];
71 double elements[16];
72 int JN1, N1N2, n1, n2, n3, i, j, k;
73 float tmp, *cubedata;
74 bool orbitalCubeFile = false;
75 int numberOfOrbitals;
76
77 // Output 0 (the default is the polydata)
78 // Output 1 will be the gridded Image data
79
80 vtkImageData *grid = this->GetGridOutput();
81
82 if (!this->FileName)
83 {
84 return 0;
85 }
86
87 if ((fp = fopen(this->FileName, "r")) == NULL)
88 {
89 vtkErrorMacro(<< "File " << this->FileName << " not found");
90 return 0;
91 }
92
93 if (!fgets(title, 256, fp))
94 {
95 vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
96 << " Premature EOF while reading title.");
97 fclose (fp);
98 return 0;
99 }
100 if(strtok(title, ":") != NULL)
101 {
102 if(strtok(NULL, ":") != NULL)
103 {
104 strcpy(data_name, strtok(NULL, ":"));
105 fprintf(stderr,"label = %s\n", data_name);
106 }
107 }
108 if (!fgets(title, 256, fp))
109 {
110 vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
111 << " Premature EOF while reading title.");
112 fclose (fp);
113 return 0;
114 }
115
116 // Read in number of atoms, x-origin, y-origin z-origin
117 //
118 if (fscanf(fp, "%d %lf %lf %lf", &(this->NumberOfAtoms), &elements[3], &elements[7], &elements[11]) != 4)
119 {
120 vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
121 << " Premature EOF while reading atoms, x-origin y-origin z-origin.");
122 fclose (fp);
123 return 0;
124 }
125 if(this->NumberOfAtoms < 0 )
126 {
127 this->NumberOfAtoms = -this->NumberOfAtoms;
128 orbitalCubeFile = true;
129 }
130
131 if (fscanf(fp, "%d %lf %lf %lf", &n1, &elements[0], &elements[4], &elements[8]) != 4)
132 {
133 vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
134 << " Premature EOF while reading elements.");
135 fclose (fp);
136 return 0;
137 }
138 if (fscanf(fp, "%d %lf %lf %lf", &n2, &elements[1], &elements[5], &elements[9]) != 4)
139 {
140 vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
141 << " Premature EOF while reading elements.");
142 fclose (fp);
143 return 0;
144 }
145 if (fscanf(fp, "%d %lf %lf %lf", &n3, &elements[2], &elements[6], &elements[10]) != 4)
146 {
147 vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
148 << " Premature EOF while reading elements.");
149 fclose (fp);
150 return 0;
151 }
152 elements[12] = 0;
153 elements[13] = 0;
154 elements[14] = 0;
155 elements[15] = 1;
156
157 vtkDebugMacro(<< "Grid Size " << n1 << " " << n2 << " " << n3);
158
159 Transform->SetMatrix(elements);
160 Transform->Inverse();
161
162 this->ReadMolecule(fp, output);
163
164 if(orbitalCubeFile)
165 {
166 if (fscanf(fp,"%d", &numberOfOrbitals) != 1)
167 {
168 vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
169 << " Premature EOF while reading number of orbitals.");
170 fclose (fp);
171 return 0;
172 }
173 for(k = 0; k < numberOfOrbitals; k++)
174 {
175 if (fscanf(fp,"%f", &tmp) != 1)
176 {
177 vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
178 << " Premature EOF while reading orbitals.");
179 fclose (fp);
180 return 0;
181 }
182 }
183 }
184
185 vtkInformation *gridInfo = this->GetExecutive()->GetOutputInformation(1);
186 gridInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),
187 0, n1-1, 0, n2-1, 0, n3-1);
188 gridInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(),
189 gridInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT()),
190 6);
191 grid->SetExtent(
192 gridInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT()));
193
194 grid->SetOrigin(0, 0, 0);
195 grid->SetSpacing(1, 1, 1);
196 grid->AllocateScalars(VTK_FLOAT, 1);
197
198 grid->GetPointData()->GetScalars()->SetName(title);
199
200 cubedata = (float *)grid->GetPointData()->GetScalars()->GetVoidPointer(0);
201 N1N2 = n1*n2;
202
203 for(i = 0; i < n1; i++)
204 {
205 JN1 = 0;
206 for(j = 0; j < n2; j++)
207 {
208 for(k = 0; k < n3; k++)
209 {
210 if (fscanf(fp,"%f", &tmp) != 1)
211 {
212 vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
213 << " Premature EOF while reading scalars.");
214 fclose (fp);
215 return 0;
216 }
217 cubedata[k*N1N2 + JN1 + i] = tmp;
218 }
219 JN1 += n1;
220 }
221 }
222 fclose(fp);
223
224 return 1;
225 }
226
227 //----------------------------------------------------------------------------
ReadSpecificMolecule(FILE * fp)228 void vtkGaussianCubeReader::ReadSpecificMolecule(FILE* fp)
229 {
230 int i, j;
231 float x[3];
232 float dummy;
233
234 for(i = 0; i < this->NumberOfAtoms; i++)
235 {
236 if (fscanf(fp, "%d %f %f %f %f", &j, &dummy, x, x+1, x+2) != 5)
237 {
238 vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
239 << " Premature EOF while reading molecule.");
240 fclose (fp);
241 return;
242 }
243 this->Transform->TransformPoint(x, x);
244 this->Points->InsertNextPoint(x);
245 this->AtomType->InsertNextValue(j-1);
246 this->AtomTypeStrings->InsertNextValue("Xx");
247 this->Residue->InsertNextValue(-1);
248 this->Chain->InsertNextValue(0);
249 this->SecondaryStructures->InsertNextValue(0);
250 this->SecondaryStructuresBegin->InsertNextValue(0);
251 this->SecondaryStructuresEnd->InsertNextValue(0);
252 this->IsHetatm->InsertNextValue(0);
253 }
254 }
255
256 //----------------------------------------------------------------------------
GetGridOutput()257 vtkImageData *vtkGaussianCubeReader::GetGridOutput()
258 {
259 if (this->GetNumberOfOutputPorts() < 2)
260 {
261 return NULL;
262 }
263 return vtkImageData::SafeDownCast(
264 this->GetExecutive()->GetOutputData(1));
265 }
266
267 //----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)268 void vtkGaussianCubeReader::PrintSelf(ostream& os, vtkIndent indent)
269 {
270 this->Superclass::PrintSelf(os,indent);
271
272 os << "Filename: " << (this->FileName?this->FileName:"(none)") << "\n";
273
274 os << "Transform: ";
275 if( this->Transform )
276 {
277 os << endl;
278 this->Transform->PrintSelf(os, indent.GetNextIndent());
279 }
280 else
281 {
282 os << "(none)\n";
283 }
284 }
285
286 //----------------------------------------------------------------------------
287 // Default implementation - copy information from first input to all outputs
RequestInformation(vtkInformation * vtkNotUsed (request),vtkInformationVector ** vtkNotUsed (inputVector),vtkInformationVector * vtkNotUsed (outputVector))288 int vtkGaussianCubeReader::RequestInformation(
289 vtkInformation *vtkNotUsed(request),
290 vtkInformationVector **vtkNotUsed(inputVector),
291 vtkInformationVector *vtkNotUsed(outputVector))
292 {
293 // the set the information for the imagedat output
294 vtkInformation *gridInfo = this->GetExecutive()->GetOutputInformation(1);
295
296 FILE *fp;
297 char title[256];
298
299 if (!this->FileName)
300 {
301 return 0;
302 }
303
304 if ((fp = fopen(this->FileName, "r")) == NULL)
305 {
306 vtkErrorMacro(<< "File " << this->FileName << " not found");
307 return 0;
308 }
309
310 if (!fgets(title, 256, fp))
311 {
312 vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
313 << " Premature EOF while reading title.");
314 fclose (fp);
315 return 0;
316 }
317 if (!fgets(title, 256, fp))
318 {
319 vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
320 << " Premature EOF while reading title.");
321 fclose (fp);
322 return 0;
323 }
324
325 // Read in number of atoms, x-origin, y-origin z-origin
326 double tmpd;
327 int n1, n2, n3;
328 if (fscanf(fp, "%d %lf %lf %lf", &n1, &tmpd, &tmpd, &tmpd) != 4)
329 {
330 vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
331 << " Premature EOF while grid size.");
332 fclose (fp);
333 return 0;
334 }
335
336 if (fscanf(fp, "%d %lf %lf %lf", &n1, &tmpd, &tmpd, &tmpd) != 4)
337 {
338 vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
339 << " Premature EOF while grid size.");
340 fclose (fp);
341 return 0;
342 }
343 if (fscanf(fp, "%d %lf %lf %lf", &n2, &tmpd, &tmpd, &tmpd) != 4)
344 {
345 vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
346 << " Premature EOF while grid size.");
347 fclose (fp);
348 return 0;
349 }
350 if (fscanf(fp, "%d %lf %lf %lf", &n3, &tmpd, &tmpd, &tmpd) != 4)
351 {
352 vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
353 << " Premature EOF while grid size.");
354 fclose (fp);
355 return 0;
356 }
357
358 vtkDebugMacro(<< "Grid Size " << n1 << " " << n2 << " " << n3);
359 gridInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),
360 0, n1-1, 0, n2-1, 0, n3-1);
361 gridInfo->Set(vtkDataObject::ORIGIN(), 0, 0, 0);
362 gridInfo->Set(vtkDataObject::SPACING(), 1, 1, 1);
363
364 fclose(fp);
365
366 vtkDataObject::SetPointDataActiveScalarInfo(gridInfo, VTK_FLOAT, -1);
367 return 1;
368 }
369
370 //----------------------------------------------------------------------------
FillOutputPortInformation(int port,vtkInformation * info)371 int vtkGaussianCubeReader::FillOutputPortInformation(int port,
372 vtkInformation* info)
373 {
374 if(port == 0)
375 {
376 return this->Superclass::FillOutputPortInformation(port, info);
377 }
378 info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkImageData");
379 return 1;
380 }
381