yade-users team mailing list archive

[William Chèvremont <question681105@xxxxxxxxxxxxxxxxxxxxx>]

Question #681105 on Yade changed:
https://answers.launchpad.net/yade/+question/681105

William Chèvremont proposed the following answer:
You can read vtk data from matlab, using this function. You need first
to compile it using mex:

/* Convert four uint8 to float64 value */

#include "mex.h"
#include "matrix.h"
#include <iostream>
#include <string>
#include <vtkSmartPointer.h>
#include <vtkXMLReader.h>
#include <vtkXMLUnstructuredGridReader.h>
#include <vtkXMLPolyDataReader.h>
#include <vtkXMLStructuredGridReader.h>
#include <vtkXMLRectilinearGridReader.h>
#include <vtkXMLHyperOctreeReader.h>
#include <vtkXMLCompositeDataReader.h>
#include <vtkXMLStructuredGridReader.h>
#include <vtkXMLImageDataReader.h>
#include <vtkDataSetReader.h>
#include <vtkDataSet.h>
#include <vtkUnstructuredGrid.h>
#include <vtkRectilinearGrid.h>
#include <vtkHyperOctree.h>
#include <vtkImageData.h>
#include <vtkPolyData.h>
#include <vtkStructuredGrid.h>
#include <vtkPointData.h>
#include <vtkCellData.h>
#include <vtkFieldData.h>
#include <vtkCellTypes.h>
#include <vtksys/SystemTools.hxx>
 
#include <map>

mxArray* readData(const char*);

/* The gateway function */
void mexFunction(int nlhs, mxArray *plhs[],
                 int nrhs, const mxArray *prhs[])
{
/* variable declarations here */
    if(nrhs != 1)
    {
        mexErrMsgIdAndTxt("VtkUncompresserToolBox:TooMuchInput","One input required");
    }
    
    if(!mxIsChar(prhs[0]))
    {
        mexErrMsgIdAndTxt("VtkUncompresserToolBox:BadInput","Input must be a filename");
    }
    
    int buflen = mxGetNumberOfElements(prhs[0]) + 1;
    char* buf = (char*)mxCalloc(buflen, sizeof(char));

    /* Copy the string data from string_array_ptr and place it into buf. */ 
    if (mxGetString(prhs[0], buf, buflen) != 0)
    mexErrMsgIdAndTxt( "VtkUncompresserToolBox:invalidStringArray",
    "Could not convert string data.");
    
    plhs[0] = readData(buf);

/* code here */
}
//
// DumpXMLFile - report on the contents of an XML or legacy vtk file
//  Usage: DumpXMLFile XMLFile1 XMLFile2 ...
//         where
//         XMLFile is a vtk XML file of type .vtu, .vtp, .vts, .vtr,
//         .vti, .vto 
//

 
template<class TReader> vtkDataSet *ReadAnXMLFile(const char*fileName)
{
  vtkSmartPointer<TReader> reader =
    vtkSmartPointer<TReader>::New();
  reader->SetFileName(fileName);
  reader->Update();
  reader->GetOutput()->Register(reader);
  return vtkDataSet::SafeDownCast(reader->GetOutput());
}
 
mxArray* readData(const char* filename)
{
    
  // Process each file on the command line

    vtkDataSet *dataSet;
    std::string extension =
      vtksys::SystemTools::GetFilenameLastExtension(filename);
    // Dispatch based on the file extension
    if (extension == ".vtu")
    {
        dataSet = ReadAnXMLFile<vtkXMLUnstructuredGridReader> (filename);
    }
    else if (extension == ".vtp")
    {
        dataSet = ReadAnXMLFile<vtkXMLPolyDataReader> (filename);
    }
    else if (extension == ".vts")
    {
        dataSet = ReadAnXMLFile<vtkXMLStructuredGridReader> (filename);
    }
    else if (extension == ".vtr")
    {
        dataSet = ReadAnXMLFile<vtkXMLRectilinearGridReader> (filename);
    }
    else if (extension == ".vti")
    {
        dataSet = ReadAnXMLFile<vtkXMLImageDataReader> (filename);
    }
    else if (extension == ".vto")
    {
        dataSet = ReadAnXMLFile<vtkXMLHyperOctreeReader> (filename);
    }
    else if (extension == ".vtk")
    {
        dataSet = ReadAnXMLFile<vtkDataSetReader> (filename);
    }
    else
    {
          mexErrMsgIdAndTxt("VtkUncompresserToolBox:UnknownExt","Invalid extension");
          return 0;
    }
    
    const char *fieldnames[5] = {"Points","CellCount","PointData","CellData","FieldData"};
 
    
    int numberOfCells = dataSet->GetNumberOfCells();    
    int numberOfPoints = dataSet->GetNumberOfPoints();
    
    mxArray* data = mxCreateStructMatrix(1,1,5,fieldnames);
    
    mxArray* points = mxCreateDoubleMatrix(numberOfPoints,3,mxREAL);
    double* pointsPtr = mxGetPr(points);
    double pos[3];
    
    for(int i=0;i<numberOfPoints;i++)
    {
        dataSet->GetPoint(i,pos);
        pointsPtr[i] = pos[0];
        pointsPtr[i+numberOfPoints] = pos[1];
        pointsPtr[i+numberOfPoints*2] = pos[2];
    }
    
    mxSetFieldByNumber(data,0,0,points);
    
    //mxSetFieldByNumber(data,0,0,mxCreateDoubleScalar(numberOfPoints));
    mxSetFieldByNumber(data,0,1,mxCreateDoubleScalar(numberOfCells));    
    
    mxArray *PD = mxCreateStructMatrix(1,1,0,0);
    vtkPointData *pd = dataSet->GetPointData();
    if (pd)
    {
        for (int i = 0; i < pd->GetNumberOfArrays(); i++)
        {
            int id = mxAddField(PD,pd->GetArrayName(i));
            
            vtkAbstractArray* d = pd->GetAbstractArray(i);
            
            if(d && d->GetDataType() == VTK_DOUBLE)
            {
                mxArray* dat = mxCreateDoubleMatrix(d->GetNumberOfTuples(),d->GetNumberOfComponents(),mxREAL);
                double *datPtr = mxGetPr(dat);
                double *vtkDat = (double*)d->GetVoidPointer(0);
                
                for(int j=0;j<d->GetNumberOfComponents();j++)
                {
                    for(int k=0;k<d->GetNumberOfTuples();k++)
                    {
                        datPtr[k+j*d->GetNumberOfTuples()] = vtkDat[j+k*d->GetNumberOfComponents()];
                    }
                }
                
                mxSetFieldByNumber(PD,0,id,dat);
               // mxSetFieldByNumber(PD,0,id,mxCreateDoubleScalar(d->GetDataType()));
            }
        }
    }
    mxSetFieldByNumber(data,0,2,PD);
    
    // Now check for cell data
    mxArray *CD = mxCreateStructMatrix(1,1,0,0);
    vtkCellData *cd = dataSet->GetCellData();
    if (cd)
    {
        for (int i = 0; i < cd->GetNumberOfArrays(); i++)
        {
            int id = mxAddField(CD,cd->GetArrayName(i));  
            vtkAbstractArray* d = cd->GetAbstractArray(i);
            
            if(d && d->GetDataType() == VTK_DOUBLE)
            {
                mxArray* dat = mxCreateDoubleMatrix(d->GetNumberOfTuples(),d->GetNumberOfComponents(),mxREAL);
                double *datPtr = mxGetPr(dat);
                double *vtkDat = (double*)d->GetVoidPointer(0);
                
                for(int j=0;j<d->GetNumberOfComponents();j++)
                {
                    for(int k=0;k<d->GetNumberOfTuples();k++)
                    {
                        datPtr[k+j*d->GetNumberOfTuples()] = vtkDat[j+k*d->GetNumberOfComponents()];
                    }
                }
                
                mxSetFieldByNumber(CD,0,id,dat);
               // mxSetFieldByNumber(PD,0,id,mxCreateDoubleScalar(d->GetDataType()));
            }
        }
    }
    mxSetFieldByNumber(data,0,3,CD);
    
    
    // Now check for field data
    mxArray *FD = mxCreateStructMatrix(1,1,0,0);
    if (dataSet->GetFieldData())
    {
        for (int i = 0; i < dataSet->GetFieldData()->GetNumberOfArrays(); i++)
        {
            mxAddField(FD,dataSet->GetFieldData()->GetArray(i)->GetName());
        }
    }
    mxSetFieldByNumber(data,0,4,FD);
    
    dataSet->Delete();
    
    return data;
}

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.