ReadImage/: ReadImage[filename_, n_] := 
   ReadList[filename, Table[Number, {n}]]
DispIt/: DispIt[t_] := ListPlot3D[t, Lighting -> False, 
    Mesh -> False, AspectRatio -> 3, 
    ViewPoint -> {0, 0, -1000}, Axes -> False]
DispIt/: DispIt[t_, a_] := ListPlot3D[t, Lighting -> False, 
    Mesh -> False, AspectRatio -> a, 
    ViewPoint -> {0, 0, -1000}, Axes -> False]
EpsOut/: EpsOut[which_] := Display["!psfix -epsf > out.ps", 
    which]
SetOptions[ListPlot, PlotJoined -> True];

nptsx = 64;
tmp = 
        ReadImage["C:/simon/math/teaching/3c72/arridge.aim", 
          nptsx];
nptsy = Length[tmp];
imr = tmp;
For[i = 1, i ≤ nptsy, i++, 
        imr[[nptsy + 1 - i]] = 
          tmp[[i]]];

nptsy  =  
        nptsx;
ims = Table[Table[80 Random[], {nptsx}], {nptsy}];
ystart = Floor[(nptsy/4)];
xstart = Floor[(nptsx/4)];
yw = Floor[(nptsy/4)];
xw = Floor[(nptsx/4)];
ycen = Floor[((3 nptsy)/4)] - 2;
xcen = Floor[((3 nptsx)/4)] - 2;
For[m = 1, m ≤ nptsy, m++, 
        For[n = 1, 
          n ≤ 
            nptsx, n++, ims[[m, 
                n]] += 
              If[(((m - ycen))^2) + (((n - xcen))^2) < 0.8 (xw^2), 200, 
                0];]];
For[m = 1, m ≤ yw, m++, 
        For[n = 1, n ≤ xw, n++, 
          ims[[ystart + m, xstart + n]] += 
            100]];

im  =  
        ims;

{nptsx, nptsy}
ListDensityPlot[im, Mesh -> False];