G[1]={{1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,-1}};
G[2]={{-(3/2),0,-1,(1/2)},{0,-1,0,0},{1,0,1,-1},{-(1/2),0,-1,-(1/2)}};
G[3]={{-(3/2),-((Sqrt[3])/2),(1/2),(1/2)},{((Sqrt[3])/2),(1/2),-((Sqrt[3])/2),-((Sqrt[3])/2)},{-(1/2),-((Sqrt[3])/2),-(1/2),(1/2)},{-(1/2),-((Sqrt[3])/2),(1/2),-(1/2)}};
G[4]={{-(3/2),((Sqrt[3])/2),(1/2),(1/2)},{-((Sqrt[3])/2),(1/2),((Sqrt[3])/2),((Sqrt[3])/2)},{-(1/2),((Sqrt[3])/2),-(1/2),(1/2)},{-(1/2),((Sqrt[3])/2),(1/2),-(1/2)}};

length=2; (*This parameter sets the number of reflections with the group elements*)
(matrices=
Nest[(
Union[#~Join~Flatten[Outer[Dot,{G[1],G[2],G[3],G[4]},#,1],1]])&,{IdentityMatrix[4]},length])//Length//Timing

Length[matrices]

ball0v[1]={1,Sqrt[(1-1/2)/2] Sin[θ]Cos[ϕ],Sqrt[(1-1/2)/2]Sin[θ]Sin[ϕ],Cos[θ] (1-1/2)/2+(1+1/2)/2};
ball0v[2]={1,Sqrt[(1-0.1428)/2] Sin[θ]Cos[ϕ],(1-0.1428)/2Sin[θ]Sin[ϕ]+(1+0.1428)/2,Sqrt[(1-0.1428)/2]Cos[θ]};
ball0v[3]={1,0.` -0.32733774606665816` Cos[ϕ] Sin[θ]+(1/2) (Sqrt[3]) (0.5714` +0.4286` Sin[θ] Sin[ϕ]),0.` -0.5669656074225313` Cos[ϕ] Sin[θ]+(1/2) (-0.5714`-0.4286` Sin[θ] Sin[ϕ]),0.` +0.6546754921333163` Cos[θ]};
ball0v[4]={1,0.` -0.32733774606665816` Cos[ϕ] Sin[θ]-(1/2) (Sqrt[3]) (0.5714` +0.4286` Sin[θ] Sin[ϕ]),0.` +0.5669656074225313` Cos[ϕ] Sin[θ]+(1/2) (-0.5714`-0.4286` Sin[θ] Sin[ϕ]),0.` +0.6546754921333163` Cos[θ]};

Balls=Drop[Flatten[Table[(matrices.G[j].ball0v[j])[[i]]/(matrices.G[j].ball0v[j])[[i,1]],{j,1,4},{i,1,Length[matrices]}],1],None,1];

Length[Balls]

graphics=ParametricPlot3D[Balls,{ϕ,  0, 2Pi},{θ,0, 2 Pi},PlotRange->1,Mesh->None,PlotStyle->Yellow,Axes->None,Boxed->False]

Export["BF2.pdf",graphics]