2

This is my first time asking a question, so please be patient with me :)

I am writing a code in Python to reduce data that is stored as P(x,y,z), where P is a scalar intensity. I'd like to choose a path through space ([x1,y1,z1] -->[x2,y2,z2]) and a maximum radius of integration (rMax), and then obtain a list vector U(x,y,z), where U is the intensity integrated on the plane perpendicular to the path from 0< r < rMax. My problem is sorting the coordinates before handling the integration. I think the way to do the integration is to determine, for a set of points in a volume around the line, which point on the line is the closest.

This question is similar to Find nearest value in numpy array, but quite a bit more complex.

Let 

qpath = [[x1, y1, z1], [x2, y2, z2], ..., [xN,yN,zN]]

represent the path.

Let 

qvol = [[x'1, y'1, z'1, P(x'1, y'1, z'1)], 
    [x'1, y'1, z'2, P(x'1, y'1, z'2)], ..., 
    [x'1, y'2, z'1, P(x'1, y'2, z'1)], ..., 
    [x'N, y'N, z'N, P(x'N, y'N, z'N)]]

represent the coordinates and intensities of ALL points within the volume sufficiently near the path (closer than some value 'R').

Here is my attempt to proceed:

qdist = np.array([[[np.linalg.norm(i[0:3]-j), j[0], j[1], j[2], i[0], i[1], i[2]]
    for j in qpath] for i in qvol])

qplanes = np.array([j for j in i if j[0]==np.min(i[:,0]) 
    and j<R for i in qdist])

qdist is a rank 3 tensor, where the ith row on the kth page contains [d_ik, xi, yi, zi, x'k, y'k, z'k], where d_ik is the distance from the ith point on the line to the kth volume pixel.

qplanes should be a rank 2 tensor, where each row contains information from each element of qdist where the line is closest to the volume pixel.

Once I have qplanes, sorting and summing should be straight-forward.

Problems:

  1. Calculating qdist takes too long (about 10 minutes for a path between [1,1,1] --> [2,2,2]). I imagine there may be a more efficient way to approach this.
  2. The calculation of qplanes produces an error:

Invalid index to scalar variable

though I'm not sure why.

Oh great wise ones of the tubes, I beseech thee to please assist in this programming quandry. Many thanks to you if you venture on this quest. I am also appreciative of any feedback related to how I asked the question. Thank you so much.

Community
  • 1
  • 1
BohemianTapestry
  • 21
  • 2
  • Thinking about your index error: can you show an example of what `qdist` ends up looking like? – Engineero Aug 06 '15 at 22:27
  • 3
    Perhaps you should first rotate to a coordinate system `(x',y',z')` such that the line of interest is now the z axis, or some appropriate axis, for existing integration routines to be applied easily. – Paul Aug 06 '15 at 22:28
  • Is the scalar intensity `P` a function you can evaluate at an arbitrary point or do you have `P` as a three-dimensional array of values? – Till Hoffmann Aug 07 '15 at 00:04
  • Thanks for the questions. @Engineero: For a path from [1,1,1] to [2,2,2] and R=0.1, qdist = [[[0.1732, 1, 1, 1, 0.9, 0.9, 0.9], [0.2078, 1.02, 1.02, 1.02, 0.9, 0.9, 0.9], ..., ], [[0.1625, 1, 1, 1, 0.9, 0.9, 0.92], ..., ], ..., [[..., [0.1732, 1.98, 1.98, 1.98, 2.08, 2.08, 2.08]]]. – BohemianTapestry Aug 07 '15 at 13:30
  • @Paul: I've thought about rotating the data, but since the line may be arbitrary, the volume pixels would end up after rotation being at "non-nice" locations for integrating. But I will think about this more. – BohemianTapestry Aug 07 '15 at 13:31
  • @Till: P is a 3D array of values. P[0,0,0] is the intensity at the bottom-front-left corner, P[-1,-1,-1] is the intensity at the top-back-right corner. – BohemianTapestry Aug 07 '15 at 13:34
  • I think @Paul is correct that I should rotate the data. I am now working to 1.) rotate the data and 2) calculate a grid interpolation for integration. – BohemianTapestry Aug 10 '15 at 13:23

0 Answers0