vtkGaussianCubeReader.cxx 6.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 "vtkTransform.h"
27 28 29

#include <ctype.h>

30
vtkCxxRevisionMacro(vtkGaussianCubeReader, "1.15");
31 32 33 34 35 36 37
vtkStandardNewMacro(vtkGaussianCubeReader);

// Construct object with merging set to true.
vtkGaussianCubeReader::vtkGaussianCubeReader()
{
  this->FileName = NULL;
  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 51 52 53 54 55 56 57
}

vtkGaussianCubeReader::~vtkGaussianCubeReader()
{
  if (this->FileName)
    {
    delete [] this->FileName;
    }
  this->Transform->Delete();
  // must delete the second output added
}

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

67 68 69 70
  FILE *fp;
  char Title[256];
  char data_name[256];
  double elements[16];
71
  int JN1, N1N2, n1, n2, n3, i, j, k;
72 73 74 75 76
  float tmp, *Cube_data;

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

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

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

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

  fgets(Title, 256, fp);
  if(strtok(Title, ":") != NULL)
92 93 94
    {
    if(strtok(NULL, ":") != NULL)
      {
95 96
      strcpy(data_name, strtok(NULL, ":"));
      fprintf(stderr,"label = %s\n", data_name);
97 98
      }
    }
99 100
  fgets(Title, 256, fp);

101 102 103 104
  // Read in number of atoms, x-origin, y-origin z-origin
  //
  fscanf(fp, "%d %lf %lf %lf", &(this->NumberOfAtoms), &elements[3], 
         &elements[7], &elements[11]);
105
  if(this->NumberOfAtoms < 0 )
106
    {
107
    this->NumberOfAtoms = -this->NumberOfAtoms;
108
    }
109 110 111 112 113 114 115 116 117 118 119 120 121 122

  fscanf(fp, "%d %lf %lf %lf", &n1, &elements[0], &elements[4], &elements[8]);
  fscanf(fp, "%d %lf %lf %lf", &n2, &elements[1], &elements[5], &elements[9]);
  fscanf(fp, "%d %lf %lf %lf", &n3, &elements[2], &elements[6], &elements[10]);
  elements[12] = 0;
  elements[13] = 0;
  elements[14] = 0;
  elements[15] = 1;
  
  vtkDebugMacro(<< "Grid Size " << n1 << " " << n2 << " " << n3);

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

123
  this->ReadMolecule(fp, output);
124

125 126 127 128 129 130 131 132
  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()));
133 134 135 136 137 138 139 140 141 142 143

  grid->SetOrigin(0, 0, 0);
  grid->SetSpacing(1, 1, 1);
  grid->SetScalarTypeToFloat();
  grid->AllocateScalars();

  grid->GetPointData()->GetScalars()->SetName("Gaussian Cube density");

  Cube_data = (float *)grid->GetPointData()->GetScalars()->GetVoidPointer(0);
  N1N2 = n1*n2;

144 145
  for(i = 0; i < n1; i++) 
    {
146
    JN1 = 0;
147 148 149 150 151 152 153
    for(j = 0; j < n2; j++) 
      {
      for(k = 0; k < n3; k++) 
        {
        fscanf(fp,"%f", &tmp);
        Cube_data[k*N1N2 + JN1 + i] = tmp;
        }
154
      JN1 += n1;
155
      }
156 157
    }
  fclose(fp);
158 159

  return 1;
160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177
}

void vtkGaussianCubeReader::ReadSpecificMolecule(FILE* fp)
{
  int i, j;
  float x[3];
  float dummy;

  for(i = 0; i < this->NumberOfAtoms; i++) 
    {
    fscanf(fp, "%d %f %f %f %f", &j, &dummy, x, x+1, x+2);
    this->Points->InsertNextPoint(x);
    this->AtomType->InsertNextValue(j-1);
    }
}

vtkImageData *vtkGaussianCubeReader::GetGridOutput()
{
178 179 180 181 182 183
  if (this->GetNumberOfOutputPorts() < 2)
    {
    return NULL;
    }
  return vtkImageData::SafeDownCast(
    this->GetExecutive()->GetOutputData(1));
184 185 186 187 188 189
}

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

190 191 192 193 194 195
  os << "Filename: " << (this->FileName?this->FileName:"(none)") << "\n";
  
  os << "Tranform: ";
  if( this->Transform )
    {
    os << endl;
Mathieu Malaterre's avatar
doh  
Mathieu Malaterre committed
196
    this->Transform->PrintSelf(os, indent.GetNextIndent());
197 198 199 200 201
    }
  else
    {
    os << "(none)\n";
    }
202
}
203

204
// Default implementation - copy information from first input to all outputs
205 206 207 208
int vtkGaussianCubeReader::RequestInformation(
  vtkInformation *vtkNotUsed(request),
  vtkInformationVector **vtkNotUsed(inputVector),
  vtkInformationVector *vtkNotUsed(outputVector))
209 210
{
  // the set the information for the imagedat output
211 212
  vtkInformation *gridInfo = this->GetExecutive()->GetOutputInformation(1);

213 214 215 216 217
  FILE *fp;
  char Title[256];
  
  if (!this->FileName)
    {
218
    return 0;
219 220 221 222 223
    }
  
  if ((fp = fopen(this->FileName, "r")) == NULL)
    {
    vtkErrorMacro(<< "File " << this->FileName << " not found");
224
    return 0;
225 226 227 228 229 230 231 232 233 234 235 236 237 238 239
    }
  
  fgets(Title, 256, fp);
  fgets(Title, 256, fp);

  // Read in number of atoms, x-origin, y-origin z-origin
  double tmpd;
  int n1, n2, n3;
  fscanf(fp, "%d %lf %lf %lf", &n1, &tmpd, &tmpd, &tmpd);
  
  fscanf(fp, "%d %lf %lf %lf", &n1, &tmpd, &tmpd, &tmpd);
  fscanf(fp, "%d %lf %lf %lf", &n2, &tmpd, &tmpd, &tmpd);
  fscanf(fp, "%d %lf %lf %lf", &n3, &tmpd, &tmpd, &tmpd);
  
  vtkDebugMacro(<< "Grid Size " << n1 << " " << n2 << " " << n3);
240 241 242 243
  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);
244 245

  fclose(fp);
246

247
  vtkDataObject::SetPointDataActiveScalarInfo(gridInfo, VTK_FLOAT, -1);
248
  return 1;
249 250
}

251 252 253 254 255 256 257 258
//----------------------------------------------------------------------------
int vtkGaussianCubeReader::FillOutputPortInformation(int port,
                                                     vtkInformation* info)
{
  if(port == 0)
    {
    return this->Superclass::FillOutputPortInformation(port, info);
    }
259
  info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkImageData");
260 261
  return 1;
}