1

I am unsure how to solve the following problem in Python:
I have a cylindrical surface which I want to mesh and then bin data to each block of the mesh. I have positional data (X,Y,Z) of rays that hit the cylindrical surface. After generating a mesh on the cylinder I have to count the number of rays (data points) in each block of my mesh.

The parameters of the cylinder are as follow:
radius = 0.1m, height = 0.15m, midpoint: [X=0,Y=0,Z=0]

To generate the mesh one could use the answer from How to generate regular points on cylindrical surface with the code copied below: `

import numpy as np
def make_cylinder(radius, length, nlength, alpha, nalpha, center, orientation):
"""
radius = radius of cylinder,
length = cylinder height, nlength = number of lengthwise divisions.
alpha = total degrees of cylinder, nalpha = number of circumferential divisions.
center = [X,Y,Z] coordinates of cylinder's midpoint.

"""
#Create the length array
I = np.linspace(0, length, nlength)
#Create alpha array avoid duplication of endpoints
#Conditional should be changed to meet your requirements
if int(alpha) == 360:
    A = np.linspace(0, alpha, num=nalpha, endpoint=False)/180*np.pi
else:
    A = np.linspace(0, alpha, num=nalpha)/180*np.pi

#Calculate X and Y
X = radius * np.cos(A)
Y = radius * np.sin(A)

#Tile/repeat indices so all unique pairs are present
pz = np.tile(I, nalpha)
px = np.repeat(X, nlength)
py = np.repeat(Y, nlength)

points = np.vstack(( pz, px, py )).T

#Shift to center
shift = np.array(center) - np.mean(points, axis=0)
points += shift

#Orient tube to new vector
#Grabbed from an old unutbu answer
def rotation_matrix(axis,theta):
    a = np.cos(theta/2)
    b,c,d = -axis*np.sin(theta/2)
    return np.array([[a*a+b*b-c*c-d*d, 2*(b*c-a*d), 2*(b*d+a*c)],
                     [2*(b*c+a*d), a*a+c*c-b*b-d*d, 2*(c*d-a*b)],
                     [2*(b*d-a*c), 2*(c*d+a*b), a*a+d*d-b*b-c*c]])

ovec = orientation / np.linalg.norm(orientation)
cylvec = np.array([1,0,0])

if np.allclose(cylvec, ovec):
    return points

#Get orthogonal axis and rotation
oaxis = np.cross(ovec, cylvec)
rot = np.arccos(np.dot(ovec, cylvec))

R = rotation_matrix(oaxis, rot)
return points.dot(R)

` What remains now is to loop through the "Ray data" to find the number of rays that hit each block of the mesh.

My initial thought process were as followed:

data = np.genfromtxt('ray_data.csv', delimiter=',') 
num_data_rows, num_data_cols = np.shape(data)

for i in range (num_data_rows): #Loop through the data

This is where I am stuck. As mentioned, the "Ray Data" is a CSV file that contains the positions (X,Y,Z) of each ray that hit the cylindrical surface. See the provided link: Ray data sample for cylindrical surface.

I just need to figure how to check where the rays fall within my mesh. The number of rays in each block will be multiplied by a constant (power per ray) to obtain the power (in Watts) in each block. This value is then divided by the area of the block to obtain the heat flux (W/m^2).

The final output I need is an array containing the centroid of each Block of the mesh with the corresponding heat flux value.

Any ideas how to solve this problem? I believe working with pandas is also an option.

Franco de Jager
  • 21
  • 7
  • 1
    You are only interested in the surface so I recommend to transform your data & cylinder surface into a 2 dimensional space – Milla Well Oct 05 '17 at 06:50
  • Transforming the data into a 2D space will still leave me with the problem of binning the rays into the grid and determining the centroid of each bin. – Franco de Jager Oct 05 '17 at 09:48
  • Your "rays" are only points. Am I correct that by "rays" you actually mean points of intersection of some ray (which I thought would be a vector) with the cylinder's surface? – Paul Brodersen Oct 05 '17 at 09:57
  • Also, if the diameter of the cylinder is much larger than your grid spacing in the major and minor axis, then I would 1) find the centroids of each bin (or simply re-interpret your bin corners as the centroids of a shifted grid) 2) find the points of intersection of your rays with the surface (which you might already have), 3) determine the corresponding bin by finding the nearest neighbor centroid for each point of intersection. The last step can be computed efficiently using KD trees (e.g. `scipy.spatial.cKDTree`). – Paul Brodersen Oct 05 '17 at 10:05

1 Answers1

0

As already proposed in the comment, I would recommend to transform your surface into a 2d space. That way, you can group your data easily.

import numpy as np
import pandas as pd

# generate some random rays (for which we will just assume they hit the surface)

rays = pd.DataFrame(
  np.random.uniform(-1,1,(8000,4)),
  columns=["x", "y", "z", "intensity"]
)

# transform x and y to polar coordinates while dropping the radius
# (again, just assuming they hit the surface)

rays["phi"] = rays.T.apply(lambda row: np.arctan2(row.y, row.x)).T

rays.head() now looks something like the following. z and phi is a representation of the rays in 2d. The intensity is what you called "power of ray".

x         y         z  intensity       phi
0 -0.237026 -0.634709 -0.889694   0.362156 -1.928199
1 -0.481137 -0.446912  0.687224   0.268080 -2.393056
2 -0.805538  0.068678  0.272009   0.990947  3.056541
3  0.549282 -0.330665  0.318683  -0.150776 -0.541886
4 -0.215676 -0.030922 -0.478929   0.408720 -2.999190

now, simply create bins and group the data. Finally sum over all intensities.

z_bins = np.arange(0, 1, .1)
phi_bins = np.arange(-np.pi, np.pi, np.pi/10)

result = rays.groupby([
  pd.cut(rays.phi, phi_bins), 
  pd.cut(rays.z, z_bins), 
]).intensity.sum()

The result.head() then looks as follows:

phi               z         
(-3.142, -2.827]  (0, 0.1]      0.719154
                  (0.1, 0.2]   -1.733479
                  (0.2, 0.3]    2.073013
                  (0.3, 0.4]    1.967453
                  (0.4, 0.5]    0.001312
Milla Well
  • 3,193
  • 3
  • 35
  • 50