pixelated artifacts in native vtk implementation compared to pyvista

[jonaspleyer]

Hey everyone,

I am having trouble, figuring out the reason for why I am seeing pixelated artifacts in one of my renders. I have a pipline which works well for my purpose and is implemented in pyvista. I wanted to recreate this pipeline in native vtk to obtain better performance. The performance is considerably better but the artifacts make the results unusable for me.

Comparison of pyvista (left) and vtk (right) results

vtk Implementation

#include <cmath>

#include <vtkActor.h>
#include <vtkCamera.h>
#include <vtkCompositeDataGeometryFilter.h>
#include <vtkCylinderSource.h>
#include <vtkExtractEdges.h>
#include <vtkGraphicsFactory.h>
#include <vtkMultiBlockDataSet.h>
#include <vtkNamedColors.h>
#include <vtkOpenGLRenderer.h>
#include <vtkNew.h>
#include <vtkPolyDataMapper.h>
#include <vtkPNGWriter.h>
#include <vtkProperty.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderer.h>
#include <vtkSphereSource.h>
#include <vtkTransform.h>
#include <vtkTransformPolyDataFilter.h>
#include <vtkWindowToImageFilter.h>
#include <vtkImageExport.h>

#include "render.h"

#define M_PI 3.14159265358979323846

namespace {
    vtkSmartPointer<vtkDataSet> CreateSphereDataSet(Vertex center, double radius)
    {
        double c[3] = {center.x, center.y, center.z};
        vtkNew<vtkSphereSource> sphere;
        sphere->SetCenter(c);
        sphere->SetRadius(radius);
        sphere->SetPhiResolution(30);
        sphere->SetThetaResolution(30);
        sphere->Update();
        return sphere->GetOutput();
    }

    vtkSmartPointer<vtkDataSet> CreateCylinderDataSet(
        double center[3],
        double direction[3],
        double radius,
        double height
    ) {
        double c[3] = {center[0], center[1], center[2]};
        vtkNew<vtkCylinderSource> cylinder;
        cylinder->SetHeight(height);
        cylinder->SetRadius(radius);
        cylinder->SetResolution(100);
        cylinder->Update();

        // Transform it
        vtkNew<vtkTransform> transform;
        transform->Identity();
        // Default orientation
        double orientation[3] = {0,1,0};
        double angle = std::acos(
            direction[0]*orientation[0] + direction[1]*orientation[1] + direction[2]*orientation[2]
        ) / 2 / M_PI * 360;
        double cross[3] = {
            orientation[1]*direction[2] - orientation[2]*direction[1],
            orientation[2]*direction[0] - orientation[0]*direction[2],
            orientation[0]*direction[1] - orientation[1]*direction[0]
        };
        transform->RotateWXYZ(angle, cross);
        transform->PostMultiply();
        transform->Translate(c);
        transform->Update();

        vtkNew<vtkTransformPolyDataFilter> transformPD;
        transformPD->SetTransform(transform);
        transformPD->SetInputData(cylinder->GetOutput());
        transformPD->SetInputConnection(cylinder->GetOutputPort());
        transformPD->Update();

        return transformPD->GetOutput();
    }
} // namespace

vtkNew<vtkMultiBlockDataSet> create_mesh(Agent agent) {
    // Create cell surfaces
    vtkNew<vtkMultiBlockDataSet> mesh;

    Vertex* vertices = agent.positions;
    double radius = agent.radius;

    mesh->SetBlock(0, CreateSphereDataSet(vertices[0], radius));
    for (int i=1; i<agent.n_vertices; i++) {
        double center[3] = {
            0.5 * (vertices[i].x + vertices[i-1].x),
            0.5 * (vertices[i].y + vertices[i-1].y),
            0.5 * (vertices[i].z + vertices[i-1].z)
        };
        double direction[3] = {
            vertices[i].x - vertices[i-1].x,
            vertices[i].y - vertices[i-1].y,
            vertices[i].z - vertices[i-1].z
        };
        double height = std::sqrt(
            (direction[0]*direction[0] + direction[1]*direction[1] + direction[2]*direction[2])
        );
        direction[0] /= height;
        direction[1] /= height;
        direction[2] /= height;
        mesh->SetBlock(2*i-1, CreateCylinderDataSet(center, direction, radius, height));
        mesh->SetBlock(2*i, CreateSphereDataSet(vertices[i], radius));
    }
    return mesh;
}

void render_img(Agent* agents, int n_agents, Camera camera, void* buffer)
{
    vtkNew<vtkRenderer> aren;
    aren->RemoveAllLights();
    aren->ResetCamera(0, camera.size_x, 0, camera.size_y, 0, camera.distance_z);
    aren->UseFXAAOff();
    aren->UseSSAOOff();
    aren->SetAmbient(0, 0, 0);
    aren->SetNearClippingPlaneTolerance(0.001);
    aren->SetAllocatedRenderTime(5000);
    aren->SetAutomaticLightCreation(false);
    aren->SetBackground(0, 0, 0);

    auto cam = aren->GetActiveCamera();
    cam->SetParallelProjection(true);
    cam->SetPosition(camera.size_x/2.0, camera.size_y/2.0, camera.distance_z);
    cam->SetFocalPoint(camera.size_x/2.0, camera.size_y/2.0, 0);
    cam->OrthogonalizeViewUp();

    cam->SetParallelScale(std::min(camera.size_x, camera.size_y) / 2);
    cam->SetClippingRange(0, 2 * camera.distance_z);
    aren->ResetCameraClippingRange(0, camera.size_x, 0, camera.size_y, 0, camera.distance_z);

    for (int i=0; i<n_agents; i++) {
        auto mesh = create_mesh(agents[i]);

        vtkNew<vtkExtractEdges> edges;
        edges->SetInputData(mesh);

        vtkNew<vtkCompositeDataGeometryFilter> polydata;
        polydata->SetInputConnection(edges->GetOutputPort());

        vtkNew<vtkPolyDataMapper> mapper;
        mapper->SetInputConnection(0, polydata->GetOutputPort(0));

        vtkNew<vtkActor> actor;
        actor->SetMapper(mapper);
        actor->GetProperty()->SetLineWidth(0);
        actor->GetProperty()->SetColor(agents[i].color);
        actor->GetProperty()->SetSpecular(0);
        actor->GetProperty()->SetSpecularPower(0);
        actor->GetProperty()->SetAmbient(1.0);
        actor->GetProperty()->SetLighting(false);
        actor->GetProperty()->SetInterpolationToFlat();
        actor->GetProperty()->SetPointSize(5);
        actor->GetProperty()->SetRepresentationToSurface();
        actor->GetProperty()->SetBaseIOR(0.0);
        actor->GetProperty()->SetMetallic(0.0);
        actor->GetProperty()->SetRoughness(0.0);
        actor->GetProperty()->SetAnisotropy(0.0);
        actor->GetProperty()->SetCoatIOR(0.0);
        actor->GetProperty()->SetCoatStrength(0.0);
        actor->GetProperty()->SetDiffuse(0.0);
        actor->GetProperty()->SetEdgeVisibility(false);
        actor->GetProperty()->SetVertexVisibility(false);
        actor->UseBoundsOff();
        aren->AddActor(actor);
    }

    int size_x = camera.size_x * camera.resolution;
    int size_y = camera.size_y * camera.resolution;

    vtkNew<vtkRenderWindow> win;
    win->AddRenderer(aren);
    win->SetOffScreenRendering(1);
    win->SetSize(size_x, size_y);
    win->SetMultiSamples(0);
    win->Render();

    vtkNew<vtkWindowToImageFilter> windowToImageFilter;
    windowToImageFilter->SetInputBufferTypeToRGB();
    windowToImageFilter->SetInput(win);
    windowToImageFilter->Update();

    vtkNew<vtkImageExport> exporter;
    exporter->SetInputConnection(windowToImageFilter->GetOutputPort());
    exporter->SetInputData(windowToImageFilter->GetOutput());
    exporter->Update();
    exporter->Export(buffer);

    vtkNew<vtkPNGWriter> writer;
    writer->SetFileName("screenshot.png");
    writer->SetInputConnection(windowToImageFilter->GetOutputPort());
    writer->Write();
}

pyvista Implementation

def __create_cell_surfaces(
    cells: dict[CellIdentifier, tuple[RodAgent, CellIdentifier | None]],
):
    for ident in cells.keys():
        meshes = []
        p = cells[ident][0].pos
        r = cells[ident][0].radius

        meshes.append(pv.Sphere(center=p[0], radius=r))
        for j in range(max(p.shape[0] - 1, 0)):
            # Add sphere at junction
            meshes.append(pv.Sphere(center=p[j + 1], radius=r))

            # Otherwise add cylinders
            pos1 = p[j]
            pos2 = p[j + 1]
            center = 0.5 * (pos1 + pos2)
            direction = pos2 - center
            height = float(np.linalg.norm(pos1 - pos2))
            cylinder = pv.Cylinder(
                center=center, direction=direction, radius=r, height=height
            )
            meshes.append(cylinder)
        # Combine all together
        mesh = pv.MultiBlock(meshes).extract_geometry()
        yield (ident, mesh)


def main_routine():
    plotter = pv.Plotter(
        off_screen=True,
        window_size=[
            int(np.round(domain_size[0] * ppm[0])),
            int(np.round(domain_size[1] * ppm[1])),
        ],
    )

    # Get camera info
    p1 = np.array([0, 0, 0])
    p2 = np.array([domain_size[0], domain_size[1], 0])
    p3 = np.array([0, domain_size[1], 0])
    rect = pv.Rectangle([p1, p2, p3])
    plotter.add_mesh(rect)
    bounds = (0.0, domain_size[0], 0.0, domain_size[1], 0.0, 0)
    pv.Plotter.view_xy(plotter, bounds=bounds)
    pv.Plotter.enable_parallel_projection(plotter)
    plotter.camera.tight(padding=0)
    plotter.camera.position = (*plotter.camera.position[:2], render_distance)
    camera = plotter.camera.copy()
    plotter.clear_actors()

    pv.Plotter.set_background(plotter, [render_settings.bg_brightness] * 3)

    for ident, cell in __create_cell_surfaces(cells):
        color = [render_settings.cell_brightness] * 3
        if colors is not None:
            color = colors[ident]
        actor = plotter.add_mesh(
            cell,
            show_edges=False,
            color=color,
            diffuse=render_settings.diffuse,
            ambient=render_settings.ambient,
            specular=render_settings.specular,
            specular_power=render_settings.specular_power,
            metallic=render_settings.metallic,
            pbr=render_settings.pbr,
            lighting=render_settings.lighting,
        )
        actor.UseBoundsOff()

    if not render_settings.render_mask:
        pv.Plotter.enable_ssao(plotter, radius=render_settings.ssao_radius)
        plotter.enable_anti_aliasing()
    else:
        plotter.disable_anti_aliasing()

    plotter.camera = camera
    img = np.array(plotter.screenshot())
    plotter.close()

[dgobbi (David Gobbi)]

The use of vtkExtractEdges in the C++ code is suspicious. Is it there intentionally? It will cause the mesh to be rendered as a wireframe, which doesn’t seem like what you want.

Regarding the mesh, you don’t need the spheres, because vtkCylinderSource has options to automatically cap the ends a cylinder with hemispheres:

        cylinder->CappingOn();
        cylinder->CapsuleCapOn();

[jonaspleyer]

Thanks that actually solved it. It also improved performance considerably.