vtkGaussianCubeReader.cxx 10.9 KB
Newer Older
1 2 3 4 5
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkGaussianCubeReader.cxx

Brad King's avatar
Brad King committed
6 7 8
  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9

Brad King's avatar
Brad King committed
10 11 12
     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.
13 14 15

=========================================================================*/
#include "vtkGaussianCubeReader.h"
16

17
#include "vtkImageData.h"
18
#include "vtkInformation.h"
19
#include "vtkInformationVector.h"
20 21 22 23 24
#include "vtkPointData.h"
#include "vtkPoints.h"
#include "vtkFloatArray.h"
#include "vtkIdTypeArray.h"
#include "vtkObjectFactory.h"
25
#include "vtkStreamingDemandDrivenPipeline.h"
26
#include "vtkStringArray.h"
27
#include "vtkTransform.h"
28

Sean McBride's avatar
Sean McBride committed
29
#include <cctype>
30 31 32

vtkStandardNewMacro(vtkGaussianCubeReader);

33
//----------------------------------------------------------------------------
34 35 36 37
// Construct object with merging set to true.
vtkGaussianCubeReader::vtkGaussianCubeReader()
{
  this->Transform = vtkTransform::New();
38
  // Add the second output for the grid data
39 40

  this->SetNumberOfOutputPorts(2);
41 42 43
  vtkImageData *grid;
  grid = vtkImageData::New();
  grid->ReleaseData();
44
  this->GetExecutive()->SetOutputData(1, grid);
45
  grid->Delete();
46 47
}

48
//----------------------------------------------------------------------------
49 50
vtkGaussianCubeReader::~vtkGaussianCubeReader()
{
51

52 53 54 55
  this->Transform->Delete();
  // must delete the second output added
}

56
//----------------------------------------------------------------------------
57 58 59 60
int vtkGaussianCubeReader::RequestData(
  vtkInformation *vtkNotUsed(request),
  vtkInformationVector **vtkNotUsed(inputVector),
  vtkInformationVector *outputVector)
61
{
62 63 64 65
  vtkInformation *outInfo = outputVector->GetInformationObject(0);
  vtkPolyData *output = vtkPolyData::SafeDownCast(
    outInfo->Get(vtkDataObject::DATA_OBJECT()));

66
  FILE *fp;
67
  char title[256];
68 69
  char data_name[256];
  double elements[16];
70
  int JN1, N1N2, n1, n2, n3, i, j, k;
71 72 73
  float tmp, *cubedata;
  bool orbitalCubeFile = false;
  int numberOfOrbitals;
74 75 76 77

  // Output 0 (the default is the polydata)
  // Output 1 will be the gridded Image data

78 79
  vtkImageData *grid = this->GetGridOutput();

80
  if (!this->FileName)
81
  {
82
    return 0;
83
  }
84

85
  if ((fp = fopen(this->FileName, "r")) == nullptr)
86
  {
87
    vtkErrorMacro(<< "File " << this->FileName << " not found");
88
    return 0;
89
  }
90

Bill Lorensen's avatar
Bill Lorensen committed
91
  if (!fgets(title, 256, fp))
92
  {
Bill Lorensen's avatar
Bill Lorensen committed
93 94 95 96
    vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
                   << " Premature EOF while reading title.");
    fclose (fp);
    return 0;
97
  }
98
  if(strtok(title, ":") != nullptr)
99
  {
100
    if(strtok(nullptr, ":") != nullptr)
101
    {
102
      strcpy(data_name, strtok(nullptr, ":"));
103
      fprintf(stderr,"label = %s\n", data_name);
104
    }
105
  }
Bill Lorensen's avatar
Bill Lorensen committed
106
  if (!fgets(title, 256, fp))
107
  {
Bill Lorensen's avatar
Bill Lorensen committed
108 109 110 111
    vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
                   << " Premature EOF while reading title.");
    fclose (fp);
    return 0;
112
  }
113

114 115
  // Read in number of atoms, x-origin, y-origin z-origin
  //
Bill Lorensen's avatar
Bill Lorensen committed
116
  if (fscanf(fp, "%d %lf %lf %lf", &(this->NumberOfAtoms), &elements[3], &elements[7], &elements[11]) != 4)
117
  {
Bill Lorensen's avatar
Bill Lorensen committed
118 119 120 121
    vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
                   << " Premature EOF while reading atoms, x-origin y-origin z-origin.");
    fclose (fp);
    return 0;
122
  }
123
  if(this->NumberOfAtoms < 0 )
124
  {
125
    this->NumberOfAtoms = -this->NumberOfAtoms;
126
    orbitalCubeFile = true;
127
  }
128

Bill Lorensen's avatar
Bill Lorensen committed
129
  if (fscanf(fp, "%d %lf %lf %lf", &n1, &elements[0], &elements[4], &elements[8]) != 4)
130
  {
Bill Lorensen's avatar
Bill Lorensen committed
131 132 133 134
    vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
                   << " Premature EOF while reading elements.");
    fclose (fp);
    return 0;
135
  }
Bill Lorensen's avatar
Bill Lorensen committed
136
  if (fscanf(fp, "%d %lf %lf %lf", &n2, &elements[1], &elements[5], &elements[9]) != 4)
137
  {
Bill Lorensen's avatar
Bill Lorensen committed
138 139 140 141
    vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
                   << " Premature EOF while reading elements.");
    fclose (fp);
    return 0;
142
  }
Bill Lorensen's avatar
Bill Lorensen committed
143
  if (fscanf(fp, "%d %lf %lf %lf", &n3, &elements[2], &elements[6], &elements[10]) != 4)
144
  {
Bill Lorensen's avatar
Bill Lorensen committed
145 146 147 148
    vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
                   << " Premature EOF while reading elements.");
    fclose (fp);
    return 0;
149
  }
150 151 152 153
  elements[12] = 0;
  elements[13] = 0;
  elements[14] = 0;
  elements[15] = 1;
154

155 156 157 158 159
  vtkDebugMacro(<< "Grid Size " << n1 << " " << n2 << " " << n3);

  Transform->SetMatrix(elements);
  Transform->Inverse();

160
  this->ReadMolecule(fp, output);
161

162
  if(orbitalCubeFile)
163
  {
Bill Lorensen's avatar
Bill Lorensen committed
164
    if (fscanf(fp,"%d", &numberOfOrbitals) != 1)
165
    {
Bill Lorensen's avatar
Bill Lorensen committed
166 167 168 169
      vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
                     << " Premature EOF while reading number of orbitals.");
      fclose (fp);
      return 0;
170
    }
171
    for(k = 0; k < numberOfOrbitals; k++)
172
    {
Bill Lorensen's avatar
Bill Lorensen committed
173
      if (fscanf(fp,"%f", &tmp) != 1)
174
      {
Bill Lorensen's avatar
Bill Lorensen committed
175 176 177 178
        vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
                       << " Premature EOF while reading orbitals.");
        fclose (fp);
        return 0;
179 180
      }
    }
181
  }
182

183 184 185 186 187 188 189 190
  vtkInformation *gridInfo = this->GetExecutive()->GetOutputInformation(1);
  gridInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),
                0, n1-1, 0, n2-1, 0, n3-1);
  gridInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(),
                gridInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT()),
                6);
  grid->SetExtent(
    gridInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT()));
191 192 193

  grid->SetOrigin(0, 0, 0);
  grid->SetSpacing(1, 1, 1);
194
  grid->AllocateScalars(VTK_FLOAT, 1);
195

196
  grid->GetPointData()->GetScalars()->SetName(title);
197

198
  cubedata = (float *)grid->GetPointData()->GetScalars()->GetVoidPointer(0);
199 200
  N1N2 = n1*n2;

201
  for(i = 0; i < n1; i++)
202
  {
203
    JN1 = 0;
204
    for(j = 0; j < n2; j++)
205
    {
206
      for(k = 0; k < n3; k++)
207
      {
Bill Lorensen's avatar
Bill Lorensen committed
208
        if (fscanf(fp,"%f", &tmp) != 1)
209
        {
Bill Lorensen's avatar
Bill Lorensen committed
210 211 212 213
          vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
                         << " Premature EOF while reading scalars.");
          fclose (fp);
          return 0;
214
        }
215
        cubedata[k*N1N2 + JN1 + i] = tmp;
216
      }
217
      JN1 += n1;
218
    }
219
  }
220
  fclose(fp);
221 222

  return 1;
223 224
}

225
//----------------------------------------------------------------------------
226 227 228 229 230 231
void vtkGaussianCubeReader::ReadSpecificMolecule(FILE* fp)
{
  int i, j;
  float x[3];
  float dummy;

232
  for(i = 0; i < this->NumberOfAtoms; i++)
233
  {
Bill Lorensen's avatar
Bill Lorensen committed
234
    if (fscanf(fp, "%d %f %f %f %f", &j, &dummy, x, x+1, x+2) != 5)
235
    {
Bill Lorensen's avatar
Bill Lorensen committed
236 237 238 239
      vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
                     << " Premature EOF while reading molecule.");
      fclose (fp);
      return;
240
    }
241
    this->Transform->TransformPoint(x, x);
242 243
    this->Points->InsertNextPoint(x);
    this->AtomType->InsertNextValue(j-1);
244 245 246 247 248 249 250
    this->AtomTypeStrings->InsertNextValue("Xx");
    this->Residue->InsertNextValue(-1);
    this->Chain->InsertNextValue(0);
    this->SecondaryStructures->InsertNextValue(0);
    this->SecondaryStructuresBegin->InsertNextValue(0);
    this->SecondaryStructuresEnd->InsertNextValue(0);
    this->IsHetatm->InsertNextValue(0);
251
  }
252 253
}

254
//----------------------------------------------------------------------------
255 256
vtkImageData *vtkGaussianCubeReader::GetGridOutput()
{
257
  if (this->GetNumberOfOutputPorts() < 2)
258
  {
259
    return nullptr;
260
  }
261 262
  return vtkImageData::SafeDownCast(
    this->GetExecutive()->GetOutputData(1));
263 264
}

265
//----------------------------------------------------------------------------
266 267 268 269
void vtkGaussianCubeReader::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os,indent);

270
  os << "Filename: " << (this->FileName?this->FileName:"(none)") << "\n";
271

Kyle Lutz's avatar
Kyle Lutz committed
272
  os << "Transform: ";
273
  if( this->Transform )
274
  {
275
    os << endl;
Mathieu Malaterre's avatar
doh  
Mathieu Malaterre committed
276
    this->Transform->PrintSelf(os, indent.GetNextIndent());
277
  }
278
  else
279
  {
280
    os << "(none)\n";
281
  }
282
}
283

284
//----------------------------------------------------------------------------
285
// Default implementation - copy information from first input to all outputs
286 287 288 289
int vtkGaussianCubeReader::RequestInformation(
  vtkInformation *vtkNotUsed(request),
  vtkInformationVector **vtkNotUsed(inputVector),
  vtkInformationVector *vtkNotUsed(outputVector))
290 291
{
  // the set the information for the imagedat output
292 293
  vtkInformation *gridInfo = this->GetExecutive()->GetOutputInformation(1);

294
  FILE *fp;
295
  char title[256];
296

297
  if (!this->FileName)
298
  {
299
    return 0;
300
  }
301

302
  if ((fp = fopen(this->FileName, "r")) == nullptr)
303
  {
304
    vtkErrorMacro(<< "File " << this->FileName << " not found");
305
    return 0;
306
  }
307

Bill Lorensen's avatar
Bill Lorensen committed
308
  if (!fgets(title, 256, fp))
309
  {
Bill Lorensen's avatar
Bill Lorensen committed
310 311 312 313
    vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
                   << " Premature EOF while reading title.");
    fclose (fp);
    return 0;
314
  }
Bill Lorensen's avatar
Bill Lorensen committed
315
  if (!fgets(title, 256, fp))
316
  {
Bill Lorensen's avatar
Bill Lorensen committed
317 318 319 320
    vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
                   << " Premature EOF while reading title.");
    fclose (fp);
    return 0;
321
  }
322 323 324 325

  // Read in number of atoms, x-origin, y-origin z-origin
  double tmpd;
  int n1, n2, n3;
Bill Lorensen's avatar
Bill Lorensen committed
326
  if (fscanf(fp, "%d %lf %lf %lf", &n1, &tmpd, &tmpd, &tmpd) != 4)
327
  {
Bill Lorensen's avatar
Bill Lorensen committed
328 329 330 331
    vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
                   << " Premature EOF while grid size.");
    fclose (fp);
    return 0;
332
  }
333

Bill Lorensen's avatar
Bill Lorensen committed
334
  if (fscanf(fp, "%d %lf %lf %lf", &n1, &tmpd, &tmpd, &tmpd) != 4)
335
  {
Bill Lorensen's avatar
Bill Lorensen committed
336 337 338 339
    vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
                   << " Premature EOF while grid size.");
    fclose (fp);
    return 0;
340
  }
Bill Lorensen's avatar
Bill Lorensen committed
341
  if (fscanf(fp, "%d %lf %lf %lf", &n2, &tmpd, &tmpd, &tmpd) != 4)
342
  {
Bill Lorensen's avatar
Bill Lorensen committed
343 344 345 346
    vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
                   << " Premature EOF while grid size.");
    fclose (fp);
    return 0;
347
  }
Bill Lorensen's avatar
Bill Lorensen committed
348
  if (fscanf(fp, "%d %lf %lf %lf", &n3, &tmpd, &tmpd, &tmpd) != 4)
349
  {
Bill Lorensen's avatar
Bill Lorensen committed
350 351 352 353
    vtkErrorMacro ("GaussianCubeReader error reading file: " << this->FileName
                   << " Premature EOF while grid size.");
    fclose (fp);
    return 0;
354
  }
355

356
  vtkDebugMacro(<< "Grid Size " << n1 << " " << n2 << " " << n3);
357 358 359 360
  gridInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),
                0, n1-1, 0, n2-1, 0, n3-1);
  gridInfo->Set(vtkDataObject::ORIGIN(), 0, 0, 0);
  gridInfo->Set(vtkDataObject::SPACING(), 1, 1, 1);
361 362

  fclose(fp);
363

364
  vtkDataObject::SetPointDataActiveScalarInfo(gridInfo, VTK_FLOAT, -1);
365
  return 1;
366 367
}

368 369 370 371 372
//----------------------------------------------------------------------------
int vtkGaussianCubeReader::FillOutputPortInformation(int port,
                                                     vtkInformation* info)
{
  if(port == 0)
373
  {
374
    return this->Superclass::FillOutputPortInformation(port, info);
375
  }
376
  info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkImageData");
377 378
  return 1;
}