<<Developer`
w1 = 600; w2 = 600; mask =
Sum[RotateRight[
Re[Exp[I 56 k] DiskMatrix[20, {w1, w2}]], {k 100, k 200}], {k,
2}] + Re[Exp[I 5] DiskMatrix[20, {w1, w2}]] +
Sum[RotateLeft[
Re[Exp[I 5 k] DiskMatrix[20, {w1, w2}]], {k 200, k}], {k,
1}]; Image[0.5 mask, Magnification -> 1] ;
ft = (Re[Fourier[mask]])^2;
ift = Image[
RotateRight[
ft, {w1/2, w2/2}]\[TensorProduct]ToPackedArray[{1.0, 0.3, 0.1}],
Magnification -> 1];
bresenham[p0_, p1_] :=
Module[{dx, dy, sx, sy, err, newp}, {dx, dy} =