CPS_sampling.m 4.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990
  1. function [sph_pnt,energy] = CPS_sampling(varargin)
  2. p = inputParser;
  3. default.pnt_num = 240; % Number of particles; for this algorithm, more points the better
  4. default.thre = 'auto'; % Convergence threshold
  5. default.RepulF = 0.1; % Repulsive force factor, the higher the stronger force
  6. default.RstF = 1.5; % Resistance force, higher stronger
  7. default.MaxIter = 500; % Maximum number of iteration if not converged
  8. default.VisIt = false; % If make visulization
  9. addOptional(p,'pnt_num',default.pnt_num);
  10. addOptional(p,'thre',default.thre);
  11. addOptional(p,'RepulF',default.RepulF,@(n) n>0&&n<=1);
  12. addOptional(p,'RstF',default.RstF,@(n) n>0&&n<=5);
  13. addOptional(p,'MaxIter',default.MaxIter,@(n) n>0&&(floor(n)==n));
  14. addOptional(p,'VisIt',default.VisIt,@(n) islogical(n));
  15. parse(p,varargin{:});
  16. pnt_num = p.Results.pnt_num;
  17. thre = p.Results.thre;
  18. RepulF = p.Results.RepulF;
  19. RstF = p.Results.RstF;
  20. MaxIter = p.Results.MaxIter;
  21. VisIt = p.Results.VisIt;
  22. if strcmp(thre,"auto")
  23. thre = max(1/pnt_num,0.0065);
  24. end
  25. sph_pnt = rand(pnt_num,3)-.5; % Randomly assign particle position on 3D space
  26. sph_pnt = sph_pnt./sqrt(sum(sph_pnt.^2,2)); % Projected particles onto sphere
  27. sph_pntQN = quaternion([zeros(pnt_num,1),sph_pnt]); % Convert 3-element vector to quaternion
  28. currentSpeed = 0; % Initial moving speed = 0
  29. energy = 1; % The energy will be lower when the particles are more evenly distributed
  30. iterN = 0; % Iteration number
  31. while energy > thre && iterN <= MaxIter
  32. % Core for charged particle system
  33. diff_vec = sph_pntQN - (sph_pntQN.conj)'; % Calculate the differnece vector between each pair of particle
  34. vec_dist = asin(diff_vec.norm./2).*2; % Calculating the angle between each pair of particles
  35. proj_diff_vec = (diff_vec+sph_pntQN.*diff_vec.*sph_pntQN)./2; % Project the difference vectors onto the tangent plane for each particle
  36. vecForce = 1./vec_dist.^2.*proj_diff_vec./proj_diff_vec.norm; % Calculating the repulsive force vector from other particles weighted by 1/(the distance)^2
  37. vecForce = quaternion(squeeze(nansum(reshape(vecForce.compact,[size(proj_diff_vec),4]),2))); % Calculate the sum of all force vectors as the overall force vector
  38. vecForce = vecForce./pnt_num.*RepulF; % Adjust the force to an appriate level so the particles won't bounce around or move very slowly
  39. currentSpeed = currentSpeed+vecForce; % Updating the moving velocity for each particle
  40. rotAxis = sph_pntQN.*currentSpeed; % Compute the rotation axis for particle movement on sphere
  41. rotAxis = (rotAxis.compact.*[0 1 1 1]); % Extract only the cross product
  42. rotAxis = quaternion((rotAxis./sqrt(sum(rotAxis.^2,2))+[1 0 0 0])... % Compute the rotation quaternion
  43. .*[cos(vecForce.norm./2),sin(vecForce.norm./2).*[1,1,1]]); % vecForce.norm is the rotation angle for the movement
  44. sph_pntQN = rotAxis.*sph_pntQN.*rotAxis.conj; % Rotate the particles by the rotation quaternion
  45. currentSpeed = currentSpeed./RstF; % "Resistance force": prevent particles move too much
  46. energy_all = sort(vec_dist); % These two lines are for computing the energy after movement
  47. energy = std(energy_all(2,:));
  48. iterN = iterN+1;
  49. if rem(iterN,50) ==0
  50. energy
  51. end
  52. if iterN>3000&&energy>.01
  53. break
  54. end
  55. % Visualization part below
  56. if VisIt
  57. sph_pnt = sph_pntQN.compact;
  58. sph_pnt = sph_pnt(:,2:4);
  59. k = delaunayTriangulation(sph_pnt);
  60. [K,~] = convexHull(k);
  61. if iterN == 1
  62. clf
  63. sct = patch('vertices',sph_pnt,'faces',K,'Facecolor',[0,.25,.5],...
  64. 'FaceAlpha',0.7,'edgecolor',[0,.1,.3],'linewidth',2); % Visualize the solid formed with these particles
  65. axis equal vis3d off
  66. view(120,30)
  67. elseif rem(iterN,3) == 0
  68. sct.Faces = K;
  69. sct.Vertices = sph_pnt;
  70. end
  71. title(iterN)
  72. pause(0.001)
  73. end
  74. end
  75. if energy > thre
  76. warning("MaxIter reached; Not Converged")
  77. end
  78. sph_pnt = sph_pntQN.compact;
  79. sph_pnt = sph_pnt(:,2:4);
  80. end