To find the largest values, say in an 8-bit grayscale image, you can do
[i,j]=find(im==max(im(:)))
or for values which are nearly as large:
m = 10
[i,j]=find(im>=max(im(:))-m)
This will give you vectors of i and j indices. And to look at pixel values around such a peak, you can always isolate a local neighbourhood with
k = 1
im(i(k)-2:i(k)+2,j(k)-2:j(k)+2)
which will give you a 5x5 neighbourhood centred on the first pixel in your list.
If you'd prefer a circular neighbourhood, try something like
rad = 3
[x,y] = meshgrid(-rad:rad);
circ = x.^2+y.^2>rad^2
im(i(k)-rad:i(k)+rad,j(k)-rad:j(k)+rad).*circ
Is this the sort of thing you want?
To get your 3 x N array, suppose your image has R rows and C columns. That is:
[R,C] = size(im)
Then:
[X,Y] = meshgrid(1:R,1:C);
Z = X*Y
[reshape(X,Z,1),reshape(Y,Z,1),reshape(im,Z,1)]
There may be an easier way than this, but if there is I can't think of it right now.