finding multiple maxima in a Gaussian distribution

Hi

I have a query and will appreciate your help.

If we try the code below,

>> c1=1;
>> for p=-50:50
c2=1;
for q=-50:50
pp1(c1,c2)=exp(- ( (p-10).^2 + (q-9).^2) /2) + exp(- ( (p+20).^2 + (q-30).^2) /2);
c2=c2+1;
end
c1=c1+1;
end
>> mesh(pp1)

you can see that the figure(matrix Po) has two peaks(maxima), stating that there are two targets, one at (10,9) and another at (-20,30).

If there is only one target(peak), I used [a,b]=find(Po==max(max(Po))), to find the locations (i,j) at the peak.

But in the above code, we have two peaks. Can you help me to find the

1. No of peaks(estimate no of targets)

2. Find the location of the peaks(i,j)

Cheers

[Information merged from duplicate Question]

Hello

I have a simple query, which was unanswered a few days ago. I am re-posting because I need a hint urgently. Please do not direct me to my previous post. It is not quite complete.

You can copy this code and see the result.

s1 = [5 4 0.5  -0.64]';
s2 = [0 -3.9 -0.5  -0.64]';

Sa = [-40:0.2:40];
Sb = [-40:0.2:40];
L = length(Sa);
m = zeros(L,L);
Pr = 2*ones(L,L);
Pr = Pr / sum(sum(Pr));

      x1 = [s1(1);s1(2)] + sqrt(R)*randn(2,1);
      x2 = [s2(1);s2(2)] + sqrt(R)*randn(2,1);

      for i = 1: length(Pr)
          for j = 1: length(Pr)
              m(i,j) = (1/sqrt(2*pi)) * (exp(-((Sa(i)-x1(1)).^2 + (Sb(j)-x1(2)).^2)/2) + exp(-((Sa(i)-x2(1)).^2 + (Sb(j)-x2(2)).^2)/2));
              m(i,j) = m(i,j) * Pr(i,j);
          end
      end
      Po = m/sum(sum(m));
      [a,b] = find(Po==max(max(Po)));
      Sest = [Sa(a);Sb(b)];
  meash(Sa,Sb,Po)

The matrix(surfing of the matrix) 'Po', will result in two peaks at the (s1(1),s1(2)) and (s2(1),s2(2)), show that there are two targets. If I had one peak, a,b would collect it and Sest would collect the location on the grid.

But I am unable to figure out how to find two peaks instead. So the last two statements do not hold now.

I look forward for suggestions on

1. finding the number of peaks(2 in this case)

2. finding the location on grid(Sest in this case) (Sa,Sb)

Will appreciate your help!

Cheers

[Walter Roberson]

http://www.mathworks.com/matlabcentral/answers/29922-why-your-question-is-not-urgent-or-an-emergency

Answer 1

[Image Analyst]

See the code I added to yours:

c1=1; 
for p=-50:50 c2=1; 
    for q=-50:50 
        pp1(c1,c2)=exp(- ( (p-10).^2 + (q-9).^2) /2) + exp(- ( (p+20).^2 + (q-30).^2) /2); 
        c2=c2+1; 
    end
    c1=c1+1;
end
mesh(pp1)

Here is code I added:

[labeledImage numberOfPeaks] = bwlabel(pp1>0.01);
measurements = regionprops(labeledImage, 'Centroid');
message = sprintf('The number of peaks is %d\n', numberOfPeaks);
uiwait(msgbox(message));

for r = 1 : numberOfPeaks
    message = sprintf('The centroid of peak %d is at (%.2f, %.2f)\n', ...
        r, measurements(r).Centroid(1),  measurements(r).Centroid(2));
    uiwait(msgbox(message));
end

Yes I know the numbers aren't the same. I'm hoping you'll be able to figure out how to shift the coordinate system by subtracting 51 from x, and y.

[Image Analyst]

No I can't - I tried. You didn't define A, B, Q, and maybe some other variables. Please, make it easy for us. Make it so we CAN simply copy and paste.

[Image Analyst]

Here we fail again. Did you read my comment about not being able to run your code because you didn't supply A, B, and Q? Well you replied and still didn't supply those. Why not. Please don't post code again with unrelated errors. Are you asking about "Undefined function or variable 'A'."? No, so don't ask us to fix errors like that first just in order to get to what you really want us to fix. Like I said before "Please, make it easy for us. Make it so we CAN simply copy and paste."

[PChoppala]

So sorry, it was a mistake(copied it wrong), unintentionally done. Very sorry. I am not asking you to fix an error. I am giving you a clear view of my prob so you might guide in the right direction. Thanks for your support.

Here you go

s1 = [5 4 0.5 -0.64]';

s2 = [0 -3.9 -0.5 -0.64]';

Sa = [-40:0.2:40];

Sb = [-40:0.2:40];

L = length(Sa);

m = zeros(L,L);

Pr = 2*ones(L,L);

Pr = Pr / sum(sum(Pr));

R = [0.04 0; 0 0.04];

x1 = [s1(1);s1(2)] + sqrt(R)*randn(2,1);

x2 = [s2(1);s2(2)] + sqrt(R)*randn(2,1);

for i = 1: length(Pr)

for j = 1: length(Pr)

m(i,j) = (1/sqrt(2*pi)) * (exp(-((Sa(i)-x1(1)).^2 + (Sb(j)-x1(2)).^2)/2) + exp(-((Sa(i)-x2(1)).^2 + (Sb(j)-x2(2)).^2)/2));

m(i,j) = m(i,j) * Pr(i,j);

end

end

Po = m/sum(sum(m));

[a,b] = find(Po==max(max(Po)));

Sest = [Sa(a);Sb(b)];

figure(2);

mesh(-40:0.2:40, -40:0.2:40, Po);

Will appreciate guidance in the direction of finding [a,b] and Sest for the two targets.

I am looking for [a1,b1], [a2,b2] and Sest1 and Sest2.

Thanks a lot again, and sorry for the previous mistake .

[Image Analyst]

I already gave you the answer. Here it is again (add this code after yours):

[labeledImage numberOfPeaks] = bwlabel(Po > 0.001);

measurements = regionprops(labeledImage, 'Centroid');

message = sprintf('The number of peaks is %d\n', numberOfPeaks);

uiwait(msgbox(message));

for r = 1 : numberOfPeaks

% Name them a and b like he wants.

a = measurements(r).Centroid(1)

b = measurements(r).Centroid(2)

message = sprintf('The centroid of peak %d is at (%.2f, %.2f)\n', ...

r, measurements(r).Centroid(1), measurements(r).Centroid(2));

uiwait(msgbox(message));

end

The only difference is you named the array Po this time, and I assign variables a and b so they have the names you want. But essentially it's the same as before.

[PChoppala]

Yep, worked. Thanks a lot.

Just to wrap up, I am not able to figure out what number I should be subtracting from x and y when the mesh is [-40:0.2:40, -40:0.2:40], can you please help me out.

I tried a few but I am still wrong.

[Image Analyst]

Well just like you subtracted -51 from x when it was -50 to send it to +1, you would want to subtract -41 (i.e. add 41) from your x and y values so that the indices go from 1 to 81 instead of -40 to +40

[PChoppala]

Sorry, but that is not working out. I tried it. The values should be near s1(1,2) and s2(1,2)

I am not able to find the indices on Sa and Sb that correspond!!

Answer 2

[Walter Roberson]

Your previous question is still here and still active. I have merged your duplicate question into it. Please do not post duplicate questions: doing that just confuses and frustrates people.

[PChoppala]

Could not agree more. But what do you expect me to do, given the fact that the continued question was unanswered for 23 hours?

I had to assume the volunteer is either not interested or busy. But you cannot expect me to keep waiting and not attempt anymore!

Anyway thanks!

[Image Analyst]

Learning how to use the debugger to step through your code is always a good course of action.

[Walter Roberson]

Re-posting doesn't make the volunteers more interested or less busy.

If you wish to remind people about your question, you can do that by adding a comment to your question. In cases where you realized your question could be improved or did not have enough information, you can edit your question.