The values are scaled so that when you view the image it's not mostly black.
I tested that the scaling works with straight integer truncation by checking that:
[A] = imread( 'myfile.pgm', 'pgm' );
p = sort(unique(A(:));
q = uint16((0:1250) * 65535 / 1251)';
all(p == q) % returns 1
So, you can restore the image like this:
map = arrayfun( @(x) uint16(x * 1251 / 65536), 0:65535 );
B = arrayfun( @(x) map(x+1), A );