source: nanovis/trunk/VtkDataSetReader.cpp @ 5704

/* -*- mode: c++; c-basic-offset: 4; indent-tabs-mode: nil -*- */
/*
 * Copyright (C) 2004-2013  HUBzero Foundation, LLC
 *
 */

#include <iostream>
#include <float.h>

#include <vtkSmartPointer.h>
#include <vtkDataSet.h>
#include <vtkImageData.h>
#include <vtkPointData.h>
#include <vtkCellData.h>
#include <vtkDataArray.h>
#include <vtkDataSetReader.h>
#include <vtkCharArray.h>

#include <vrmath/Vector3f.h>

#include "ReaderCommon.h"
#include "DataSetResample.h"

#include "config.h"
#include "nanovis.h"
#include "VtkDataSetReader.h"
#include "Volume.h"
#include "Trace.h"

using namespace nv;

static void print(vtkDataSet *ds)
{
    if (ds == NULL)
        return;

    TRACE("DataSet class: %s", ds->GetClassName());

    TRACE("DataSet memory: %g MiB", ds->GetActualMemorySize()/1024.);

    double bounds[6];
    ds->GetBounds(bounds);

    // Topology
    TRACE("DataSet bounds: %g %g %g %g %g %g",
          bounds[0], bounds[1],
          bounds[2], bounds[3],
          bounds[4], bounds[5]);
    TRACE("Points: %d Cells: %d", ds->GetNumberOfPoints(), ds->GetNumberOfCells());

    double dataRange[2];
    if (ds->GetPointData() != NULL) {
        TRACE("PointData arrays: %d", ds->GetPointData()->GetNumberOfArrays());
        for (int i = 0; i < ds->GetPointData()->GetNumberOfArrays(); i++) {
            if (ds->GetPointData()->GetArray(i) != NULL) {
                ds->GetPointData()->GetArray(i)->GetRange(dataRange, -1);
                TRACE("PointData[%d]: '%s' comp:%d, (%g,%g)", i,
                      ds->GetPointData()->GetArrayName(i),
                      ds->GetPointData()->GetAbstractArray(i)->GetNumberOfComponents(),
                      dataRange[0], dataRange[1]);
            } else {
                TRACE("PointData[%d]: '%s' comp:%d", i,
                      ds->GetPointData()->GetArrayName(i),
                      ds->GetPointData()->GetAbstractArray(i)->GetNumberOfComponents());
            }
        }
        if (ds->GetPointData()->GetScalars() != NULL) {
            TRACE("Active point scalars: %s", ds->GetPointData()->GetScalars()->GetName());
        }
        if (ds->GetPointData()->GetVectors() != NULL) {
            TRACE("Active point vectors: %s", ds->GetPointData()->GetVectors()->GetName());
        }
    }
    if (ds->GetCellData() != NULL) {
        TRACE("CellData arrays: %d", ds->GetCellData()->GetNumberOfArrays());
        for (int i = 0; i < ds->GetCellData()->GetNumberOfArrays(); i++) {
            if (ds->GetCellData()->GetArray(i) != NULL) {
                ds->GetCellData()->GetArray(i)->GetRange(dataRange, -1);
                TRACE("CellData[%d]: '%s' comp:%d, (%g,%g)", i,
                      ds->GetCellData()->GetArrayName(i),
                      ds->GetCellData()->GetAbstractArray(i)->GetNumberOfComponents(),
                      dataRange[0], dataRange[1]);
            } else {
                TRACE("CellData[%d]: '%s' comp:%d", i,
                      ds->GetCellData()->GetArrayName(i),
                      ds->GetCellData()->GetAbstractArray(i)->GetNumberOfComponents());
            }
        }
        if (ds->GetCellData()->GetScalars() != NULL) {
            TRACE("Active cell scalars: %s", ds->GetCellData()->GetScalars()->GetName());
        }
        if (ds->GetCellData()->GetVectors() != NULL) {
            TRACE("Active cell vectors: %s", ds->GetCellData()->GetVectors()->GetName());
        }
    }
    if (ds->GetFieldData() != NULL) {
        TRACE("FieldData arrays: %d", ds->GetFieldData()->GetNumberOfArrays());
        for (int i = 0; i < ds->GetFieldData()->GetNumberOfArrays(); i++) {
            if (ds->GetFieldData()->GetArray(i) != NULL) {
                ds->GetFieldData()->GetArray(i)->GetRange(dataRange, -1);
                TRACE("FieldData[%d]: '%s' comp:%d, tuples:%d (%g,%g)", i,
                      ds->GetFieldData()->GetArrayName(i),
                      ds->GetFieldData()->GetAbstractArray(i)->GetNumberOfComponents(),
                      ds->GetFieldData()->GetAbstractArray(i)->GetNumberOfTuples(),
                      dataRange[0], dataRange[1]);
            } else {
                TRACE("FieldData[%d]: '%s' comp:%d, tuples:%d", i,
                      ds->GetFieldData()->GetArrayName(i),
                      ds->GetFieldData()->GetAbstractArray(i)->GetNumberOfComponents(),
                      ds->GetFieldData()->GetAbstractArray(i)->GetNumberOfTuples());
            }
        }
    }
}

Volume *
nv::load_vtk_volume_stream(const char *tag, const char *bytes, int nBytes)
{
    TRACE("Enter tag:%s", tag);

    vtkSmartPointer<vtkImageData> resampledDataSet;
    {
        vtkSmartPointer<vtkDataSet> dataSet;
        {
            vtkSmartPointer<vtkDataSetReader> reader = vtkSmartPointer<vtkDataSetReader>::New();
            vtkSmartPointer<vtkCharArray> dataSetString = vtkSmartPointer<vtkCharArray>::New();

            dataSetString->SetArray(const_cast<char *>(bytes), nBytes, 1);
            reader->SetInputArray(dataSetString);
            reader->ReadFromInputStringOn();
            //reader->ReadAllNormalsOn();
            //reader->ReadAllTCoordsOn();
            reader->ReadAllScalarsOn();
            //reader->ReadAllColorScalarsOn();
            reader->ReadAllVectorsOn();
            //reader->ReadAllTensorsOn();
            reader->ReadAllFieldsOn();
            reader->SetLookupTableName("");
            reader->Update();

            dataSet = reader->GetOutput();
        }

        if (dataSet == NULL)
            return NULL;

        print(dataSet);

        int maxDim = 64;

        resampledDataSet = vtkImageData::SafeDownCast(dataSet.GetPointer());
        if (resampledDataSet != NULL) {
            // Have a uniform grid, check if we need to resample
#ifdef DOWNSAMPLE_DATA
            if (resampledDataSet->GetDimensions()[0] > maxDim ||
                resampledDataSet->GetDimensions()[1] > maxDim ||
                resampledDataSet->GetDimensions()[2] > maxDim) {
                resampledDataSet = resampleVTKDataSet(dataSet, maxDim);
            }
#endif
        } else {
            resampledDataSet = resampleVTKDataSet(dataSet, maxDim);
        }
    }

    int nx, ny, nz, npts;
    // origin
    double x0, y0, z0;
    // deltas (cell dimensions)
    double dx, dy, dz;
    // length of volume edges
    double lx, ly, lz;

    nx = resampledDataSet->GetDimensions()[0];
    ny = resampledDataSet->GetDimensions()[1];
    nz = resampledDataSet->GetDimensions()[2];
    npts = nx * ny * nz;
    resampledDataSet->GetSpacing(dx, dy, dz);
    resampledDataSet->GetOrigin(x0, y0, z0);

    lx = (nx - 1) * dx;
    ly = (ny - 1) * dy;
    lz = (nz - 1) * dz;

    Volume *volume = NULL;
    double vmin = DBL_MAX;
    double nzero_min = DBL_MAX;
    double vmax = -DBL_MAX;

    float *data = new float[npts * 4];
    memset(data, 0, npts * 4);

    bool isVectorData = (resampledDataSet->GetPointData()->GetVectors() != NULL);

    int ix = 0;
    int iy = 0;
    int iz = 0;
    vtkDataArray *mask = resampledDataSet->GetPointData()->GetArray("vtkValidPointMask");
    int numValidPoints = 0;
    for (int p = 0; p < npts; p++) {
        int nindex = p * 4;
        double val;
        int loc[3];
        loc[0] = ix; loc[1] = iy; loc[2] = iz;
        vtkIdType idx = resampledDataSet->ComputePointId(loc);
        bool valid = (mask == NULL) ? true : (mask->GetComponent(idx, 0) != 0.0);
        if (isVectorData) {
            double vec[3];
            resampledDataSet->GetPointData()->GetVectors()->GetTuple(idx, vec);
            val = sqrt(vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2]);
            data[nindex] = (float)val;
            data[nindex+1] = (float)vec[0];
            data[nindex+2] = (float)vec[1];
            data[nindex+3] = (float)vec[2];
        } else {
            //val = resampledDataSet->GetScalarComponentAsDouble(ix, iy, iz, 0);
            val = resampledDataSet->GetPointData()->GetScalars()->GetComponent(idx, 0);
            data[nindex] = valid ? (float)val : -FLT_MAX;
        }
        if (valid) {
            numValidPoints++;
            if (val < vmin) {
                vmin = val;
            } else if (val > vmax) {
                vmax = val;
            }
            if (val != 0.0 && val < nzero_min) {
                nzero_min = val;
            }
        }

        if (++ix >= nx) {
            ix = 0;
            if (++iy >= ny) {
                iy = 0;
                ++iz;
            }
        }
    }

    if (numValidPoints < 1) {
        ERROR("No valid points in resample volume");
        vmin = 0.0;
        vmax = 1.0;
    } else {
        TRACE("Valid points: %d of %d", numValidPoints, resampledDataSet->GetNumberOfPoints());
    }

    if (isVectorData) {
        // Normalize magnitude [0,1] and vector components [0,1]
        normalizeVector(data, npts, vmin, vmax);
    } else {
        // Normalize all values [0,1], -1 => out of bounds
        normalizeScalar(data, npts, 4, vmin, vmax);
        computeSimpleGradient(data, nx, ny, nz,
                              dx, dy, dz);
    }

    TRACE("isVectorData: %s", isVectorData ? "yes" : "no");
    TRACE("nx = %i ny = %i nz = %i", nx, ny, nz);
    TRACE("x0 = %g y0 = %g z0 = %g", x0, y0, z0);
    TRACE("lx = %g ly = %g lz = %g", lx, ly, lz);
    TRACE("dx = %g dy = %g dz = %g", dx, dy, dz);
    TRACE("dataMin = %g dataMax = %g nzero_min = %g",
          vmin, vmax, nzero_min);

    volume = NanoVis::loadVolume(tag, nx, ny, nz, 4, data,
                                 vmin, vmax, nzero_min);
    volume->xAxis.setRange(x0, x0 + lx);
    volume->yAxis.setRange(y0, y0 + ly);
    volume->zAxis.setRange(z0, z0 + lz);
    volume->updatePending = true;

    delete [] data;

    TRACE("Leave");
    return volume;
}