将速度场可视化为定向箭头

Posted 2021-04-28

我有一个点云，对于每个点，我都有代表11个属性的11个数据数组，包括vx，vy，vz。我创建了一个新的数据数组，将vx，vy，vz存储为每个点的元组，并将其命名为v（即速度）。我想通过使用vtkGlyph3D将每个点的速度可视化为定向箭头，但它们都朝一个方向移动。我尝试了所有相关示例，但它们仍在同一方向上。我是vtk和C ++的新手。如果有人可以帮助我，我将不胜感激。

提前感谢！

这是我的代码：

//get vx, vy, vz
    vtkDataArray* vx0 = polyData->GetPointData()->GetArray("vx");
    vtkDataArray* vy0 = polyData->GetPointData()->GetArray("vy");
    vtkDataArray* vz0 = polyData->GetPointData()->GetArray("vz");

//define and initialize v and magnitude
    vtkDataArray* v = polyData->GetPointData()->GetArray("vz") ;
    vtkDataArray* magnitude = polyData->GetPointData()->GetArray("vz");
    v->SetName("v");

//insert v and magnitude to each point
    int numPoints = polyData->GetNumberOfPoints();  
    for (vtkIdType i = 0; i < numPoints; i++) { //define i; set loop condition
        double VX[3]; 
        double VY[3];
        double VZ[3];
        vx0->GetTuple(i, VX);
        vy0->GetTuple(i, VY);
        vz0->GetTuple(i, VZ);
        v->InsertTuple3(i, VX[0], VY[0], VZ[0]);

        double m = 0.0;
        m = sqrt(VX[0]*VX[0]+VY[0]*VY[0]+VZ[0]*VZ[0]);
        magnitude->InsertTuple1(i, m);
    }

    polyData->GetPointData()->AddArray(v);
    polyData->GetPointData()->SetVectors(v);
    polyData->GetPointData()->AddArray(magnitude);
    polyData->GetPointData()->SetScalars(magnitude);
    polyData->GetPointData()->SetActiveVectors("v");

//visualize
    vtkSmartPointer<vtkArrowSource> arrowSource = vtkSmartPointer<vtkArrowSource>::New();
    vtkSmartPointer<vtkGlyph3D> glyph3D = vtkSmartPointer<vtkGlyph3D>::New();
    glyph3D->SetSourceConnection(arrowSource->GetOutputPort());
    glyph3D->SetInputData(polyData);
    glyph3D->SetScaleFactor(0.01);
    glyph3D->OrientOn();
    glyph3D->SetVectorModeToUseVector();
    glyph3D->SetColorModeToColorByScalar();
    glyph3D->Update();


    vtkSmartPointer<vtkPolyDataMapper> sMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
    sMapper->SetInputConnection(glyph3D->GetOutputPort());
    sMapper->ScalarVisibilityOn();

    vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
    actor->SetMapper(sMapper);
    vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
    vtkSmartPointer<vtkRenderWindow> renderWindow = 
vtkSmartPointer<vtkRenderWindow>::New();
    renderWindow->AddRenderer(renderer);
    vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor = vtkSmartPointer<vtkRenderWindowInteractor>::New();
    renderWindowInteractor->SetRenderWindow(renderWindow);

    renderer->AddActor(actor);
    renderer->SetBackground(.3, .6, .3); // Background color green

    renderWindow->Render();
    renderWindowInteractor->Start(); 

答案

我认为您正在弄乱点数据数组：

1）在您的代码中，v和magnitude指向相同的对象：

vtkDataArray* v = polyData->GetPointData()->GetArray("vz") ;
vtkDataArray* magnitude = polyData->GetPointData()->GetArray("vz");

2）然后向该数组写入两次-首先是一个三元组，然后是一个元组：

v->InsertTuple3(i, VX[0], VY[0], VZ[0]);
magnitude->InsertTuple1(i, m);

相反，您应该使用两个不同的数据数组。这是一些python中的代码来说明用法。

import vtk
import numpy as np

# Create input.
poly, vel, mag = createData()

# Create point data.
velocity = vtk.vtkDoubleArray()
velocity.SetName("velocity")
velocity.SetNumberOfComponents(3)
velocity.SetNumberOfTuples(nPoints)
for i in range(nPoints):
    velocity.SetTuple(i, list(vel[i]))
magnitude = vtk.vtkDoubleArray()
# Similar as SetNumberOfValues(1) + SetNumberOfTuples(nPoints):
magnitude.SetNumberOfValues(nPoints)
magnitude.SetName("magnitude")
for i in range(nPoints):
    magnitude.SetValue(i, mag[i])

# Add to point data array.
poly.GetPointData().AddArray(velocity)
poly.GetPointData().AddArray(magnitude)
poly.GetPointData().SetActiveScalars("magnitude")
poly.GetPointData().SetActiveVectors("velocity")

# Create glyph.
glyph = vtk.vtkGlyph3D()
glyph.SetInputData(poly)
glyph.SetSourceConnection(arrow.GetOutputPort())
glyph.SetScaleFactor(0.1)
glyph.OrientOn()
glyph.SetVectorModeToUseVector()
glyph.SetColorModeToColorByScalar()
glyph.Update()

visualize()

为了完整起见，我还将提供其余代码。

def createData():
    source = vtk.vtkSphereSource()
    source.SetPhiResolution(30)
    source.SetThetaResolution(30)
    source.Update()
    poly = source.GetOutput()
    arrow = vtk.vtkArrowSource()
    arrow.Update()
    # Compute some velocity field. This data you already got.
    # I'm creating here some data myself.
    nPoints = poly.GetNumberOfPoints()
    points = [poly.GetPoint(i) for i in range(nPoints)]
    points = np.asarray(points)
    x,y,z = points.T
    vel = np.c_[x*np.cos(y*np.pi), y*np.sin(2*x*np.pi), np.cos(2*z*np.pi)]
    mag = np.linalg.norm(vel, axis=1)
    return poly, vel, mag

def visualize():
    # Visualization (my own toolbox)
    from vtkutils.item_renderer import ItemRenderer
    ren = ItemRenderer()
    ren.backgroundColor = "black"
    ren.addObject(glyph, showScalars=True, showArray="magnitude")
    ren.show()