I dont think it is possible to fully do this workflow strictly with array-based operations in IDL. However, you can eliminate at least one loop by using a 1D array. For example, the bottom right pixel of 100x100 image can be referenced in IDL as either image[99,99] or by image[9999], the latter being the single vector approach which may should run your code ~50% faster utilizing only one loop.

With that said, I was not able to easily get this working as it then becomes tricky about how you move the 3x3 window through a linear vector of values. This would be akin to moving three separate 3-value-wide bands through the array simultaneously.

Perhaps that helps at least get you started?

You may also find algorithms for this online. Just search moving max and moving mins of 2D arrays. That may give you an idea how to code this up in IDL.

While not completely vectorized, here is an alternative that does reduce the number of loops significantly by increasing the vectorized work done by IDL per loop. (reduces the number of iterations in the nested loops from 260100 to 9). The nested loops in the example code above took approximately 0.15 seconds to run, the loops below 0.0021 seconds, a speed up of over 50x.

Image = BINDGEN(512, 512) ; simulate the image for this example
ReturnImage = BYTARR(512, 512) ; define the data array

TempImage = MEDIAN(Image, 9) ; remove salt and pepper noise

boxsize = 3 ; size of entropy box (boxsize is always odd (1, 3, 5...))
if not boxsize mod 2 then boxsize += 1 ;Ensure boxsize is odd by adding 1 if the value is even
n = (boxsize -1) / 2

;set up some temporary arrays that will be used for finding the min and max values
;If we don't fill temp1 and temp2 with some original data then the min/max comparisons may fail.
;For instance, if we fill temp2 with all zeros, then the minimum value of the data will always be zero

temp1 = TempImage[n:(-1*n)-1,n:(-1*n)-1] ;The working arrays are smaller than the original image in each dimension by a value of 2 * n
temp2 = temp1
dat = bytarr(temp1.dim,/nozero)

;For each pixel there is set of boxsize^2 pixels to be checked for max and min values.
;We can take the entire set of data and shift the values boxsize^2 times, and for each iteration grab the highest and lowest values over the entire data set at once.
;For each pixel there are 8 surrounding pixels, plus the original, to be checked
;Here we shift the data 9 times with a set of vertial and horizontal shifts so that each of the 9 pixels of interest are compared
;and the min and max values saved
tic
for i = 0, boxsize - 1 do begin
for j = 0, boxsize - 1 do begin
dat[0,0] = TempImage[i:i-boxsize,j:j-boxsize] ;fill dat with image data that has been shifted
temp1 >= dat ;If there is a value temp1 that is less than the value in dat, replace the value in temp1 with the value in dat
temp2 <= dat ;If there is a value temp2 that is greater than the value in dat, replace the value in temp2 with the value in dat
endfor
endfor
ReturnImage[n,n] = temp1 - temp2 ;Put the data into ReturnImage shifted by n indices in the x and y. The preserves the zero values at the border from the original example
toc
ReturnImage = BYTSCL(ReturnImage)