get rid of loop and large speed up in octave

Moved away from "for" loops, by using matrix notation.
Code is faster now and more importantly faster at least in octave
which is bad with loops.

I used to have in octave 17 seconds execution time now it is only 0.3
seconds. Matlab was reasonably fast even before.

1 files changed, 12 insertions, 13 deletions

diff --git a/stripeeraser.m b/stripeeraser.m
index 3198b1f..c0d48a6 100644
--- a/stripeeraser.m
+++ b/stripeeraser.m
@@ -1,4 +1,4 @@
-function [img_cleaned,img_fourier]=stripeeraser(img,radius,threshold_power)
+function [img_cleaned,img_fourier,img_mask]=stripeeraser(img,radius,threshold_power)
 % Cleans an image from the stripes produced by interference fringes.
 % Such parallel fringes appears like a strong peak in the Fourier transformed image.
 % so it easy to clean them, especially if their spatial frequency is high.
@@ -11,16 +11,15 @@ img_fourier=fftshift(fft2(img)); % move to Fourier space
 max1=max(abs(img_fourier(:)));
 threshold=max1*(10^(-threshold_power));
 [Ny, Nx] = size (img); % image size
-for x=1:Nx
-    for y=1:Ny
-        d=sqrt((x-Nx/2)^2+(y-Ny/2)^2);
-	if d>radius && abs(img_fourier(y,x))>threshold
-		% For high enough frequency components
-		% suppress this frequency components which are mostly stripes.
-		% But leave untouched small sharp features of the beam 
-		% since they are spectrally weak
-		img_fourier(y,x)=img_fourier(y,x)*exp(-(d/radius)^2);
-        end;
-    end;
-end;
+[x,y]=meshgrid(1:Nx,1:Ny); % matrix xy coordinates
+
+d=sqrt((x-(Nx+1)/2).^2+(y-(Ny+1)/2).^2); % distance from center of the image
+% For high enough frequency components
+% suppress this frequency components which are mostly stripes.
+% But leave untouched small sharp features of the beam 
+% since they are spectrally weak
+img_mask = d>radius & abs(img_fourier)>threshold;
+
+img_fourier(img_mask)=img_fourier(img_mask).*exp(-(d(img_mask)/radius).^2);
+
 img_cleaned = abs( ifft2(ifftshift(img_fourier)) );