User:Loisel

From Wikipedia, the free encyclopedia

I am Sébastien Loisel and I work at Temple University.

[edit] MATLAB code

Image:Blurry-star.png Image:Blurry-star-fourier.png

If you put power=0 in the code below, the blurry star image looks like complete noise, but somehow it still is able to compute the star radius. It's uncanny.

Here is the code..

function [radius,star,starft,blurrystar,blurryft,thresholded,thresholded1]=make_star(n,r0,power,threshold);
% try radius=make_star(150,15,6,25)
coords=(-n/2):(n/2-1);
% create a random filter that tapers off near the edge
filter=random('unif',0,1,n,n);
dampen=(cos((coords)*2*pi/n)+1)/2;
dampen2=(ones(n,1)*dampen).^power.*(dampen'*ones(1,n)).^power;
filter3=filter.*dampen2;
filter4=real(filter3/sum(sum(filter3)));
% create an off-center star
x0=random('unif',-n/20,n/20);
y0=random('unif',-n/20,n/20);
xc=ones(n,1)*coords;
yc=(coords)'*ones(1,n);
star=(((xc-x0).^2+(yc-y0).^2)<=(r0^2))+0;
% compute a blurry star
blurrystar=conv2(star,filter4,'same');
% compute the FFT of the blurry star
blurryft2=fft2(blurrystar);
% shift fft origin to center of matrix
modcoords=mod(coords,n)+1;
blurryft(modcoords,modcoords)=blurryft2;
starft2=fft2(star);
starft(modcoords,modcoords)=starft2;

%threshold and compute radius of star
thresholded=((abs(starft)<threshold)+0);
radiuses=sqrt(xc.^2+yc.^2);
radius0=min(min(radiuses.*thresholded+(1-thresholded)*10000));
thresholded1=thresholded.*((radiuses<(radius0+2))+0);
radius1=sum(sum(radiuses.*thresholded1))/sum(sum(thresholded1));
J1zero=3.83170597;
radius=n/pi/2*J1zero/radius1;
r1=n/pi/2*J1zero/r0;
close all;
img=star; [m,n]=size(img);
figure('Units','pixels','Position',[100 100 2*n 2*m],'Name','Star');
imagesc(img); colormap(gray); set(gca,'Position',[0 0 1 1]);
line([n/2+1,n/2+1],[1,m],'Linestyle',':','Color','r');
line([1,n],[m/2+1,m/2+1],'Linestyle',':','Color','r');
img=blurrystar; [m,n]=size(img);
figure('Units','pixels','Position',[400 100 2*n 2*m],'Name','Blurry star');
imagesc(img); colormap(gray); set(gca,'Position',[0 0 1 1]);
line([n/2+1,n/2+1],[1,m],'Linestyle',':','Color','r');
line([1,n],[m/2+1,m/2+1],'Linestyle',':','Color','r');
img=log(abs(starft)); [m,n]=size(img);
figure('Units','pixels','Position',[100 400 2*n 2*m],'Name','Star Fourier transform');
imagesc(img); colormap(gray); set(gca,'Position',[0 0 1 1]);
line([n/2+1,n/2+1],[1,m],'Linestyle',':','Color','r');
line([1,n],[m/2+1,m/2+1],'Linestyle',':','Color','r');
line(r1*cos(pi*(0:.1:2))+m/2+1,r1*sin(pi*(0:.1:2))+n/2+1,'Linestyle',':','Color','b');
img=log(abs(blurryft)); [m,n]=size(img);
figure('Units','pixels','Position',[400 400 2*n 2*m],'Name','Blurry star Fourier transform');
imagesc(img); colormap(gray); set(gca,'Position',[0 0 1 1]);
line([n/2+1,n/2+1],[1,m],'Linestyle',':','Color','r');
line([1,n],[m/2+1,m/2+1],'Linestyle',':','Color','r');
line(radius1*cos(pi*(0:.1:2))+m/2+1,radius1*sin(pi*(0:.1:2))+n/2+1,'Linestyle',':','Color','b');

This code is public domain, released by its author, Wikipedia user Loisel.