function h = imagehistogram(imago)
  if ( isgray(imago) )
    for j = 0:255
      h(j+1) = numel(imago( imago == j ));
    endfor
  else
    error("histogram on gray img only");
  endif
endfunction

% test
im = jpgread("Lenna100.jpg");
img = rgb2gray(im);
h = imagehistogram(img);
% let's try to show the histogram
bar(h);
pause;

% in order to obtain the requested filtering, we
% can use median directly on the img, and then
% use that value, this way:
m = median(reshape(img, 1, numel(img)));
disp(m);
ibw = img;
ibw( img > m ) = 255;
ibw( img <= m ) = 0;
jpgwrite("lennamed_.jpg", ibw, 100);
% which disagree (128) with the m computed with histog_med (130).
% If we compute it this way:
% m = sort(reshape(img, 1, numel(img)))(ceil(numel(img)/2));
% we obtain 130... but builtin median works as expected, since
% N (number of pixel of Lenna) is even, not odd.

% but let's use our histogram h instead
function m = histog_med(histog)
  from = 0; to = 255;
  left = histog(from + 1); right = histog(to+1);
  while ( from != to )
    if ( left < right )
      from++; left += histog(from+1);
    else
      to--; right += histog(to+1);
    endif
  endwhile
  m = from;
endfunction

m = histog_med(h);
disp(m);
ibw( img > m ) = 255;
ibw( img <= m ) = 0;
jpgwrite("lennamed.jpg", ibw, 100);