finding multiple maxima in a Gaussian distribution

3 次查看（过去 30 天）
PChoppala
PChoppala 2012-3-9
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
  1 个评论
Walter Roberson
Walter Roberson 2012-3-10
http://www.mathworks.com/matlabcentral/answers/29922-why-your-question-is-not-urgent-or-an-emergency

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Image Analyst
Image Analyst 2012-3-9
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.
  7 个评论
显示 5更早的评论
Image Analyst
Image Analyst 2012-3-12
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
PChoppala 2012-3-12
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!!

请先登录，再进行评论。

Walter Roberson
Walter Roberson 2012-3-10
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.
  3 个评论
显示 1更早的评论
Image Analyst
Image Analyst 2012-3-10
Learning how to use the debugger to step through your code is always a good course of action.
Walter Roberson
Walter Roberson 2012-3-10
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Logical 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by