A severe problem is that Table generates an unpacked array, a thing that is often very annoying. When converting to an image, it is packed automatically (one can check that, e.g., with Developer`PackedArrayQ[ImageData[Image[dat]]]). And because many functions work faster on packed arrays than on unpacked ones, this increases the performance.
A vectorized implementation that exploits the very nature of the function (and is thus not very generalizable) is the following:
nx = ny = 1000;
dat = KroneckerProduct[
(0.1 + Cos[Subdivide[0., 2. Pi, ny]]),
Sin[Subdivide[0., 2. Pi, nx]]
];
0.003691
The maxima can be found by comparing each value to the maximum of its neighbors. The latter can be done efficiently with MaxFilter:
getindex = {Quotient[#, Dimensions[dat][[1]]] + 1, Mod[#, Dimensions[dat][[2]], 1]} &;
idx = getindex /@
Random`Private`PositionsOf[
Flatten@UnitStep[dat - MaxFilter[dat, {1, 1}]],
1
]; // AbsoluteTiming // First
idx
MatrixPlot@SparseArray[idx -> 1, Dimensions[dat]]
0.00108
{{1, 26}, {29, 1}, {30, 1}, {31, 1}, {32, 1}, {33, 1}, {34, 1}, {35, 1}, {36, 1}, {37, 1}, {38, 1}, {39, 1}, {40, 1}, {41, 1}, {42, 1}, {43, 1}, {44, 1}, {45, 1}, {46, 1}, {47, 1}, {48, 1}, {49, 1}, {50, 1}, {51, 1}, {51, 76}, {52, 1}, {53, 1}, {54, 1}, {55, 1}, {56, 1}, {57, 1}, {58, 1}, {59, 1}, {60, 1}, {61, 1}, {62, 1}, {63, 1}, {64, 1}, {65, 1}, {66, 1}, {67, 1}, {68, 1}, {69, 1}, {70, 1}, {71, 1}, {72, 1}, {73, 1}, {101, 26}}
This is about 10 times faster than MaxDetect>. But it does not distiguishdistinguish between local maxima and strict local maxima. However, a second sweep over the local maxima (which is much less expensive) could filter out the strict local maxima.
