pythonvtkparaview

Replace values in VTKCompositeDataArray


I use a programmable filter in Paraview to export a field named metrique for the purpose of remeshing. metrique is a VTKCompositeDataArray consisting in 9 values for each point of my computation results (3*3 matrix). I would like to replace values on the diagonal term of this VTKCompositeDataArray if passing under a threshold but I get the error "'VTKCompositeDataArray' object does not support item assignement".

Here is my python script (programmable filter in Paraview):

import numpy as np
import vtk
#Parameters
nb_points=395699
replacement_value=0.5
lower_threshold=0.5

#Reading quantities
input0=inputs[0]
Velocity=input0.PointData["Velocity"]

gradU=gradient(Velocity)
U=Velocity[:,0]

# Repeat the array along a new axis to create the desired shape
B = np.array([[[0.1,0,0],[0,0.1,0],[0,0,0.1]]])
B_tensorial = np.tile(B, (nb_points, 1, 1))

# Get indices of diagonal elements
diagonal_indices = np.diag_indices(B_tensorial.shape[1])

metrique=gradU*B_tensorial*U
# Apply the condition and replace values for diagonal elements
metrique[..., diagonal_indices[0], diagonal_indices[1]] = np.where(abs(metrique[..., diagonal_indices[0], diagonal_indices[1]]) > lower_threshold, replacement_value, metrique[..., diagonal_indices[0], diagonal_indices[1]])
output.PointData.append(abs(metrique), "metrique")

This code works well if I apply it to a Numpy array but not with the VTK array metrique. How can I replace the diagonal terms in the VTK array metrique following the wanted conditions ?

In response to one answer, I tried to create a function to convert the VTKCompositeDataArray into a Numpy array to work on it, but get the error File "", line 18, in vtk_composite_data_array_to_numpy AttributeError: 'VTKCompositeDataArray' object has no attribute 'GetNumberOfArrays'

def vtk_composite_data_array_to_numpy(composite_data_array): # Get the number of arrays in the composite data array num_arrays = composite_data_array.GetNumberOfArrays()

# Initialize an empty list to store the arrays
arrays = []

# Iterate over each array in the composite data array
for i in range(num_arrays):
    # Get the i-th array
    array = composite_data_array.GetArray(i)

    # Convert the vtkDataArray to a NumPy array
    num_values = array.GetNumberOfTuples() * array.GetNumberOfComponents()
    numpy_array = np.zeros(num_values, dtype=array.GetDataType())
    array.ExportToVoidPointer(numpy_array)

    # Reshape the 1D NumPy array into the appropriate shape
    numpy_array = numpy_array.reshape((array.GetNumberOfTuples(), array.GetNumberOfComponents()))

    # Append the NumPy array to the list
    arrays.append(numpy_array)

# Concatenate the arrays into a single NumPy array
numpy_array = np.concatenate(arrays, axis=0)

return numpy_array

Thanks a lot for your help!


Solution

  • VTKCompositeDataArray is a thin layer around a list of actual array, due to some memory management.

    You may iterate over each actual arrays through the metrique.Arrays list.

    See this doc for more info about this wrapping.

    edit: more details about Composite arrays

    CompositeArray is a collection of arrays by nature and because its related dataset is itself a collection of datasets under the hood. As each of the sub-dataset should have its own array defined in memory, you should not / cannot manipulate a single numpy array to cover all of them.

    The best way, if possible, is to use the algorithm.py module that wrap some of numpy. Otherwise, iterate over the arrays, do your numpy stuff, append numpy array in a list a fill back a composite object.

    import numpy
    
    output.ShallowCopy(inputs[0].VTKObject)
    
    composite_dataset = output
    composite_array = composite_dataset.PointData["RandomPointScalars"]
    
    # shortest way if available
    #res = add(composite_array, 2)
    
    # iterate and recreate composite array
    res_list = []
    for subarray in composite_array.Arrays:
      res_list.append(numpy.add(subarray, 2))
    res = dataset_adapter.VTKCompositeDataArray(res_list)
    
    # append to dataset attributes
    output.PointData.append(res, "res")