Fourier transform possible on non-rectangular part of an image

Question

Dear Signal Processing readers,

I want to introduce 'noise' into parts of images.

Until now, I worked with rectangluar images for a similar purpose and did the following, using (inverse) Fourier transform (2D) with MATLAB:

fft = fft2(image);    % run fourier transform on image
amplitude = abs(fft); % extract power spectrum>
newPhase = ( (rand(size(amplitude)) * 2) - 1) * pi;          %generate random phase
newImage = real(ifft2(amplitude. *  exp(1i * newPhase)));    %use both to generate noised image

This randomizes the phase all over the image, preserving the power spectrum and should guarantee the same spatial frequency spectra in input and output.

What I want to do now, is applying this effect to only a part of an image, which is not rectangular (its the "inner" parts of a face and will fit in an ellipse). So I want to "cut out" this ellipse from an image, turn it to "phase noise" (with power spectrum/spatial frequencies constant) as above, and but it back into the "hole" I cut out of the source image before.

My "prime solution" of course would be that FFT/iFFT could operate on ellipses as well, so I could keep using my method shown above. After some amount of google-ing I am quite pessimistic about that option - is there any hope?

I found the following traces to possibly promising approaches, but as I am not very familiar with the math behind FFT, I'd appreciate any hint/explanation (why) which of these ideas go into a good direction:

• Somehow "stretch" the ellipse to a rectangle shape, apply my method as is, strech the resulting rectangle back to the input ellipse shape. Here I fear that the rectangle/ellipse-transformation might destroy the matched power spectrum. Do you know good way of transforming without loosing this information? (This led me to this idea)

Two more hints came from this post:

• Using a 3D-fourier transform, in which (as far as I got it) the 3rd dimension gives information about which part of the image shall be used. Any idea / experience how to implement this?

• Using "Elliptic Discrete Fourier Transform" - but I do not understand what (and how) it does and how it can be implemented.

• In this forum, someone suggested to a related question some sort of window tapering (which, as far as I got it, is made to "filter" the influence that only comes into the image because the ellipse of interest is on a black background) that I did not understand as well, plus i am not sure if it could help me because it did not sound like that question would require an ellipse as output.

Thanks in advance for any suggestion!

Cheers. Jan

Answer 1

Even when the image part is rectangular, the FFT result will be contaminated by boundary effects. The FFT of the part is the Fourier transform of an infinite tiling of the rectangular part. (An elliptical or circular part cannot be tiled. A hexagonal tiling would get pretty close though.) If there is a discontinuity going from one tile to another, this will create high-frequency components that are not descriptive of the insides of the part. Randomizing the phases will spatially translate these components so that they appear also inside the part, which is probably not what you want. Mirroring at the boundaries does not eliminate discontinuities of the derivative, and requires doubling of the FFT dimensions.

This is where a smoothly tapering 2-d window function can help. For a circular part you might like the Tukey window with the distance from the part center as its argument.

For an ellipse, stretch the window function accordingly. Multiply the data by the window function to select the area you want to analyze, and FFT the result. Then do the phase randomization. Mask the noise by multiplying it with an ellipitical disc-like mask that has value 0 outside the disc and value 1 inside the disc. You can also use the window function as the mask to get a fuzzy boundary. I did not find yet a perfect way to normalize the noise, but matching the mean and the standard deviation both weighted by the mask works quite OK:

Multiply the input image by the complementary mask: 1 - mask, and sum the two together. Octave source:

pkg load signal;
r = 120; # Circle radius
n = 512; # Image width and height
tukeyr = 0.5; # Window transition relative size
lena = double(imread("lena.gif"))/255;
tukeyn = round(r/(1-tukeyr));
x = repmat([1:n], n, 1);
y = repmat([1:n]', 1, n);
dist = sqrt((x-n/2).^2 + (y-n/2).^2);
tukey = tukeywin(tukeyn*2, tukeyr)(tukeyn+1:tukeyn*2);
win = reshape(interp1(1:tukeyn, tukey, reshape(dist, n*n, 1), 0), n, n);
randphase = 2*pi*rand(n, n);
noisefft = abs(fft2(lena.*win)).*exp(j*randphase);
noise = real(ifft2(noisefft));

#Circular mask
weight = r-dist;
weight = (weight >= 1)*1 + (weight > 0 & weight < 1).*weight;
#Or use window as mask
#weight = win;

weightsum = sum(sum(weight));
noiseweightedmean = sum(sum(noise.*weight))/weightsum;
noiseweightedstdev = sqrt(sum(sum((noise-noiseweightedmean).^2.*weight))/weightsum);
lenaweightedmean = sum(sum(lena.*weight))/weightsum;
lenaweightedstdev = sqrt(sum(sum((lena-lenaweightedmean).^2.*weight))/weightsum);
normalizednoise = ((noise - noiseweightedmean) / noiseweightedstdev * lenaweightedstdev) + lenaweightedmean;
imshow(lena.*(1-weight).+normalizednoise.*weight);

P.S. I found that highpass filtering the copy of the Lena image used for noise generation, by setting some of its lowest-frequency bins to zero, will reduce the error in the magnitudes of the fft of the final image compared to fft of the source image. I think the reason is that low-frequency amplitude fluctuations will be reduced. Within the circle (ellipse) those appear in practice as zero frequency, and in their absence less normalization is needed. So that's one thing to investigate.