pro blurring

  ; leggiamo l'immagine
fits_read,'D:\Ba\Ei1\phantom.fits',img

  ; dimensione dell'immagine
s=size(img)
print,'size',s
n_pix=s[2]

  ; definizione degli array
blur=make_array(n_pix,n_pix,/double,value=0.)
psf=make_array(n_pix,n_pix,/double,value=0.)

  ;creazione psf (si puņ anche leggere da file)
for i =n_pix/2,n_pix/2+19 do begin
psf[i,n_pix/2]=1./20.
end
print,'min psf',min(psf)
print,'massimo psf', max(psf)
print,'somma sui pixel' , total(psf)

  ;blurring, creazione della funzione di trasferimento
tf=fft(shift(psf,n_pix/2,n_pix/2),-1,/double)*n_pix*n_pix
print,'max TF',max(abs(tf))
print,'min TF', min(abs(tf))
print,'cond', max(abs(tf))/ min(abs(tf))

  ;visualizzazione funzione di trasferimento
window,0,xsize=n_pix,ysize=n_pix
tvscl,alog10(abs(tf))

  ;traformata di Fourier dell'immagine
tfimg=fft(img,-1,/double)*n_pix*n_pix

  ;convoluzione tra immagine e funzione di trasferimento
blur=(fft((tf*tfimg),1,/double))/n_pix/n_pix

  ;scrittura immagine blurrata
fits_write,'blur_nonoise.fts',float(blur)

  ;visualizzazione del risultato
window,1,xsize=n_pix,ysize=n_pix
tvscl,float(blur)

  ;aggiunta del noise
noise=randomu(1.,n_pix,n_pix,normal=0.2*max(img))
blurnoise=blur+noise

  ;scrittura dell'immagine blurrata
fits_write,'blur.fts',float(blurnoise)

  ;visualizzazione del risualtato
window,2,xsize=n_pix,ysize=n_pix
tvscl,float(blurnoise)-float(blur)

  end