I want to develop a shperical histogram for 3D object. Given a 3D object, I put it in a cylinder where its axis pass through the centroid of the object. The cylinder is sampled into a number N of control points. Than, for each a control point Pn, each point of the considered object is encoded in a spherical frame of reference centered in Pn with dimensions ρ
(from 0 to a suitable value), θ
(form 0 to 180 degree) and ϕ
(from 0 to 360 degree). Each polar coordinate is uniformly sampled into ten parts, obtaining a set of 1000 elements {(ρi, θj, ϕk) : 0 ≤ i , j , k ≤ 9}. This image illustrates the reference cylinder and some controls points.
data3D={{607.587, 401.119, 1140}, {607.587, 481.119, 1100}, {607.587,
561.119, 1020}, {647.587, 361.119, 1140}, {647.587, 361.119,
1300}, {647.587, 401.119, 1180}, {647.587, 401.119, 1340}, {647.587,
441.119, 1140}, {647.587, 481.119, 1020}, {647.587, 521.119,
980}, {647.587, 561.119, 980}, {647.587, 601.119, 980}, {647.587,
641.119, 980}, {687.587, 321.119, 1180}, {687.587, 361.119,
1060}, {687.587, 361.119, 1220}, {687.587, 361.119, 1380}, {687.587,
401.119, 1140}, {687.587, 401.119, 1300}, {687.587, 441.119,
1100}, {687.587, 441.119, 1260}, {687.587, 481.119, 1060}, {687.587,
481.119, 1220}, {687.587, 561.119, 1020}, {687.587, 601.119,
1020}, {687.587, 681.119, 940}, {727.587, 321.119, 1060}, {727.587,
321.119, 1220}, {727.587, 321.119, 1380}, {727.587, 361.119,
1100}, {727.587, 361.119, 1260}, {727.587, 361.119, 1420}, {727.587,
401.119, 1020}, {727.587, 401.119, 1180}, {727.587, 401.119,
1340}, {727.587, 441.119, 1060}, {727.587, 441.119, 1220}, {727.587,
441.119, 1380}, {727.587, 481.119, 1180}, {727.587, 481.119,
1340}, {727.587, 521.119, 1220}, {727.587, 721.119, 980}, {767.587,
281.119, 1220}, {767.587, 281.119, 1380}, {767.587, 321.119,
1180}, {767.587, 321.119, 1340}, {767.587, 361.119, 860}, {767.587,
361.119, 1020}, {767.587, 361.119, 1180}, {767.587, 361.119,
1340}, {767.587, 401.119, 740}, {767.587, 401.119, 900}, {767.587,
401.119, 1060}, {767.587, 401.119, 1220}, {767.587, 401.119,
1380}, {767.587, 441.119, 740}, {767.587, 441.119, 900}, {767.587,
441.119, 1060}, {767.587, 441.119, 1220}, {767.587, 441.119,
1380}, {767.587, 481.119, 820}, {767.587, 481.119, 980}, {767.587,
481.119, 1140}, {767.587, 481.119, 1300}, {767.587, 521.119,
1060}, {767.587, 521.119, 1220}, {767.587, 521.119, 1380}, {767.587,
561.119, 1220}, {807.587, 241.119, 1180}, {807.587, 241.119,
1340}, {807.587, 281.119, 1180}, {807.587, 281.119, 1340}, {807.587,
321.119, 980}, {807.587, 321.119, 1140}, {807.587, 321.119,
1300}, {807.587, 321.119, 1460}, {807.587, 361.119, 820}, {807.587,
361.119, 980}, {807.587, 361.119, 1140}, {807.587, 361.119,
1300}, {807.587, 361.119, 1460}, {807.587, 401.119, 620}, {807.587,
401.119, 780}, {807.587, 401.119, 940}, {807.587, 401.119,
1100}, {807.587, 401.119, 1260}, {807.587, 401.119, 1420}, {807.587,
441.119, 580}, {807.587, 441.119, 740}, {807.587, 441.119,
900}, {807.587, 441.119, 1060}, {807.587, 441.119, 1220}, {807.587,
441.119, 1380}, {807.587, 441.119, 1580}, {807.587, 481.119,
660}, {807.587, 481.119, 820}, {807.587, 481.119, 980}, {807.587,
481.119, 1140}, {807.587, 481.119, 1300}, {807.587, 481.119,
1460}, {807.587, 481.119, 1700}, {807.587, 521.119, 1020}, {807.587,
521.119, 1180}, {807.587, 521.119, 1340}, {807.587, 521.119,
1620}, {807.587, 561.119, 1060}, {807.587, 561.119, 1220}, {807.587,
561.119, 1380}, {807.587, 601.119, 1220}, {847.587, 201.119,
1260}, {847.587, 241.119, 1140}, {847.587, 241.119, 1300}, {847.587,
281.119, 1020}, {847.587, 281.119, 1180}, {847.587, 281.119,
1340}, {847.587, 321.119, 340}, {847.587, 321.119, 500}, {847.587,
321.119, 780}, {847.587, 321.119, 980}, {847.587, 321.119,
1140}, {847.587, 321.119, 1300}, {847.587, 321.119, 1460}, {847.587,
361.119, 380}, {847.587, 361.119, 540}, {847.587, 361.119,
700}, {847.587, 361.119, 860}, {847.587, 361.119, 1020}, {847.587,
361.119, 1180}, {847.587, 361.119, 1340}, {847.587, 361.119,
1500}, {847.587, 401.119, 420}, {847.587, 401.119, 580}, {847.587,
401.119, 740}, {847.587, 401.119, 900}, {847.587, 401.119,
1060}, {847.587, 401.119, 1220}, {847.587, 401.119, 1380}, {847.587,
401.119, 1540}, {847.587, 441.119, 140}, {847.587, 441.119,
300}, {847.587, 441.119, 460}, {847.587, 441.119, 620}, {847.587,
441.119, 780}, {847.587, 441.119, 940}, {847.587, 441.119,
1100}, {847.587, 441.119, 1260}, {847.587, 441.119, 1420}, {847.587,
441.119, 1580}, {847.587, 441.119, 1740}, {847.587, 481.119,
620}, {847.587, 481.119, 780}, {847.587, 481.119, 940}, {847.587,
481.119, 1100}, {847.587, 481.119, 1260}, {847.587, 481.119,
1420}, {847.587, 481.119, 1580}, {847.587, 481.119, 1740}, {847.587,
521.119, 820}, {847.587, 521.119, 980}, {847.587, 521.119,
1140}, {847.587, 521.119, 1300}, {847.587, 521.119, 1460}, {847.587,
521.119, 1660}, {847.587, 561.119, 1020}, {847.587, 561.119,
1180}, {847.587, 561.119, 1340}, {847.587, 561.119, 1620}, {847.587,
601.119, 1060}, {847.587, 601.119, 1220}, {847.587, 641.119,
1100}, {887.587, 201.119, 1100}, {887.587, 201.119, 1260}, {887.587,
201.119, 1420}, {887.587, 241.119, 1140}, {887.587, 241.119,
1300}, {887.587, 241.119, 1460}, {887.587, 281.119, 380}, {887.587,
281.119, 860}, {887.587, 281.119, 1020}, {887.587, 281.119,
1180}, {887.587, 281.119, 1340}, {887.587, 321.119, 100}, {887.587,
321.119, 260}, {887.587, 321.119, 420}, {887.587, 321.119,
580}, {887.587, 321.119, 740}, {887.587, 321.119, 980}, {887.587,
321.119, 1140}, {887.587, 321.119, 1300}, {887.587, 321.119,
1460}, {887.587, 361.119, 180}, {887.587, 361.119, 340}, {887.587,
361.119, 500}, {887.587, 361.119, 660}, {887.587, 361.119,
900}, {887.587, 361.119, 1060}, {887.587, 361.119, 1220}, {887.587,
361.119, 1380}, {887.587, 361.119, 1540}, {887.587, 401.119,
300}, {887.587, 401.119, 460}, {887.587, 401.119, 620}, {887.587,
401.119, 780}, {887.587, 401.119, 980}, {887.587, 401.119,
1140}, {887.587, 401.119, 1300}, {887.587, 401.119, 1460}, {887.587,
401.119, 1620}, {887.587, 441.119, 60}, {887.587, 441.119,
420}, {887.587, 441.119, 580}, {887.587, 441.119, 780}, {887.587,
441.119, 940}, {887.587, 441.119, 1100}, {887.587, 441.119,
1260}, {887.587, 441.119, 1420}, {887.587, 441.119, 1580}, {887.587,
441.119, 1740}, {887.587, 481.119, 140}, {887.587, 481.119,
700}, {887.587, 481.119, 860}, {887.587, 481.119, 1020}, {887.587,
481.119, 1180}, {887.587, 481.119, 1340}, {887.587, 481.119,
1500}, {887.587, 481.119, 1660}, {887.587, 521.119, 100}, {887.587,
521.119, 820}, {887.587, 521.119, 980}, {887.587, 521.119,
1140}, {887.587, 521.119, 1300}, {887.587, 521.119, 1460}, {887.587,
521.119, 1620}, {887.587, 561.119, 700}, {887.587, 561.119,
860}, {887.587, 561.119, 1020}, {887.587, 561.119, 1180}, {887.587,
561.119, 1340}, {887.587, 561.119, 1540}, {887.587, 561.119,
1700}, {887.587, 601.119, 1100}, {887.587, 601.119, 1260}, {887.587,
641.119, 1060}, {927.587, 161.119, 1100}, {927.587, 161.119,
1260}, {927.587, 201.119, 1140}, {927.587, 201.119, 1300}, {927.587,
241.119, 340}, {927.587, 241.119, 1060}, {927.587, 241.119,
1220}, {927.587, 241.119, 1380}, {927.587, 281.119, 140}, {927.587,
281.119, 300}, {927.587, 281.119, 900}, {927.587, 281.119,
1060}, {927.587, 281.119, 1220}, {927.587, 281.119, 1380}, {927.587,
321.119, 60}, {927.587, 321.119, 220}, {927.587, 321.119,
380}, {927.587, 321.119, 580}, {927.587, 321.119, 980}, {927.587,
321.119, 1140}, {927.587, 321.119, 1300}, {927.587, 321.119,
1460}, {927.587, 361.119, 100}, {927.587, 361.119, 260}, {927.587,
361.119, 420}, {927.587, 361.119, 580}, {927.587, 361.119,
780}, {927.587, 361.119, 980}, {927.587, 361.119, 1140}, {927.587,
361.119, 1300}, {927.587, 361.119, 1460}, {927.587, 401.119,
100}, {927.587, 401.119, 340}, {927.587, 401.119, 500}, {927.587,
401.119, 700}, {927.587, 401.119, 900}, {927.587, 401.119,
1060}, {927.587, 401.119, 1220}, {927.587, 401.119, 1380}, {927.587,
401.119, 1580}, {927.587, 441.119, 20}, {927.587, 441.119,
540}, {927.587, 441.119, 740}, {927.587, 441.119, 980}, {927.587,
441.119, 1140}, {927.587, 441.119, 1300}, {927.587, 441.119,
1460}, {927.587, 441.119, 1620}, {927.587, 481.119, 20}, {927.587,
481.119, 780}, {927.587, 481.119, 980}, {927.587, 481.119,
1140}, {927.587, 481.119, 1300}, {927.587, 481.119, 1460}, {927.587,
481.119, 1620}, {927.587, 521.119, 60}, {927.587, 521.119,
820}, {927.587, 521.119, 980}, {927.587, 521.119, 1140}, {927.587,
521.119, 1300}, {927.587, 521.119, 1460}, {927.587, 521.119,
1620}, {927.587, 561.119, 700}, {927.587, 561.119, 860}, {927.587,
561.119, 1020}, {927.587, 561.119, 1180}, {927.587, 561.119,
1340}, {927.587, 561.119, 1620}, {927.587, 601.119, 740}, {927.587,
601.119, 900}, {927.587, 601.119, 1140}, {927.587, 601.119,
1300}, {967.587, 121.119, 1180}, {967.587, 161.119, 1140}, {967.587,
201.119, 980}, {967.587, 201.119, 1140}, {967.587, 201.119,
1300}, {967.587, 241.119, 140}, {967.587, 241.119, 1180}, {967.587,
241.119, 1340}, {967.587, 281.119, 60}, {967.587, 281.119,
220}, {967.587, 281.119, 980}, {967.587, 281.119, 1140}, {967.587,
281.119, 1300}, {967.587, 321.119, 20}, {967.587, 321.119,
180}, {967.587, 321.119, 340}, {967.587, 321.119, 980}, {967.587,
321.119, 1140}, {967.587, 321.119, 1300}, {967.587, 321.119,
1460}, {967.587, 361.119, 140}, {967.587, 361.119, 300}, {967.587,
361.119, 460}, {967.587, 361.119, 620}, {967.587, 361.119,
1020}, {967.587, 361.119, 1180}, {967.587, 361.119, 1340}, {967.587,
361.119, 1500}, {967.587, 401.119, 140}, {967.587, 401.119,
700}, {967.587, 401.119, 900}, {967.587, 401.119, 1060}, {967.587,
401.119, 1220}, {967.587, 401.119, 1380}, {967.587, 441.119,
20}, {967.587, 441.119, 780}, {967.587, 441.119, 980}, {967.587,
441.119, 1140}, {967.587, 441.119, 1300}, {967.587, 441.119,
1460}, {967.587, 481.119, 60}, {967.587, 481.119, 860}, {967.587,
481.119, 1060}, {967.587, 481.119, 1220}, {967.587, 481.119,
1380}, {967.587, 481.119, 1620}, {967.587, 521.119, 740}, {967.587,
521.119, 940}, {967.587, 521.119, 1100}, {967.587, 521.119,
1260}, {967.587, 521.119, 1420}, {967.587, 561.119, 620}, {967.587,
561.119, 780}, {967.587, 561.119, 1020}, {967.587, 561.119,
1180}, {967.587, 561.119, 1340}, {967.587, 601.119, 540}, {967.587,
601.119, 700}, {967.587, 601.119, 860}, {967.587, 641.119,
420}, {967.587, 641.119, 580}, {967.587, 641.119, 780}, {1007.59,
161.119, 980}, {1007.59, 161.119, 1180}, {1007.59, 201.119,
980}, {1007.59, 201.119, 1140}, {1007.59, 241.119, 820}, {1007.59,
241.119, 1020}, {1007.59, 241.119, 1220}, {1007.59, 241.119,
1380}, {1007.59, 281.119, 1060}, {1007.59, 281.119, 1220}, {1007.59,
281.119, 1380}, {1007.59, 321.119, 100}, {1007.59, 321.119,
980}, {1007.59, 321.119, 1140}, {1007.59, 321.119, 1300}, {1007.59,
321.119, 1460}, {1007.59, 361.119, 140}, {1007.59, 361.119,
1020}, {1007.59, 361.119, 1180}, {1007.59, 361.119, 1340}, {1007.59,
401.119, 20}, {1007.59, 401.119, 580}, {1007.59, 401.119,
780}, {1007.59, 401.119, 1020}, {1007.59, 401.119, 1180}, {1007.59,
401.119, 1340}, {1007.59, 401.119, 1500}, {1007.59, 441.119,
820}, {1007.59, 441.119, 1100}, {1007.59, 441.119, 1260}, {1007.59,
441.119, 1420}, {1007.59, 481.119, 700}, {1007.59, 481.119,
940}, {1007.59, 481.119, 1100}, {1007.59, 481.119, 1260}, {1007.59,
481.119, 1420}, {1007.59, 521.119, 660}, {1007.59, 521.119,
820}, {1007.59, 521.119, 1020}, {1007.59, 521.119, 1180}, {1007.59,
521.119, 1340}, {1007.59, 561.119, 540}, {1007.59, 561.119,
700}, {1007.59, 561.119, 900}, {1007.59, 601.119, 380}, {1007.59,
601.119, 540}, {1007.59, 601.119, 700}, {1007.59, 641.119,
220}, {1007.59, 641.119, 380}, {1007.59, 641.119, 540}, {1007.59,
641.119, 700}, {1007.59, 641.119, 860}, {1007.59, 681.119,
220}, {1007.59, 681.119, 380}, {1007.59, 681.119, 540}, {1007.59,
681.119, 820}, {1007.59, 721.119, 180}, {1007.59, 721.119,
340}, {1007.59, 761.119, 60}, {1007.59, 841.119, 60}, {1047.59,
201.119, 1100}, {1047.59, 241.119, 1060}, {1047.59, 241.119,
1220}, {1047.59, 281.119, 1060}, {1047.59, 281.119, 1300}, {1047.59,
321.119, 1100}, {1047.59, 321.119, 1260}, {1047.59, 321.119,
1420}, {1047.59, 361.119, 1020}, {1047.59, 361.119, 1180}, {1047.59,
361.119, 1340}, {1047.59, 401.119, 820}, {1047.59, 401.119,
1060}, {1047.59, 401.119, 1220}, {1047.59, 401.119, 1380}, {1047.59,
441.119, 700}, {1047.59, 441.119, 860}, {1047.59, 441.119,
1100}, {1047.59, 441.119, 1260}, {1047.59, 441.119, 1420}, {1047.59,
481.119, 700}, {1047.59, 481.119, 860}, {1047.59, 481.119,
1100}, {1047.59, 481.119, 1260}, {1047.59, 481.119, 1420}, {1047.59,
521.119, 620}, {1047.59, 521.119, 780}, {1047.59, 521.119,
1020}, {1047.59, 521.119, 1180}, {1047.59, 561.119, 380}, {1047.59,
561.119, 540}, {1047.59, 561.119, 700}, {1047.59, 561.119,
860}, {1047.59, 601.119, 380}, {1047.59, 601.119, 540}, {1047.59,
601.119, 700}, {1047.59, 601.119, 860}, {1047.59, 641.119,
180}, {1047.59, 641.119, 340}, {1047.59, 641.119, 500}, {1047.59,
641.119, 660}, {1047.59, 681.119, 60}, {1047.59, 681.119,
220}, {1047.59, 681.119, 380}, {1047.59, 681.119, 540}, {1047.59,
721.119, 100}, {1047.59, 721.119, 260}, {1047.59, 721.119,
420}, {1047.59, 761.119, 140}, {1047.59, 801.119, 20}, {1047.59,
841.119, 20}, {1047.59, 881.119, 60}}
{cx, cy, cz} = Round[Mean[data3D]];
p1 = {cx, cy + 500, 880};
p2 = {cx, cy - 500, 880};
p3 = {cx + 500, cy, 880};
p4 = {cx - 500, cy, 880};
p5 = {cx + 354, cy + 354, 880};
p6 = {cx - 354, cy - 354, 880};
p7 = {cx + 354, cy - 354, 880};
p8 = {cx - 354, cy + 354, 880};
p9 = {cx + 468, cy + 177, 880};
p10 = {cx - 468, cy + 177, 880};
p11 = {cx + 468, cy - 177, 880};
p12 = {cx - 468, cy - 177, 880};
p13 = {cx + 177, cy + 468, 880};
p14 = {cx + 177, cy - 468, 880};
p15 = {cx - 177, cy + 468, 880};
p16 = {cx - 177, cy - 468, 880};
Show[Graphics3D[{PointSize[.015], Blue, Point[#] & /@ data3D}],
Graphics3D[{Opacity[.3],
Cylinder[{{cx, cy, 880}, {cx, cy, 881}}, 500]}],
Graphics3D[{PointSize[.02], Red, Point[p1]}],
Graphics3D[{PointSize[.02], Red, Point[p2]}],
Graphics3D[{PointSize[.02], Red, Point[p3]}],
Graphics3D[{PointSize[.02], Red, Point[p4]}],
Graphics3D[{PointSize[.02], Red, Point[p1]}],
Graphics3D[{PointSize[.02], Red, Point[p5]}],
Graphics3D[{PointSize[.02], Red, Point[p6]}],
Graphics3D[{PointSize[.02], Red, Point[p7]}],
Graphics3D[{PointSize[.02], Red, Point[p8]}],
Graphics3D[{PointSize[.02], Red, Point[p9]}],
Graphics3D[{PointSize[.02], Red, Point[p10]}],
Graphics3D[{PointSize[.02], Red, Point[p11]}],
Graphics3D[{PointSize[.02], Red, Point[p12]}],
Graphics3D[{PointSize[.02], Red, Point[p13]}],
Graphics3D[{PointSize[.02], Red, Point[p14]}],
Graphics3D[{PointSize[.02], Red, Point[p15]}],
Graphics3D[{PointSize[.02], Red, Point[p16]}]]
Now I need to get the spherical histogram corresponding to each control point Pn.
I start to code this using these 3 functions but I am not sure about the result because I get similar histograms for the different control points that having 0 values between 300 and 700.
1) convert data3D in spherical coordinate corresponding to each control points. 2) compute number of points that laying on each bins (ρi, θj, ϕk). (this function is to convert from cartesian to spherical coordinates)
cartesian2spherical[{x0_, y0_, z0_}, {x_, y_, z_}] :=
Module[{r, θ, ϕ}, (
(*0<=θ≤π,
0≤ϕ≤2π*)
r = N[Norm[{x, y, z} - {x0, y0, z0}]];
θ = Mod[N[ArcCos[(z - z0)/r]/Degree] + 180, 180];
If[(x - x0) == 0, ϕ =
Mod[N[ArcTan[(y - y0)/(0.00000000000000000000000001)]/Degree] +
360, 360], ϕ =
Mod[N[ArcTan[(y - y0)/(x - x0)]/Degree] + 360, 360]];
{r, θ, ϕ}
)]
newBinCounts[angles_, bins_] := Module[{hist, sectorIndex}, (
hist = BinCounts[angles, {bins}];
sectorIndex =
Table[Flatten[
Union[Position[angles, #] & /@
Select[angles, bins[[i]] <= # < bins[[i + 1]] &]]], {i, 1,
Length[bins] - 1}];
sectorIndex
)]
(histogram corresponding to each control point)
histogramPoint[p_, voxelset_] :=
Module[{coordinates, anglesϕ, anglesθ, raduis,
binsθ, binsϕ, binsr, sectorIndexϕ, sectorIndexr,
listofraduis, listofθ, hist, histogram, maxraduis}, (
(*Convertir les points de systeme cartisien vers le systeme \
spherique*)
coordinates = cartesian2spherical[p, #] & /@ voxelset;
anglesϕ = coordinates[[All, 3]];
anglesθ = coordinates[[All, 2]];
raduis = coordinates[[All, 1]];
maxraduis = Round[Max[raduis]];
(*pick bins for each parameters θ, ϕ,r*)
binsϕ = Range[0, 360, 360/10];
binsθ = Range[0, 180, 180/10];
binsr = Range[0, 1900, 1900/10];
(*Now for each ϕ for each θ for each r calculate the \
number of voxel in the correspondant volume*)
(*1. return the \
histogram according to ϕ and the index of each voxel laying in \
each sector ϕ*)
sectorIndexϕ = newBinCounts[anglesϕ, binsϕ];
(*recuperate raduis laying on each ϕ*)
listofraduis = raduis[[#]] & /@ sectorIndexϕ;
(*for each list of raduis laying on each ϕ,
we calculate index laying on each r*)
sectorIndexr = newBinCounts[#, binsr] & /@ listofraduis;
(*for each list of raduis on each ϕ,
we calculate the index laying on each θ*)
listofθ =
anglesθ[[#]] & /@ Flatten[sectorIndexr, 1];
hist = BinCounts[#, {binsθ}] & /@ listofθ;
histogram = Flatten[hist /. {} -> ConstantArray[0, 10]])]
listOfhistograms1 =
histogramPoint[#, data3D] & /@ {p1, p2, p3, p4, p5, p6, p7, p8, p9,
p10, p11, p12, p13, p14, p15, p16};
totalhist = Total[listOfhistograms1];
normalizedhistogram1 = N[totalhist/Max[totalhist]]
ListPlot[normalizedhistogram1, Joined -> True, DataRange -> All,
PlotRange -> All, AspectRatio -> 1/2, AxesOrigin -> 0,
Ticks -> {Range[0, 1000, 100], Automatic}]