How to slice a sphere using spherical coordinates

I have scored the dose in spherical sectors created by dividing evenly the solid sphere along the radius, the azimuthal and polar angle. I need to slice the sphere selecting the sub-radius, sub-azimuthal angle, and sub-polar angle. The slice should not necessarily be parallel to the cartesian axes. I know command slice does that in cartesian coordinates. I would like to work in spherical coordinates to keep consistent with the equipment I am simulating. Is there a command or application to do that ?

Please answer to [e-mail]

Thank you so much maura

Answer 1

[John BG]

have you noticed that the points belonging to a sphere of radius 1 (surface points) can be defined with

    [x,y,z]=sphere

Then to gradually 'peal' the slices, shrink the surface with:

    h1=surf(.9*x,.9*y,.9*z)

the points of the spheric slice are contained in

[h1.XData h1.YData h1.ZData]

So, for your spherical processing from R1 inwards R2<R1 you could do something like:

R1=10
[x0 y0 z0]=sphere
[xR,yR,zR]=surf(R1*x0,R1*y0,R1*z0)  

dR =-.1
R=R1

while R>R2

       % do whatever process using xR yR zR
       % if here you need spheric coordinates 
       % use [azimuth,elevation,r] = cart2sph(xR,yR,zR)

       R=R-dR
       [xR,yR,zR]=surf(R*x0,R*y0,R*z0)
  end

Change R1 R2 and dR to your values of interest.

If you find this answer of any help solving your question, please click on the thumbs-up vote link,

thanks in advance

John

[e-mail]

Answer 2

[John BG]

1.- you want

R: radius [0 .8] cm R_step=.2

phi: azimuth (horizontal), range [0 2*pi] phi_step=36*pi/180

theta: polar (vertical), range [0 pi] theta_step=36*pi/180

2.- just in case, from

http://uk.mathworks.com/help/phased/ug/spherical-coordinates.html

the pair azimuth elevation and pair theta phi are commonly used in slightly different ways.

Also, there is the polar coordinates system, sometimes referred as cylindrical coordinates

the theta in polar is the phi in azimuth-elevation and theta-phi.

3.- slice is not constrained to slicing planes, you can also use spheres or the approximation of a sphere with the command sphere and surf

R1=2
R2=.2
dR=.2

[x0,y0,z0] = sphere(10);
[x,y,z] = meshgrid(-2:.2:2,-2:.2:2,-2:.2:2)
v = x.*exp(-x.^2-y.^2-z.^2);                         % test field, change to .csv inputs

slice(x,y,z,v,[-2,2],2,-2)

% colormap hot % default colormap is parula
pause(1)

for k= R1:-dR:R2
   hsp = surface(k*x0,k*y0,k*z0);
   xd = hsp.XData;   yd = hsp.YData;   zd = hsp.ZData;
   delete(hsp)

     hold on
     hslicer = slice(x,y,z,v,xd,yd,zd);
      axis tight;xlim([-3,3]);view(-10,35)
     drawnow

     pause(1)
     delete(hslicer)

     hold off
  end

modified example from <http://uk.mathworks.com/help/matlab/ref/slice.html?searchHighlight=slice>

as the sphere shrinks,

left and right values of the sphere surface show color according to 3D values of test field v.

Now that you see that slice can be used with other than flat surfaces let's swap the test values v with the readings in your .csv file

format long
data_in(:,5)=[]   % you said so
r=data_in(:,1)

4.- r=0 is only one point, yet your input file gives many readings for such point.

If you could get just one more reading for r=0 and for any r>0 then is when you use the angles currently assigned r=0.

That would either be an additional measurement no supplied in your .csv or one of the r=0 remains origin measurement, and the rest are really measurements already away from center.

A possible quick fix (assuming no r=0 done)

r=r+1

5.- From the azimuth-elevation theta-phi and polar screen shots above attached, it seems as the elevation (or phi) is positive above XY plane and negative below.

I don't really know the test instrument you have used to get the .csv supplied measurement, but just in case let me correct the angle range with

th=th-2

6.- the variable ranges are:

dr=.2;
r_range=dr*range(r)            % [min(r_range):dr:max(range)]=[0:dr:dr*4]

dph=36*pi/180;
ph_range=dph*range(ph)   % [min(ph_range),dph:max(ph_range)]=   [0:dph:2*pi]

dth=36*pi/180;
th_range= dth*range(th)     % [max(th_range),dth:min(th_range)]= [pi/2:dth:-pi/2]

L1=length(data_in)

7.- from spherical coordinates to cartesian. including both az-el and phi-theta just in case you find out that your readings are not reall az-el but phi-theta as defined above.

    using_phitheta=0

    switch using_phitheta

        case 1   % just in case when you said theta phi you referred to azimuth and elevation

        [az el]=phitheta2azel([ph*dph th*dth])
        [x1 y1 z1]=sph2cart(az,el,r*dr)
        data_in2=[x1 y1 z1 v_csv]

        case 0
            az=th;el=ph
            [x1 y1 z1]=sph2cart(az,el,r*dr)
            data_in2=[x1 y1 z1 v_csv]

    end

8.- You have supplied 250 measurements, and claim spherical sectors of solid 
angles 36º*36º Because you have divided both angle ranges by 10. 

This means your sectors look like this

   [x0,y0,z0] = sphere(10);
    surf(x0,y0,z0);

insert 30 36x36 100 solid angles

If all your r=0 measurements are same point, you have 10*3+1 sectors ~ 1000 the amount of expected measurements should be 1000.

If there is no r=0 measurement and all supplied measurements belong to a sector then you have 10*4 sectors = 1e4 expected amount of measurements.

Yet you have supplied 250 measurements only.

Replacing the 3D points of the test field v with .csv measurements is not difficult yet because there are less measurements than sectors, unless you explain how to interpolate, copy, or any other procedure to fill the gaps, any further step requires you supplying details.

I am copying this answer in the forum because other people may be interested to know how to use slice with concentric spheres rather than planes.

If you supply more details i can replace v with measurements and guessed values, but there has to be a clear procedure to do so.

Since your question is about the slice command and how to work with spherical coordinates, and considering there's further details needed for an interpolation policy, i respectfully consider that my explanation deserves, if you don't mind, accepting my answer by clicking on the thumbs-up vote in the MATLAB forum where you originally posted the question.

Thanks in advance,

let me know any further detail to refine question and answer.

John

[e-mail]

[Maura Monville]

The measurements are dose data obtained by a Monte Carlo code. I chose to score the dose in spherical sectors rather than in voxels of a rectilinear grid to get easily dose profiles through the central axis generated by slicing the sphefical grid.

In the attached output file the scoring spherical grid consists of 5 radius segments, 10 azimuthal segments, and 5 polar segments. You only showed how to slice the sphere along the radius therefore generating shrinking concentrical spheres.

I would like to slice the scoring grid with different planes as explained in the following. Assuming the spherical grid is centered at cartesian coordinates (X=0,Y=0,Z=0) and its diameter is 2*R then I would like to slice the spherical grid

along some transaxial planes: Z = 0, Z = +-0.5*R, Z = +-R, Z = +- 0.75*R

along some coronal planes: X = 0, X = +-0.5*R, X = +-R, X = +- 0.75*R

along some sagittal planes: Y = 0, Y = +-0.5*R, Y = +-R, Y = +- 0.75*R

along planes passing through (X=0, Y=0) incrementally rotated by Delta_Phi along the azimuthal angle like slicing an apple.

Furthermore, I do NOT know how to transform my data (dose data in spherical sectors like the one contained in the attached file) into a volume (v) as requested by the MatLab functions you called. Can you please explain ?

Please, feel free to create a data volume in spherical sectors (as many as you want) provided the the data structure is the same as the one in the attached file (Please, ignore the 5th column. The dose data is in the 4th column).

Thank you so much.

Maura

email: [e-mail] * * * bold***

Answer 3

[John BG]

Ok, it's a Montecarlo simulation: 5 elevation (polar) sectors, 10 azimuth sectors, 5 radius steps

How to translate the your .csv spherical coordinates:

1.- acquire data

format shortE
load data_1.mat
data_in(:,5)=[]           % no variance
data_in(2:5,:)=[]        % all initial 5 points are same measurement point, sphere centre

the relevant .csv file contents

data_in
 %    (  n_R      n_phi      n_theta    v   ) 
            0          0            0            0
            0   1.0000e+00            0   0
            0   1.0000e+00   1.0000e+00   2.6200e-14
            0   1.0000e+00   2.0000e+00            0
            0   1.0000e+00   3.0000e+00            0
            0   1.0000e+00   4.0000e+00            0
            0   2.0000e+00            0   3.6700e-14
            0   2.0000e+00   1.0000e+00            0
            ..
            0   9.0000e+00   3.0000e+00            0
            0   9.0000e+00   4.0000e+00   1.7800e-15
   1.0000e+00            0            0   4.7000e-15
   1.0000e+00            0   1.0000e+00   4.4500e-15
   ..
   1.0000e+00   9.0000e+00   3.0000e+00   1.7500e-14
   1.0000e+00   9.0000e+00   4.0000e+00            0
   2.0000e+00            0            0            0
   2.0000e+00            0   1.0000e+00   6.5800e-15
   2.0000e+00            0   2.0000e+00   4.0600e-15
   ..
   4.0000e+00   9.0000e+00   3.0000e+00   5.7900e-15
   4.0000e+00   9.0000e+00   4.0000e+00   2.1900e-15

get the different coordinates into variables

    r=data_in(:,1)
    ph=data_in(:,2)
    th=data_in(:,3)
    v_csv=data_in(:,4)

2.- scaling and corrections

v_csv=v_csv*1e15
    r=r+1
    th=th-2
    th=-th
    th(1)=0
    r(1)=0

3.- from numeral steps to spherical radius and angles

dr=.2     % unit: cm
dph=36  % unit: degree. 36*pi/180 rad
dth=36  % unit: degree. 36*pi/180 rad

r=r*dr
ph=ph*dph
th=th*dth

your data in spherical coordinates is this:

[r ph th v_csv]

4.- to avoid confusion, renaming ph (your phi) to el (elevation) and th (your theta) to az (azimuth)

el=th
az=ph
clear th ph data_sph

5.- your start data in elevation azimuth spherical coordinates has units:

% {cm degree degree 1e-15*[concentration of whatever malignant simulating]}

% start data in Cartesian coordinates

% R=(x^2+y^2+z^2)^.5
% az=atan(y/x)
% el=atan(z/(x^2+y^2)^.5)
%
% x=R*cos(el)*cos(az)
% y=R*cos(el)*sin(az)
% z=R*sin(el)
format bank
data_spheric_elaz=[r el az v_csv]   % DATA IN SPHERICAL COORDINATES HERE

6.- now your data looks like this (note that elevation (your theta) and azimuth (your phi) columns are swapped

     0             0             0             0
          0.20         72.00         36.00             0
          0.20         36.00         36.00         26.20
          0.20             0         36.00             0
          0.20        -36.00         36.00             0
          0.20        -72.00         36.00             0
          0.20         72.00         72.00         36.70
          0.20         36.00         72.00             0
        ..
          0.20        -36.00        324.00             0
          0.20        -72.00        324.00          1.78
          0.40         72.00             0          4.70
          0.40         36.00             0          4.45
         ..
          0.40        -36.00        324.00         17.50
          0.40        -72.00        324.00             0
          0.60         72.00             0             0
          0.60         36.00             0          6.58
          0.60             0             0          4.06
         ..
          1.00        -36.00        324.00          5.79
          1.00        -72.00        324.00          2.19

probably, where you now read r=1 you would like to read r=.8, where r=.8 you would like it to be r=.6 and so on but you gave r=0 to both origin, and to group of different measurements that could not possibly be same point, so the 'crust' around [0 0 0], the smallest radius different than 1 with different measurement values should have r>0

7.- and the [x y z] you need are calculated here:

format long
[x y z]=sph2cart(az,el,r)
data_cart=[x y z v_csv]   % START PLANE SLICING HERE

8.- in your second clarification you define the planes you are after, all you have to do is to select the data belonging to the planes you want.

right now you don't need anything else to proceed

% try the following, if you get stuck i will help you, but try first:

% along some transaxial planes: Z = 0, Z = +-0.5*R, Z = +-R, Z = +- 0.75*R

% along some coronal planes: X = 0, X = +-0.5*R, X = +-R, X = +- 0.75*R

% along some sagittal planes: Y = 0, Y = +-0.5*R, Y = +-R, Y = +- 0.75*R

%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following, you don't really need it right now yet you may later on:

polar coordinates

% data in polar coordinates [pol_rho pol_theta z], 
% please note this 'polar' is different than your 'polar'
% at least from the way you initially use this term to refer to elevation.

[pol_theta pol_rho]=cart2pol(x,y)
data_polar=[pol_rho pol_theta z v_csv]

phi theta spherical coordinates, not the phi theta initially used

    % data in phi theta spherical coordinates
    % sin(el)=sin(ϕ)*sin(θ)
    % tan(az)=cos(ϕ)*tan(θ)
    % cosθ=cos(el)*cos(az)
    % tanϕ=tan(el)/sin(az)

      % probably not used that often so the input format slightly different
    % draft az2=az-180;azel=[az2'; el']
    % draft phitheta1=azel2phitheta(azel)
    % draft th=phitheta1(1,:)  % THIS NEEDS CHECK, haven't so far
    % draft ph=phitheta1(2,:)  % THIS NEEDS CHECK, haven't so far

data_shperic_phitheta=[r ph th v_csv]

u v spherical coordinates

    % data in u v coordinates 
    % u=sin(θ)*cos(ϕ)
    % v=sin(θ)*sin(ϕ)
    % u=cos(el)*sin(az)
    % v=sin(el)
  azel2uv()

If you find this answer of any help solving your question, please click on the thumbs-up vote link,

thanks in advance

John

[e-mail]

[John BG]

Just to close my answer. either use thick planes for the slicing or interpolate

Since you just translated spherical points, but they are only aligned for the planes where Z=0, or Y=0 or X=0.

In the following figure your measurements / simulations are on the vertices. This would be for instance the view of any plane containing axis Z:

You have 2 choices

1.- define the planes, not as z= .. but as z>=z1 && z<=z2

or

2.- interpolate: There may be advanced literature / medical research that uses interpolation in similar case, if so it's your job to find the algorithm used to extend data points without adding noise, or too much of it.

Ask the ones involved in relevant research: Interpolating (aka guessing) may be a problem: To guess the dose on hundreds points of Polonium 210 intoxicated tissue, if you haven't really

Measured them, you might end up telling a patient in death row he/she is ok when they are nor, or telling patients that are not that bad that they don't have much left.

You really have to find out what medical research is on regarding whether interpolation is currently considered for the kind of material dosage and tissue you are simulating. Interpolating on certain tissue may be ok and accepted practice, but wrong on another type of tissue.

[Maura Monville]

Thank you. Some questions and explanations:

(1) Why is scaling and corrections necessary ? I understand dose values are too small to be visualized with colors. This is because the output file was generated by a MC simulation transporting only a few particles. Moreover, the MC code I used is written in C++ so indices start from 0 rather than from 1.

(2)Dose has units of Gy = [Joule/Kg] this is the energy deposited in the patient or phantom being irradiated.

(3)Sorry. I omitted an important detail about how the MC code outputs data in a CSV file. Consider the spherical indices (R,Phi,Theta) then Theta is the fastest coordinate, followed by Phi, followed by R. Likewise if the MC had scored the dose in the voxels of a rectilinear grid (X,Y,Z) then Z would be the fastest index, followed by Y, followed by Z. Like in a nested loop: for R=0:1:4 for Phi=0:36:260 for Theta=0:36:180 fprintf(out,'%f %f %f %f\n',R,Phi,Theta, Dose) end end end

(4) In order to slice the sphere through any plane having the data in data_cart I have (I think) to reshape the v_csv into a 3D matrix indexed through x,y,z because v_csv is just a 1D vector. Here I am stuck. I would possibly guess how to reshape v_csv if indexed by r,el,az but I am confused about how to reshape it after transforming from spherical to cartesian coordinates. Can you please help ?

(5) How can I slice the sphre through its meridians ? Which coordinates system should I use in this case ?

Thank you so much. Best regards maura

[e-mail]