Extracting normals from point data only

[fchirono]

Description

Hi, new PyVista user here. I have a mesh of points on the surface of an airplane engine blade, and I am trying to use PyVista to extract the outwards normal vectors at each point. I feel like this should be a fairly straightforward and well-posed problem, but it's proving surprisingly difficult to solve. I have tried looking at previous issues here and on StackOverflow, but couldn't find any questions similar to mine.

I can happily create a PolyData object with my points, as shown below with a toy dataset:

import numpy as np
import pyvista as pv

# %% *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
# build a generic airfoil using 'airfoils' package (available on PyPI)
# ---> https://airfoils.readthedocs.io/en/latest/index.html

# from airfoils import Airfoil
# foil = Airfoil.NACA4('4812')

# Nx = 102
# x_foil = np.linspace(0, 1, Nx//2)
# y_foil_up = foil.y_upper(x_foil)
# y_foil_lo = foil.y_lower(x_foil)

# # concatenate top and bottom
# x = np.concatenate((x_foil, x_foil[::-1]))
# y = np.concatenate((y_foil_up, y_foil_lo[::-1]))

# *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
# Alternatively, define data points directly:

x = np.array([ 0.  , 0.02, 0.04, 0.06, 0.08, 0.1 , 0.12, 0.14, 0.16, 0.18, 0.2 ,
               0.22, 0.24, 0.26, 0.28, 0.3 , 0.32, 0.34, 0.36, 0.38, 0.4 , 0.42,
               0.44, 0.46, 0.48, 0.5 , 0.52, 0.54, 0.56, 0.58, 0.6 , 0.62, 0.64,
               0.66, 0.68, 0.7 , 0.72, 0.74, 0.76, 0.78, 0.8 , 0.82, 0.84, 0.86,
               0.88, 0.9 , 0.92, 0.94, 0.96, 0.98, 1.  , 1.  , 0.98, 0.96, 0.94,
               0.92, 0.9 , 0.88, 0.86, 0.84, 0.82, 0.8 , 0.78, 0.76, 0.74, 0.72,
               0.7 , 0.68, 0.66, 0.64, 0.62, 0.6 , 0.58, 0.56, 0.54, 0.52, 0.5 ,
               0.48, 0.46, 0.44, 0.42, 0.4 , 0.38, 0.36, 0.34, 0.32, 0.3 , 0.28,
               0.26, 0.24, 0.22, 0.2 , 0.18, 0.16, 0.14, 0.12, 0.1 , 0.08, 0.06,
               0.04, 0.02, 0.  ])

y = np.array([  0.        ,  0.03805971,  0.05349321,  0.06597081,  0.07664938,
                0.08600053,  0.09427433,  0.10162394,  0.10815206,  0.11393193,
                0.11901819,  0.12345299,  0.1272697 ,  0.13049523,  0.13315164,
                0.13525721,  0.13682722,  0.13787454,  0.13841001,  0.13844282,
                0.13797963,  0.13710874,  0.13595412,  0.13452314,  0.13282156,
                0.13085448,  0.12862644,  0.12614146,  0.12340309,  0.12041439,
                0.11717802,  0.11369624,  0.10997091,  0.10600351,  0.10179516,
                0.09734664,  0.09265833,  0.0877303 ,  0.08256226,  0.07715356,
                0.07150319,  0.06560979,  0.05947162,  0.05308659,  0.04645218,
                0.03956552,  0.03242328,  0.02502174,  0.0173567 ,  0.00942353,
                0.00121706, -0.0012992 ,  0.00117944,  0.00354851,  0.00580537,
                0.00794761,  0.00997303,  0.01187967,  0.01366579,  0.01532991,
                0.01687078,  0.01828745,  0.01957923,  0.0207457 ,  0.02178678,
                0.02270265,  0.02349386,  0.02416127,  0.02470611,  0.02512998,
                0.02543486,  0.02562315,  0.02569767,  0.02566172,  0.02551906,
                0.02527399,  0.02493133,  0.02449652,  0.02397561,  0.02337537,
                0.0227033 ,  0.0219667 ,  0.02103652,  0.01982814,  0.0183556 ,
                0.01663352,  0.01467847,  0.01250933,  0.01014766,  0.00761818,
                0.00494959,  0.00217557, -0.00066369, -0.00351907, -0.00632877,
               -0.00901202, -0.01145764, -0.01350184, -0.01487876, -0.01508707,
               -0.01287165,  0.        ])

Nx = x.shape[0]

# *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-

# z-direction coordinates (span)
Nz = 51
z_max = 5
z = np.linspace(0, z_max, Nz)

# build 3D coordinate array
XYZ = np.zeros((Nx, Nz, 3))
XYZ[:, :, 0] = x[:, np.newaxis]
XYZ[:, :, 1] = y[:, np.newaxis]
XYZ[:, :, 2] = z[np.newaxis, :]

XYZ = XYZ.reshape((-1, 3))

# apply spanwise-varying twist - max 45deg at tip
theta_twist = -np.pi/4

XYZ[:, 0] = (np.cos(theta_twist*XYZ[:, 2]/z_max)*XYZ[:, 0] - np.sin(theta_twist*XYZ[:, 2]/z_max)*XYZ[:, 1])
XYZ[:, 1] = (np.sin(theta_twist*XYZ[:, 2]/z_max)*XYZ[:, 0] + np.cos(theta_twist*XYZ[:, 2]/z_max)*XYZ[:, 1])

# %% *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
# pass the numpy points to the PolyData constructor, plot points for visual
# inspection

airfoil = pv.PolyData(XYZ)
airfoil.plot()

I'm hoping to obtain the outwards-pointing normal vectors at each point, for both the top and the bottom of the airfoil. However, straight-up calling plot_normals on this airfoil object gives a KeyError: 'Normals', which is unfortunately not very descriptive of what is going wrong. I suppose it might be related to this object not having any lines or faces associated with it yet, but that's just a guess...

Here are some other things I have tried with this data:

1) delaunay_2d is very slow and/or gets stuck with my real data, and returns really weird results with this toy data. Here, it did a strange zig-zag pattern, connecting points at the top at one z-value to points at the bottom at the next z-value - hard to explain, but when you zoom in it should be fairly clear it's not what I want:

surface = airfoil.delaunay_2d()
surface.plot_normals(mag=0.1)

2) extract_surface outputs an identical PolyData obj to the input obj - i.e. it also has only points, and no lines or faces. I guess the input object doesn't have any surface to be extracted? In any case, it also leads to KeyError: 'Normals' when trying to call plot_normals on it:

surface = airfoil.extract_surface()
surface.plot_normals(mag=0.1)

3) reconstruct_surface unfortunately results in absolute gibberish for the surface :(

surface = airfoil.reconstruct_surface()
surface.plot()

4) delaunay_3d results in an UnstructuredGrid, which doesn't have the plot_normals attribute/method.

surface = airfoil.delaunay_3d(alpha=0.005, progress_bar=True)
surface.plot()

I'm getting frustrated because to me this sounds like it should be a "simple" problem, but it's proving surprisingly difficult. I'm currently trying to manually define the line connectivity list and see if that can "nudge" PyVista where to create the surfaces, but I'm not quite there yet - and I'm not sure if that's the right direction.

Does anyone has any suggestions on how to approach this problem, or can point me to something I might have missed? Thank you all so much!

[akaszynski]

Cool example! This is a perfect use case for using a pyvista.StructuredGrid

import numpy as np
import pyvista as pv

# %% *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
# build a generic airfoil using 'airfoils' package (available on PyPI)
# ---> https://airfoils.readthedocs.io/en/latest/index.html

# from airfoils import Airfoil
# foil = Airfoil.NACA4('4812')

# Nx = 102
# x_foil = np.linspace(0, 1, Nx//2)
# y_foil_up = foil.y_upper(x_foil)
# y_foil_lo = foil.y_lower(x_foil)

# # concatenate top and bottom
# x = np.concatenate((x_foil, x_foil[::-1]))
# y = np.concatenate((y_foil_up, y_foil_lo[::-1]))

# *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
# Alternatively, define data points directly:

x = np.array([ 0.  , 0.02, 0.04, 0.06, 0.08, 0.1 , 0.12, 0.14, 0.16, 0.18, 0.2 ,
               0.22, 0.24, 0.26, 0.28, 0.3 , 0.32, 0.34, 0.36, 0.38, 0.4 , 0.42,
               0.44, 0.46, 0.48, 0.5 , 0.52, 0.54, 0.56, 0.58, 0.6 , 0.62, 0.64,
               0.66, 0.68, 0.7 , 0.72, 0.74, 0.76, 0.78, 0.8 , 0.82, 0.84, 0.86,
               0.88, 0.9 , 0.92, 0.94, 0.96, 0.98, 1.  , 1.  , 0.98, 0.96, 0.94,
               0.92, 0.9 , 0.88, 0.86, 0.84, 0.82, 0.8 , 0.78, 0.76, 0.74, 0.72,
               0.7 , 0.68, 0.66, 0.64, 0.62, 0.6 , 0.58, 0.56, 0.54, 0.52, 0.5 ,
               0.48, 0.46, 0.44, 0.42, 0.4 , 0.38, 0.36, 0.34, 0.32, 0.3 , 0.28,
               0.26, 0.24, 0.22, 0.2 , 0.18, 0.16, 0.14, 0.12, 0.1 , 0.08, 0.06,
               0.04, 0.02, 0.  ])

y = np.array([  0.        ,  0.03805971,  0.05349321,  0.06597081,  0.07664938,
                0.08600053,  0.09427433,  0.10162394,  0.10815206,  0.11393193,
                0.11901819,  0.12345299,  0.1272697 ,  0.13049523,  0.13315164,
                0.13525721,  0.13682722,  0.13787454,  0.13841001,  0.13844282,
                0.13797963,  0.13710874,  0.13595412,  0.13452314,  0.13282156,
                0.13085448,  0.12862644,  0.12614146,  0.12340309,  0.12041439,
                0.11717802,  0.11369624,  0.10997091,  0.10600351,  0.10179516,
                0.09734664,  0.09265833,  0.0877303 ,  0.08256226,  0.07715356,
                0.07150319,  0.06560979,  0.05947162,  0.05308659,  0.04645218,
                0.03956552,  0.03242328,  0.02502174,  0.0173567 ,  0.00942353,
                0.00121706, -0.0012992 ,  0.00117944,  0.00354851,  0.00580537,
                0.00794761,  0.00997303,  0.01187967,  0.01366579,  0.01532991,
                0.01687078,  0.01828745,  0.01957923,  0.0207457 ,  0.02178678,
                0.02270265,  0.02349386,  0.02416127,  0.02470611,  0.02512998,
                0.02543486,  0.02562315,  0.02569767,  0.02566172,  0.02551906,
                0.02527399,  0.02493133,  0.02449652,  0.02397561,  0.02337537,
                0.0227033 ,  0.0219667 ,  0.02103652,  0.01982814,  0.0183556 ,
                0.01663352,  0.01467847,  0.01250933,  0.01014766,  0.00761818,
                0.00494959,  0.00217557, -0.00066369, -0.00351907, -0.00632877,
               -0.00901202, -0.01145764, -0.01350184, -0.01487876, -0.01508707,
               -0.01287165,  0.        ])

Nx = x.shape[0]

# z-direction coordinates (span)
Nz = 51
z_max = 5
z = np.linspace(0, z_max, Nz)

# build 3D coordinate array
XYZ = np.zeros((Nx, Nz, 3))
XYZ[:, :, 0] = x[:, np.newaxis]
XYZ[:, :, 1] = y[:, np.newaxis]
XYZ[:, :, 2] = z[np.newaxis, :]

# XYZ = XYZ.reshape((-1, 3))

# apply spanwise-varying twist - max 45deg at tip
theta_twist = -np.pi/4

XYZ[:, :, 0] = (np.cos(theta_twist*XYZ[:, :, 2]/z_max)*XYZ[:, :, 0]
                - np.sin(theta_twist*XYZ[:, :, 2]/z_max)*XYZ[:, :, 1])
XYZ[:, :, 1] = (np.sin(theta_twist*XYZ[:, :, 2]/z_max)*XYZ[:, :, 0]
             + np.cos(theta_twist*XYZ[:, :, 2]/z_max)*XYZ[:, :, 1])

# %% *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
# pass the numpy points to the PolyData constructor, plot points for visual
# inspection

airfoil = pv.StructuredGrid(XYZ[:, :, 0], XYZ[:, :, 1], XYZ[:, :, 2])
airfoil.plot(background='w', show_edges=True, pbr=True, metallic=1.0, roughness=0.4)

Plot using Physically Based Rendering

Mind if I add this as an example to our gallery?

[fchirono]

Ooh, that looks great! Thank you for your help - and no, I don't mind at all if you add this as an example :)

However, what I'm looking for are the point normals, and from the API reference it seems that the StructuredGrid class does not have any methods for retrieving the point normals - is that correct?

I'm wondering whether it's possible to convert this StructuredGrid object to PolyData in order to have access to the point_normals method. Could that be done?

[fchirono]

Ha, I think I found out how to get the vectors!

# extract PolyData representation of mesh, and use that to find normals
airfoil2 = airfoil.extract_surface()

# extract normal vectors as numpy array
normals = airfoil2.point_normals

# plot mesh with normal vectors and save a screenshot
plotter = pv.Plotter()
plotter.add_mesh(airfoil2) 
plotter.add_arrows(airfoil2.points, airfoil2.point_normals, mag=0.02)
plotter.show(auto_close=False)
plotter.show(screenshot='my_normal_vectors.png')

This is the data I was looking for, so my issue can now be considered formally closed! 😃 Thank you so much!

[One final related questions, if that's ok: is there a reason why the normals are pointing inwards? I see that I could use airfoil2.compute_normals with flip_normals=True or auto_orient_normals=True as arguments, but is there a general way to determine a priori whether the normals will be inward- or outward-facing?]

[akaszynski]

Ha, I think I found out how to get the vectors!

Nice work. Yes, you have to convert to a PolyData first and then work from there.

[One final related questions, if that's ok: is there a reason why the normals are pointing inwards? I see that I could use airfoil2.compute_normals with flip_normals=True or auto_orient_normals=True as arguments, but is there a general way to determine a priori whether the normals will be inward- or outward-facing?]

You might want to reverse the direction of the points when creating the structured grid. In regards to pointing inward and outward, there is some logic that allows you to perform automatic orientation of normals, but that requires a manifold surface and isn't always successful. I'd just reverse the order of points supplied to the StructuredGrid. That or use flip_normals.

[fchirono]

Thank you for the suggestions. I did succeed in getting outwards-pointing normals by calling either

airfoil2 = airfoil.extract_surface()
airfoil2 = airfoil2.compute_normals(flip_normals=True)

or

airfoil2 = airfoil.extract_surface()
airfoil2 = airfoil2.compute_normals(auto_orient_normals=True)

Both options work, so that's the happy end of my issue here :)

For completeness' sake, unfortunately PolyData.flip_normals() requires a triangular mesh, so it doesn't work in this case. I'm also not 100% sure what you meant by

You might want to reverse the direction of the points when creating the structured grid

I tried using pv.StructuredGrid(XYZ[::-1, ::-1, 0], XYZ[::-1, ::-1, 1], XYZ[::-1, ::-1, 2]) but that didn't work (normals are still pointing inwards). Did I misunderstand you?