title: Find maxima on a surf plot using Particle Swarm Optimization
body: <p>I'm new to Particle swarm optimization.</p>
  
  <p>I have a surf plot of a 2D image as shown below:
  <img src="https://i.stack.imgur.com/ZXykg.jpg" alt="enter image description here"></p>
  
  <p>And is it possible to fire up the particles randomly and make them find the maxima (global)?</p>
  
  <p>If yes, Can anyone please provide an algorithm or a sample code.</p>
  
  <p>Thank you.</p>
  
  <p>EDIT:
  Did some coding and got the answer, for a small number of iterations like 10-15 the particles swarm at the maxima but if the iterations are made more than 15 like 50 etc, the particles move away from the maxima.</p>
  
  <p>Can someone give a solution for it, I've included the matlab code and a sample image:</p>
  
  <pre><code>img = imread('gray_scale_img.jpg');
  %img is grayscale
  
  a=10; b=0.2; %a is convergence and b is convergence rate
  % range=[xmin xmax ymin ymax zmin zmax];
  range=[0 440 0 228 0 255]; %&lt;&lt;&lt;---size of the image
  
  Num_iterations = 15;%&lt;&lt;&lt;----problem exists if more iterations
  n = 30; % number of particles
  best=zeros(Num_iterations,3);
  
  % ----- Start Particle Swarm Optimization -----------
  % generating the initial locations of n particles
  xrange=range(2)-range(1);
  yrange=range(4)-range(3);
  zrange=range(6)-range(5);
  xn=(rand(1,n)*xrange+range(1));
  yn=(rand(1,n)*yrange+range(3));
  
  % Start iterations
  for k=1:Num_iterations,
  % Show the contour of the function
  surfc(double(img)); hold on;
  % Find the current best location (xo,yo)
  
  for d = 1:n
  xn = ceil(xn);
  yn = ceil(yn);
  zn(d)=img(yn(d),xn(d));
  end
  
  zn_min=max(max(zn));
  xo=max(max(xn(zn==zn_min)));
  yo=max(max(yn(zn==zn_min)));
  zo=max(max(zn(zn==zn_min)));
  
  % Trace the paths of all roaming particles
  % Display these roaming particles
  plot3(xn,yn,zn,'.',xo,yo,zo,'rp'); axis(range);
  
  % Move all the particles to new locations
  nn=size(yn,2);  %a=alpha, b=beta
  xn=xn.*(1-b)+xo.*b+a.*(rand(1,nn)-0.5);
  yn=yn.*(1-b)+yo.*b+a.*(rand(1,nn)-0.5);
  %zn=zn.*(1-b)+zo.*b+a.*(rand(1,nn)-0.5);
  
  %To make sure the particles are well within the range
  nn=length(yn);
  for l=1:nn,
    if xn(l)&lt;=range(1), xn(l)=range(1); end
    if xn(l)&gt;=range(2), xn(l)=range(2); end
    if yn(l)&lt;=range(3), yn(l)=range(3); end
    if yn(l)&gt;=range(4), yn(l)=range(4); end
    %if zn(l)&lt;=range(5), yn(l)=range(5); end
    %if zn(l)&gt;=range(6), yn(l)=range(6); end
  end
  drawnow;
  hold off;
  best(k,1)=xo; best(k,2)=yo; best(k,3)=zo;pause(0.2);
  end   %%%%% end of iterations
  
  x_best = best(k,1);
  y_best = best(k,2);
  z_best = best(k,3);
  global_maxima_coordinate =[x_best,y_best]
  
  figure(3);
  imshow(A,[]);
  hold on
  plot(x_best, y_best, 'bx') 
  % Draw circle of required radius around the global maxima
  theta = 0:pi/50:2*pi;
  radius_1 = 15;
  xunit_circle_1 = radius_1 * cos(theta) + x_best;
  yunit_circle_1 = radius_1 * sin(theta) + y_best;
  h_1 = plot(xunit_circle_1, yunit_circle_1,'r');
  hold off
  title('Global maxima');
  </code></pre>