02 — Scalar fields¶

Meshes carry data on them: a value per cell, or a value per point. This notebook covers attaching those fields, moving between the two, and getting them into NumPy.

You will learn:

• the difference between cell data and point data, and why it matters
• how to attach a NumPy array to a mesh as a named field
• how to convert cell -> point and point -> cell
• the one trap that catches everyone (these functions return a new mesh)

Self-contained: builds its own mesh, no data files.

Setup¶

In [ ]:

import numpy as np
import vtk

from imatools.core import mesh as cmesh

src = vtk.vtkSphereSource()
src.SetRadius(10.0)
src.SetThetaResolution(12)
src.SetPhiResolution(12)
src.Update()
sphere = src.GetOutput()

print(f"{sphere.GetNumberOfPoints()} points, {sphere.GetNumberOfCells()} cells")

import numpy as np import vtk from imatools.core import mesh as cmesh src = vtk.vtkSphereSource() src.SetRadius(10.0) src.SetThetaResolution(12) src.SetPhiResolution(12) src.Update() sphere = src.GetOutput() print(f"{sphere.GetNumberOfPoints()} points, {sphere.GetNumberOfCells()} cells")

Cell data vs point data¶

A triangulated sphere here has 240 cells but only 122 points. A field is defined on one or the other, so the array length tells you which:

Field	Length	Typical meaning
cell data	one per triangle	a measurement over a patch of surface
point data	one per vertex	a measurement at a location

imatools is cell-first: fibrosis_score defaults to type="cell".

Attaching a field¶

np_to_vtk_array wraps a NumPy array as a named VTK array. Adding it is not enough — you must also mark it active for functions that read "the scalars" to find it.

In [ ]:

rng = np.random.default_rng(42)
values = rng.uniform(0, 2, sphere.GetNumberOfCells()).astype(np.float32)

arr = cmesh.np_to_vtk_array(values, "scalars")
sphere.GetCellData().AddArray(arr)
sphere.GetCellData().SetActiveScalars("scalars")  # <- easy to forget

print("active cell scalars:", sphere.GetCellData().GetScalars().GetName())
print("tuples:", sphere.GetCellData().GetScalars().GetNumberOfTuples())

rng = np.random.default_rng(42) values = rng.uniform(0, 2, sphere.GetNumberOfCells()).astype(np.float32) arr = cmesh.np_to_vtk_array(values, "scalars") sphere.GetCellData().AddArray(arr) sphere.GetCellData().SetActiveScalars("scalars") # <- easy to forget print("active cell scalars:", sphere.GetCellData().GetScalars().GetName()) print("tuples:", sphere.GetCellData().GetScalars().GetNumberOfTuples())

Back to NumPy¶

In [ ]:

as_np = np.asarray(cmesh.convertCellDataToNpArray(sphere, "scalars"))

print("shape:", as_np.shape, "| mean:", round(float(as_np.mean()), 4))

as_np = np.asarray(cmesh.convertCellDataToNpArray(sphere, "scalars")) print("shape:", as_np.shape, "| mean:", round(float(as_np.mean()), 4))

Cell -> point¶

The trap. set_cell_to_point_data returns a new mesh. It does not modify the one you pass in. Ignore the return value and nothing appears to happen — the single most common mistake with this API.

In [ ]:

on_points = cmesh.set_cell_to_point_data(sphere, "scalars")  # capture the return!

ps = on_points.GetPointData().GetScalars()
print("new mesh point scalars:", ps.GetName(), "->", ps.GetNumberOfTuples(), "values")
print("original still has no point scalars:", sphere.GetPointData().GetScalars() is None)

on_points = cmesh.set_cell_to_point_data(sphere, "scalars") # capture the return! ps = on_points.GetPointData().GetScalars() print("new mesh point scalars:", ps.GetName(), "->", ps.GetNumberOfTuples(), "values") print("original still has no point scalars:", sphere.GetPointData().GetScalars() is None)

The 240 cell values are interpolated down to 122 point values — each vertex averages the cells touching it. This is lossy and not round-trippable.

In [ ]:

pn = np.asarray(cmesh.convertPointDataToNpArray(on_points, "scalars"))

print("cell values: ", as_np.shape, "mean", round(float(as_np.mean()), 4))
print("point values:", pn.shape, "mean", round(float(pn.mean()), 4))

pn = np.asarray(cmesh.convertPointDataToNpArray(on_points, "scalars")) print("cell values: ", as_np.shape, "mean", round(float(as_np.mean()), 4)) print("point values:", pn.shape, "mean", round(float(pn.mean()), 4))

Point -> cell¶

The inverse exists too, with the same return-a-new-mesh contract.

In [ ]:

back = cmesh.point_to_cell_data(on_points, "scalars")

cs = back.GetCellData().GetScalars()
print("round-tripped to cells:", cs.GetNumberOfTuples())
print("note: averaging twice does not recover the original values")

back = cmesh.point_to_cell_data(on_points, "scalars") cs = back.GetCellData().GetScalars() print("round-tripped to cells:", cs.GetNumberOfTuples()) print("note: averaging twice does not recover the original values")

Recap¶

Task	Call
NumPy -> VTK array	cmesh.np_to_vtk_array(a, "scalars")
Mark active	mesh.GetCellData().SetActiveScalars("scalars")
Cell field -> NumPy	cmesh.convertCellDataToNpArray(mesh, "scalars")
Point field -> NumPy	cmesh.convertPointDataToNpArray(mesh, "scalars")
Cell -> point	new = cmesh.set_cell_to_point_data(mesh, "scalars")
Point -> cell	new = cmesh.point_to_cell_data(mesh, "scalars")

Next: 03_scar_and_fibrosis — turning a scalar field into a scar score.