Thread Subject: Embedded Image Resizing using FFT

I'm currently looking for a way to rescale/resize images using Embedded Matlab.
I've read it should be possible using the Frequency domain. And since the fft2 and ifft2 functions are supported by the Embedded Matlab subset I thought I'd go this way. If you have any other suggestions they are also welcome

However I haven't found any actual information on the subject or how to proceed. Can you offer any help or examples?

I'd like to achieve something like:

factor = 0.5;
image = imread('test.bmp');
%Each of the RGB planes is 0 to 255
imageR(:,:,1) = fftResize(image(:,:,1), factor);
imageR(:,:,2) = fftResize(image(:,:,2), factor);
imageR(:,:,3) = fftResize(image(:,:,3), factor);
imwrite(imageR, 'test_resized.bmp', 'bmp');

function img_re = fftResize(img, factor)
img = im2double(img);
fft = fft2(img);

%Part where I need help%

img_re = ifft2(img_resized);
img_re = uint8(floor(img_re));
end

"Kevin Voet" wrote:
> I'm currently looking for a way to rescale/resize images using Embedded Matlab.
> I've read it should be possible using the Frequency domain. And since the fft2 and ifft2 functions are supported by the Embedded Matlab subset I thought I'd go this way.

This has often struck me as an excellent way to resize images. You crop the image in the frequency domain, getting rid of the high frequencies that cannot be seen in the smaller image, then transform back to the spatial domain. It should produce the 'ideal' resized image, i.e. one that retains the maximum amount of information possible in the smaller size.

I've tried to do this myself, cropping off various different parts of the frequency image, but never got back anything that looked right. The question is, can you just cut out bits of the frequency image, then transform back, or is there some reason why this cannot work? If it is possible, which parts do you cut out?

"Kevin Voet" <kevin.voet@gmail.com> wrote in message <gg4co5$2cu$1@fred.mathworks.com>...
> I'm currently looking for a way to rescale/resize images using Embedded Matlab.
> I've read it should be possible using the Frequency domain. And since the fft2 and ifft2 functions are supported by the Embedded Matlab subset I thought I'd go this way. If you have any other suggestions they are also welcome
>
> However I haven't found any actual information on the subject or how to proceed. Can you offer any help or examples?
>
> I'd like to achieve something like:
>
> factor = 0.5;
> image = imread('test.bmp');
> %Each of the RGB planes is 0 to 255
> imageR(:,:,1) = fftResize(image(:,:,1), factor);
> imageR(:,:,2) = fftResize(image(:,:,2), factor);
> imageR(:,:,3) = fftResize(image(:,:,3), factor);
> imwrite(imageR, 'test_resized.bmp', 'bmp');
>
> function img_re = fftResize(img, factor)
> img = im2double(img);
> fft = fft2(img);
>
> %Part where I need help%
>
> img_re = ifft2(img_resized);
> img_re = uint8(floor(img_re));
> end



Presumably the method exploits the property that f(a*t) Fourier transforms into
F(w/a)/|a|.

However, I can't see why this would be better than resizing in the non-Fourier domain. You end up doing this same operation, pretty much, in either domain.

"Matt" <mjacobson.removethis@xorantech.com> wrote in message <gg4jc2$bck$1@fred.mathworks.com>...
> "Kevin Voet" <kevin.voet@gmail.com> wrote in message <gg4co5$2cu$1@fred.mathworks.com>...
> > I'm currently looking for a way to rescale/resize images using Embedded Matlab.
> > I've read it should be possible using the Frequency domain. And since the fft2 and ifft2 functions are supported by the Embedded Matlab subset I thought I'd go this way. If you have any other suggestions they are also welcome
> >
> > However I haven't found any actual information on the subject or how to proceed. Can you offer any help or examples?
> >
> > I'd like to achieve something like:
> >
> > factor = 0.5;
> > image = imread('test.bmp');
> > %Each of the RGB planes is 0 to 255
> > imageR(:,:,1) = fftResize(image(:,:,1), factor);
> > imageR(:,:,2) = fftResize(image(:,:,2), factor);
> > imageR(:,:,3) = fftResize(image(:,:,3), factor);
> > imwrite(imageR, 'test_resized.bmp', 'bmp');
> >
> > function img_re = fftResize(img, factor)
> > img = im2double(img);
> > fft = fft2(img);
> >
> > %Part where I need help%
> >
> > img_re = ifft2(img_resized);
> > img_re = uint8(floor(img_re));
> > end
>
>
>
> Presumably the method exploits the property that f(a*t) Fourier transforms into
> F(w/a)/|a|.
>
> However, I can't see why this would be better than resizing in the non-Fourier domain. You end up doing this same operation, pretty much, in either domain.
>

Performance for one, instead of using bicubic or bilineair transforms on the image you should be able to resize the image using i.e. a lowpass filter to make the image smaller and add zero-padding to make the image larger (or so I've read).

What technique would you suggest to resize an image in the Embedded Matlab Subset?

PS: Yes, this will be implemented in hardware eventually as part of a larger system

> > Presumably the method exploits the property that f(a*t) Fourier transforms into
> > F(w/a)/|a|.
> >
> > However, I can't see why this would be better than resizing in the non-Fourier domain. You end up doing this same operation, pretty much, in either domain.
> >
>
> Performance for one, instead of using bicubic or bilineair transforms on the image you should be able to resize the image using i.e. a lowpass filter to make the image smaller and add zero-padding to make the image larger (or so I've read).

I don't really follow that. Lowpass filtering an image will make it smoother, not smaller.

Zero-padding in frequency space will cause the image to be upsampled by sinc-interpolation. You could then use these upsamplings to resize your image.
 
However, it would give you many more interpolated samples than you need and than you would use if you just re-interpolated the image directly in image space using say, interp2(). I can't see how that ends up being more efficient.


> What technique would you suggest to resize an image in the Embedded Matlab Subset?

I don't see why you wouldn't just write your own fast interp2.c function and use that...

"Matt" <mjacobson.removethis@xorantech.com> wrote in message <gg4mdj$sjg$1@fred.mathworks.com>...
>
> > > Presumably the method exploits the property that f(a*t) Fourier transforms into
> > > F(w/a)/|a|.
> > >
> > > However, I can't see why this would be better than resizing in the non-Fourier domain. You end up doing this same operation, pretty much, in either domain.
> > >
> >
> > Performance for one, instead of using bicubic or bilineair transforms on the image you should be able to resize the image using i.e. a lowpass filter to make the image smaller and add zero-padding to make the image larger (or so I've read).
>
> I don't really follow that. Lowpass filtering an image will make it smoother, not smaller.
>
> Zero-padding in frequency space will cause the image to be upsampled by sinc-interpolation. You could then use these upsamplings to resize your image.
>
> However, it would give you many more interpolated samples than you need and than you would use if you just re-interpolated the image directly in image space using say, interp2(). I can't see how that ends up being more efficient.
>
>
> > What technique would you suggest to resize an image in the Embedded Matlab Subset?
>
> I don't see why you wouldn't just write your own fast interp2.c function and use that...

This was my source for the lowpass comment:
http://www.watermarkingworld.org/WMMLArchive/0302/msg00055.html

I haven't got a clue how to use the interp2() function to resize an image and certainly not in C.

> I haven't got a clue how to use the interp2() function to resize an image and certainly not in C.

I'm pretty sure the bilinear/bicubic transformation operations you were refering to involved re-interpolating the image

> This was my source for the lowpass comment:
> http://www.watermarkingworld.org/WMMLArchive/0302/msg00055.html

OK.
 Well, he isn't talking about making an image smaller by low-pass filtering it. He's talking about making it smaller by shortening the FFT array of the signal, i.e., you actually throwing away the high frequency components of the array, not set them to zero.
This has the effect of decreasing the sampling interval of your image (in effect shortening it). You get the converse effect by zero-padding the FFT array.

However, I can't see how any of this is more efficient than just interpolating your image to obtain more or less densely spaced samples
Doing geometric transformations this way on an NxN image is an O(N^2) operation.
Taking FFTs of an NxN image is an O( N^2*log(N) ) operation and so is more intensive.

This is made worse by

(1) If you zero-pad the FFT array, so as to dilate the image, then N increases, requiring even more computation.

(2) FFTs involve complex arithmetic, i.e. twice as much numeric data to manipulate as in the image domain where only real arithmetic is required.

"Matt" wrote:
> Doing geometric transformations this way on an NxN image is an O(N^2) operation.
> Taking FFTs of an NxN image is an O( N^2*log(N) ) operation and so is more intensive.

Have you tried this? There are scale factors to consider. In practice fft2 isn't too bad:
>> A = rand(1000,1000,3);
>> tic; imresize(A, 1/2); toc
Elapsed time is 0.606730 seconds.
>> tic; fft2(A); toc
Elapsed time is 0.519417 seconds.

Forgetting speed for a minute, I think a benefit of the frequency domain approach---simply cropping the unwanted high frequencies out---is that it gives you the best possible resizing (if someone can work out whether it's possible and how) because it keeps all the frequencies which will be visible in the smaller (assuming this is for downsampling) image. You don't need to work out any optimal smoothing filter for the spatial domain.

> Have you tried this? There are scale factors to consider. In practice fft2 isn't too bad:
> >> A = rand(1000,1000,3);
> >> tic; imresize(A, 1/2); toc
> Elapsed time is 0.606730 seconds.
> >> tic; fft2(A); toc
> Elapsed time is 0.519417 seconds.
#############################

Here's what I tried.

First the data prep.
>>A=rand(1000,1000); A=single(A);
>>[xx,yy]=ndgrid(1:.5:1000,1:.5:1000); xx=single(xx); yy=single(yy);

These are the operations you would do to upsample by a factor of 2 using the frequency based method.

>>tic;B=fft2(A); C=ifft2(B,2000,2000);toc
Elapsed time is 1.184620 seconds.

Conversely, these are the operations you would do to upsample by straight linear interpolation. I'm using my own C-coded linear interpolation routine because MATLAB's interp2 is handicapped (it is implemented as an mfile).

>> >> tic; z=linterpcpp(A,xx,yy);toc
Elapsed time is 0.142840 seconds.


#####################
> Forgetting speed for a minute, I think a benefit of the frequency domain approach---simply cropping the unwanted high frequencies out---is that it gives you the best possible resizing (if someone can work out whether it's possible and how) because it keeps all the frequencies which will be visible in the smaller (assuming this is for downsampling) image. You don't need to work out any optimal smoothing filter for the spatial domain.
#######################


I don't think so. As far as I can tell, cropping frequencies (say by an integer factor 2) is equivalent to decimating the image by the same factor f(n,m)-->f(2n,2m).
This is known to give bad results.

> I don't think so. As far as I can tell, cropping frequencies (say by an integer factor 2) is equivalent to decimating the image by the same factor f(n,m)-->f(2n,2m).
> This is known to give bad results.

We don't agree so we should do it and see.

I've found something that seems to work:

>> A = imread('cameraman.tif');
>> subplot(121); imshow(A);
>> B = fftshift(fft2(im2double(A)));
>> B = ifft2(ifftshift(B(71:end-70,71:end-70)));
>> subplot(122); imagesc(real(B)); colormap(gray(255)); axis image off

It reduces the image by 140 pixels in each dimension. You can enlarge by padding with zeros. You can avoid the shifts by working out what they do and cropping the correct regions directly.

You'll note that I take the real component - the output is complex, which is a bad sign. I also use imagesc to display the output image, as it isn't scaled in the [0, 1] range any more, again a bad sign.

Looking at image quality, it is clear cropping high frequencies doesn't simply subsample the image (i.e. take every ith pixel). However, there are some high frequency artifacts. Overall it seems like it won't work well.

Changing the last line of my code to:
>> subplot(122); imagesc(abs(B)); colormap(gray(255)); axis image off
gives a slightly more pleasing result, though still not perfect.

"Oliver Woodford" <o.j.woodford.98@cantab.net> wrote in message <gg64up$sn8$1@fred.mathworks.com>...
> > I don't think so. As far as I can tell, cropping frequencies (say by an integer factor 2) is equivalent to decimating the image by the same factor f(n,m)-->f(2n,2m).
> > This is known to give bad results.
>
> We don't agree so we should do it and see.
>
> I've found something that seems to work:
>
> >> A = imread('cameraman.tif');
> >> subplot(121); imshow(A);
> >> B = fftshift(fft2(im2double(A)));
> >> B = ifft2(ifftshift(B(71:end-70,71:end-70)));
> >> subplot(122); imagesc(real(B)); colormap(gray(255)); axis image off
>
> It reduces the image by 140 pixels in each dimension. You can enlarge by padding with zeros. You can avoid the shifts by working out what they do and cropping the correct regions directly.

> You'll note that I take the real component - the output is complex, which is a bad sign. I also use imagesc to display the output image, as it isn't scaled in the [0, 1] range any more, again a bad sign.
>

I don't have the Image Processing Toolbox, so I can't try this myself.

However, I take it that the dimensions of A are even? If so, the reason you are getting a complex output is because you are not truncating frequencies symmetrically. To get symmetry in an even-sized array, you have to truncate the upper range of the array by 1 less than the lower.

 B = ifft2(ifftshift(B(71:end-69,71:end-69)));

My guess is that if you fix this, the scaling issue will go away too, because you'll no longer have to throw away imaginary components.


> Looking at image quality, it is clear cropping high frequencies doesn't simply subsample the image (i.e. take every ith pixel).

Yes, I agree I was wrong about that. It should, however, end up being nearly equivalent to subsampling a filtered version of the image.

When you throw away high frequencies, it is equivalent to first doing low pass rect-windowing of the spectrum, or equivalently of convolving the image with an (aliased) sinc kernel. Perhaps this is what you were refering to earlier by frequency-preservation...

When you then throw away the components outside the rect-window, that is the frequency-domain dual of subsampling the sinc filtered image.

So, anyway, I'm becoming more convinced of the image-quality potential of the technique. Computationally, however, I still wonder if you could find more efficient image-domain equivalents. Sinc-filtering for example can be approximated very cheaply in the non-Fourier domain with cubic B-splines.

Code which doesn't use toolboxes:

>> A = sum(imread('street1.jpg'), 3);
>> subplot(121); imagesc(A); colormap(gray(255)); axis image off
>> B = fftshift(fft2(A));
>> B = ifft2(ifftshift(B(116:end-114,171:end-169)));
>> subplot(122); imagesc(B); colormap(gray(255)); axis image off

Matt, your point on symmetry was correct and removes the imaginary component. However, the intensity scale still changes. Quality is ok, but there are still some high frequency artifacts - most visible at the top of this image, but more visible still on the cameraman image.

> Matt, your point on symmetry was correct and removes the imaginary component. However, the intensity scale still changes. Quality is ok, but there are still some high frequency artifacts - most visible at the top of this image, but more visible still on the cameraman image.

Sinc-ripple, I suspect. As I was saying, the technique is a combintion of sinc-interpolation and subsampling operations...

"Matt" wrote:
> As I was saying, the technique is a combintion of sinc-interpolation and subsampling operations...

I agree. However, sinc filters have long tails, so using them can be slow, and slower still the more you want to reduce your image (as the filters get larger). This method gets faster (the ifft2 bit, anyway) the more you downsize. Shame about the ripples...

Anyway, Kevin, I hope we've helped.

"Oliver Woodford" <o.j.woodford.98@cantab.net> wrote in message <gg6qrh$afk$1@fred.mathworks.com>...
> "Matt" wrote:
> > As I was saying, the technique is a combintion of sinc-interpolation and subsampling operations...
>
> I agree. However, sinc filters have long tails, so using them can be slow, and slower still the more you want to reduce your image (as the filters get larger).

That's true, however, it is known that linear combinations of sincs can be well approximated by linear combinations of cubic B-splines, which are of small support and faster to manipulate. This is one reason why cubic B-splines are so popular -- they provide a computationally efficient way of aproximating sinc interpolation.

"Oliver Woodford" <o.j.woodford.98@cantab.net> wrote in message <gg6qrh$afk$1@fred.mathworks.com>...
> "Matt" wrote:
> > As I was saying, the technique is a combintion of sinc-interpolation and subsampling operations...
>
> I agree. However, sinc filters have long tails, so using them can be slow, and slower still the more you want to reduce your image (as the filters get larger).

That's true, however, it is known that linear combinations of sincs can be well approximated by linear combinations of cubic B-splines, which are of small support and faster to manipulate. This is one reason why cubic B-splines are so popular -- they provide a computationally efficient way of aproximating sinc interpolation.

"Matt" <mjacobson.removethis@xorantech.com> wrote in message <gg71ql$57t$1@fred.mathworks.com>...
> "Oliver Woodford" <o.j.woodford.98@cantab.net> wrote in message <gg6qrh$afk$1@fred.mathworks.com>...
> > "Matt" wrote:
> > > As I was saying, the technique is a combintion of sinc-interpolation and subsampling operations...
> >
> > I agree. However, sinc filters have long tails, so using them can be slow, and slower still the more you want to reduce your image (as the filters get larger).
>
> That's true, however, it is known that linear combinations of sincs can be well approximated by linear combinations of cubic B-splines, which are of small support and faster to manipulate. This is one reason why cubic B-splines are so popular -- they provide a computationally efficient way of aproximating sinc interpolation.

Would blockprocessing make a difference in performance?
e.g. img_dft = blkproc(img,[8 8],@fft2);

Too bad block processing isn't supported in the Embedded MATLAB subset.
And also not sure how I'd add a scaling factor.

If I were to start working on a cubic B-splines algorithm would it be possible to create one in the Embedded MATLAB subset?

> If I were to start working on a cubic B-splines algorithm would it be possible to create one in the Embedded MATLAB subset?

Should be possible. Most B-Spline manipulation can be posed as simple matrix operations. You should also be able to find 3rd party B-Spline manipulation code fairly easily online.