diff --git a/Examples/Test_SLIT.py b/Examples/Test_SLIT.py index 34112e2..4c40004 100644 --- a/Examples/Test_SLIT.py +++ b/Examples/Test_SLIT.py @@ -9,6 +9,13 @@ import warnings warnings.simplefilter("ignore") +#Example of a run of the SLIT algorithm on simulated images. +#Here the first part of the file shows how simulations are generated. +#For users intereseted only in seeing the code run, have a look at the running SLIT section. +#The command line that builds the Fkappa operator is also of outmost importance. + + +###############################Simulation############################### def SIE(x0,y0,n1,n2,b,beta,q,xc,theta): kappa = np.zeros((n1,n2)) x,y = np.where(kappa == 0) @@ -25,78 +32,77 @@ def SIE(x0,y0,n1,n2,b,beta,q,xc,theta): count += 1 return kappa - +#Source light profile newsource = pf.open('../Files/source.fits')[0].data ##N1,N2 are the numbers of pixels in the image plane. nt1= 100 nt2 = 100 +#Size ratio of the source to image number of pixels size = 1 -na1 = nt1 -na2 = nt2 +#PSF +PSF0 = pf.open('../Files/PSF.fits')[0].data +PSF = PSF0[1:,1:] +PSFconj = np.real(np.fft.ifft2(np.conjugate(np.fft.fft2(PSF0[:-1,:-1])))) +PSFconj=PSFconj/np.sum(PSFconj) +PSF = PSF/np.sum(PSF) + ## Lens mass distribution. b = 1.53/0.05 xc = 0.95 q = 0.71 betata = 2.1 thetata = 25.2 - -kappa = SIE(na1/2.+50,na2/2.+50,na1+100,na2+100,b,betata,q,xc, thetata) -theta = SLIT.Lens.F(kappa, nt1,nt2, size,na1/2.,na2/2.) - -x00 =3744 -y00 =2028 - - -PSF0 = pf.open('../Files/PSF.fits')[0].data -PSF = PSF0[1:,1:] -PSFconj = np.real(np.fft.ifft2(np.conjugate(np.fft.fft2(PSF0[:-1,:-1])))) -PSFconj=PSFconj/np.sum(PSFconj) -PSF = PSF/np.sum(PSF) +kappa = SIE(nt1/2.+50,nt2/2.+50,nt1+100,nt2+100,b,betata,q,xc, thetata) +#Mapping between lens and source IMPORTANT +Fkappa = SLIT.Lens.F(kappa, nt1,nt2, size,nt1/2.,nt2/2.) -I2 = SLIT.Lens.source_to_image(newsource, nt1 ,nt2 , theta) -PSF0 = np.abs(PSF0/np.sum(PSF0)) -PSF = np.abs(PSF/np.sum(PSF)) +#Generation of lensed source +I2 = SLIT.Lens.source_to_image(newsource, nt1 ,nt2 , Fkappa) +#Noise levels SNR = 500 sigma = np.sqrt(np.sum(I2**2)/SNR/(nt1*nt2*size**2)) - -I3 = np.copy(I2) +#Convolution by the PSF and generation of the final image I2 = scp.fftconvolve(I2, PSF, mode = 'same') + +#Final simulated image Image = I2+np.random.randn(nt1,nt2)*sigma + +################################Running SLIT############################ +#Parameters kmax = 5 niter =50 -pos_real = 1 -pos_wave = 0 -weight = [0] +levels = [0] +#Comment the following to have the level estimation routine run (takes more time) +levels = pf.open('../Files/Noise_levels_SLIT.fits')[0].data + +#Start clock start = time.clock() -sourcesl, Imsl = SLIT.SLIT(Image, theta, kmax, niter, size, mymy = 0, - pos = pos_real, posit = pos_wave, decrease = 0, soft =1, - weight = weight, PSF = PSF, PSFconj = PSFconj) +#Running SLIT +sourcesl, Imsl = SLIT.SLIT(Image, Fkappa, kmax, niter, size, PSF, PSFconj, levels = levels) +#Stop clock elapsed = (time.clock()-start) - print('execution time:', elapsed, 'seconds') +#Reconstruction goodness real_source = newsource - source_error = np.sum(np.abs(real_source[np.where(real_source!=0)] -sourcesl[np.where(real_source!=0)])**2 /real_source[np.where(real_source!=0)]**2)/(np.size( np.where(real_source!=0))/2.) - - image_chi2 = np.std(Image-Imsl)**2/sigma**2 print('Residuals in source space', source_error) print('Residuals in image space',image_chi2) - +#Display of results for i in [1]: plt.figure(2) # plt.suptitle('FISTA: error per pixel on the source: '+str(source_error)+' image chi2:'+str(image_chi2)) diff --git a/Examples/Test_SLIT_HR.py b/Examples/Test_SLIT_HR.py new file mode 100644 index 0000000..4ee2e7b --- /dev/null +++ b/Examples/Test_SLIT_HR.py @@ -0,0 +1,152 @@ +import pyfits as pf +import matplotlib.pyplot as plt +import numpy as np +import matplotlib.cm as cm +from scipy import signal as scp +import SLIT +import time +from scipy import signal as scp +import warnings +warnings.simplefilter("ignore") + +#Example of a run of the SLIT algorithm on simulated images. +#Here the first part of the file shows how simulations are generated. +#For users intereseted only in seeing the code run, have a look at the running SLIT section. +#The command line that builds the Fkappa operator is also of outmost importance. + + +###############################Simulation############################### +def SIE(x0,y0,n1,n2,b,beta,q,xc,theta): + kappa = np.zeros((n1,n2)) + x,y = np.where(kappa == 0) + eps = (1-q**2)/(1+q**2) + up = b**(beta-1) + pre = up/(2*(1-eps)**((beta-1)/2)) + count = 0 + theta = theta*np.pi/180. + for i in x: + Xr = (x[count]-x0)*np.cos(theta)-(y[count]-y0)*np.sin(theta) + Yr = (x[count]-x0)*np.sin(theta)+(y[count]-y0)*np.cos(theta) + kappa[x[count],y[count]] = pre/((xc**2.)/(1.-eps)+(Xr)**2.+((Yr)**2.)/q**2.)**((beta-1.)/2.) + + count += 1 + return kappa + +#Source light profile +newsource = pf.open('../Files/Source_HR.fits')[0].data +##N1,N2 are the numbers of pixels in the image plane. +nt1= 100 +nt2 = 100 +#Size ratio of the source to image number of pixels +size = 2 + +#PSF +PSF0 = pf.open('../Files/PSF.fits')[0].data +PSF = PSF0[1:,1:] +PSFconj = np.real(np.fft.ifft2(np.conjugate(np.fft.fft2(PSF0[:-1,:-1])))) +PSFconj=PSFconj/np.sum(PSFconj) +PSF = PSF/np.sum(PSF) + +## Lens mass distribution. +b = 1.53/0.05 +xc = 0.95 +q = 0.71 +betata = 2.1 +thetata = 25.2 +kappa = SIE(nt1/2.+50,nt2/2.+50,nt1+100,nt2+100,b,betata,q,xc, thetata) +#Mapping between lens and source IMPORTANT +Fkappa = SLIT.Lens.F(kappa, nt1,nt2, size,nt1/2.,nt2/2.) + + +#Generation of lensed source +I2 = SLIT.Lens.source_to_image(newsource, nt1 ,nt2 , Fkappa) + +#Noise levels +SNR = 500 +sigma = np.sqrt(np.sum(I2**2)/SNR/(nt1*nt2*size**2)) + +#Convolution by the PSF and generation of the final image +I2 = scp.fftconvolve(I2, PSF, mode = 'same') + +#Final simulated image +Image = I2+np.random.randn(nt1,nt2)*sigma + + +################################Running SLIT############################ +#Parameters +kmax = 5 +niter =50 +levels = [0] + +#Comment the following to have the level estimation routine run (takes more time) +levels = pf.open('../Files/Noise_levels_SLIT_HR.fits')[0].data + +#Start clock +start = time.clock() + +#Running SLIT +sourcesl, Imsl = SLIT.SLIT(Image, Fkappa, kmax, niter, size, PSF, PSFconj, levels = levels) + +#Stop clock +elapsed = (time.clock()-start) +print('execution time:', elapsed, 'seconds') + +#Reconstruction goodness +real_source = newsource +source_error = np.sum(np.abs(real_source[np.where(real_source!=0)] + -sourcesl[np.where(real_source!=0)])**2 + /real_source[np.where(real_source!=0)]**2)/(np.size( + np.where(real_source!=0))/2.) + +image_chi2 = np.std(Image-Imsl)**2/sigma**2 +print('Residuals in source space', source_error) +print('Residuals in image space',image_chi2) + +#Display of results +for i in [1]: + plt.figure(2) + # plt.suptitle('FISTA: error per pixel on the source: '+str(source_error)+' image chi2:'+str(image_chi2)) + # plt.subplot(2,3,1) + plt.title('Source from SLIT') + plt.imshow((sourcesl), vmin = np.min(real_source), vmax = np.max(real_source),cmap = cm.gist_stern, interpolation = 'nearest') + plt.axis('off') + plt.colorbar() +# plt.subplot(2,3,2) + plt.figure(3) + plt.title('Original source') + plt.imshow(real_source, cmap = cm.gist_stern, interpolation = 'nearest') + plt.axis('off') + plt.colorbar() + # plt.subplot(2,3,3) + plt.figure(4) + plt.title('Lensed source') + plt.imshow(Image, cmap = cm.gist_stern, interpolation = 'nearest') + plt.axis('off') + plt.colorbar() + plt.figure(41) + plt.title('Reconstructed lensed source') + plt.imshow(Imsl, vmin = np.min(Image), vmax = np.max(Image), cmap = cm.gist_stern, interpolation = 'nearest') + plt.axis('off') + plt.colorbar() + # plt.subplot(2,3,4) + plt.figure(5) + plt.title('relative difference') + diff = (real_source-sourcesl)/real_source + diff[np.where(real_source==0)] = 0 + diff[np.where(diff>1)]= np.log(0.) + plt.imshow(np.abs(diff), vmax = 1., vmin = 0., cmap = cm.gist_stern, interpolation = 'nearest' ) + plt.axis('off') + plt.colorbar() + # plt.subplot(2,3,5) + plt.figure(6) + plt.title('difference with reconstructed lensed image') + plt.imshow(Image-Imsl, vmin = -5*sigma, vmax = 5*sigma, cmap = cm.gist_stern, interpolation = 'nearest') + plt.axis('off') + plt.colorbar() + # plt.subplot(2,3,6) + plt.figure(7) + plt.title('difference with true source') + plt.imshow((np.abs(real_source-sourcesl)), cmap = cm.gist_stern, interpolation = 'nearest') + plt.axis('off') + plt.colorbar() +plt.show() diff --git a/Examples/Test_SLIT_MCA.py b/Examples/Test_SLIT_MCA.py index 7ed074a..879fa38 100644 --- a/Examples/Test_SLIT_MCA.py +++ b/Examples/Test_SLIT_MCA.py @@ -1,18 +1,22 @@ - -import Lens_better as Lens +from SLIT import Lens import pyfits as pf import matplotlib.pyplot as plt import numpy as np import matplotlib.cm as cm from scipy import signal as scp -import gaussian as gs import SLIT as slit import time from scipy import signal as scp import warnings warnings.simplefilter("ignore") +#Example of a run of the SLIT_MCA algorithm on simulated images. +#Here the first part of the file shows how simulations are generated. +#For users intereseted only in seeing the code run, have a look at the running SLIT_MCA section. +#The command line that builds the Fkappa operator is also of outmost importance. + +###############################Simulation############################### def SIE(x0,y0,n1,n2,b,beta,q,xc,theta): kappa = np.zeros((n1,n2)) x,y = np.where(kappa == 0) @@ -29,141 +33,73 @@ def SIE(x0,y0,n1,n2,b,beta,q,xc,theta): count += 1 return kappa - -real_source = pf.open('/Users/remy/Desktop/These/MuSCADeT/HFF/Abell2744/Colour_images_recombined.fits')[0].data#pf.open('imp_lens.fits')[0].data +#Source light profile +newsource = pf.open('../Files/source.fits')[0].data ##N1,N2 are the numbers of pixels in the image plane. nt1= 100 nt2 = 100 +#Size ratio of the source to image number of pixels +size = 1 -na1 = nt1 -na2 = nt2 +#PSF +PSF0 = pf.open('../Files/PSF.fits')[0].data +PSF = PSF0[1:,1:] +PSFconj = np.real(np.fft.ifft2(np.conjugate(np.fft.fft2(PSF0[:-1,:-1])))) +PSFconj=PSFconj/np.sum(PSFconj) +PSF = PSF/np.sum(PSF) + ## Lens mass distribution. b = 1.53/0.05 xc = 0.95 q = 0.71 betata = 2.1 thetata = 25.2 +kappa = SIE(nt1/2.+50,nt2/2.+50,nt1+100,nt2+100,b,betata,q,xc, thetata) +#Mapping between lens and source IMPORTANT +Fkappa = slit.Lens.F(kappa, nt1,nt2, size,nt1/2.,nt2/2.) -gal0 = real_source[0,2871-nt1/2:2871+nt1/2,1878-nt2/2:1878+nt2/2]#gs.sersic(nt1,nt2,nt1/2.,nt2/2.,2,0.08,0.1,0.2,1) - -gal2 = gs.gaussian(nt1,nt2,nt1/2.+3,3-nt2/2.,4,0.05,0.12,0) - - - -#gal0[-3:,-3:] = 0 -#gal2[-3:,-3:]=0 -extra = 100 -kappa = SIE(na1/2.+extra/2.,na2/2.+extra/2.,na1+extra,na2+extra,b,betata,q,xc, thetata) -#kappa2 = SIE(na1/2.+extra/2.,na2/2.+extra/2.,na1+extra,na2+extra,b+0.1,betata,q,xc, thetata) -#kappa = barkana(na1/2.,na2/2.,na1,na2,b,betata,q,thetata) - -hdus = pf.PrimaryHDU(kappa) -lists = pf.HDUList([hdus]) -lists.writeto('kappa.fits', clobber=True) - -x00 =3743+1#4843-35 -y00 =2030-2#6420+15 -size = 1 -newsource =real_source[0,x00-nt1*size/2.:x00+nt1*size/2.,y00-nt2*size/2.:y00+nt2*size/2.] -s = wine.MCA.MAD(newsource) -newsource,c = wine.MCA.mr_filter(newsource, 20,5, s, lvl = 5) -gal0,c = wine.MCA.mr_filter(gal0, 20,7, s, lvl = 5) -gal0 = (gal0-np.min(gal0))*4 -newsource = newsource-np.min(newsource) - - - -plt.imshow(gal0, cmap = cm.gist_stern, interpolation = 'nearest') -plt.axis('off') -plt.title('Lens galaxy') -plt.colorbar(); plt.show() - -start = time.clock() -theta = Lens.F(kappa, nt1,nt2, size,na1/2.,na2/2.) -#theta2 = Lens.F(kappa2, nt1,nt2, size,na1/2.,na2/2.) -elapsed = (time.clock()-start) - - -npsf1 = 65. -npsf2 = 65. -PSF0 = pf.open('PSF.fits')[0].data#('HST_1430/psf1.fits')[0].data## -PSF = PSF0[1:,1:] -PSFconj = np.real(np.fft.ifft2(np.conjugate(np.fft.fft2(PSF0[:-1,:-1])))) -PSFconj=PSFconj/np.sum(PSFconj) -PSF = PSF/np.sum(PSF) - - -##PSF0 = gs.gaussian(npsf1,npsf2,33.,33.,1.,0.95,0.95,0) -##PSF = gs.gaussian(npsf1,npsf2,32.,32.,1.,0.95,0.95,0) -##PSF[npsf2-1,npsf1-2] = 0 -##PSF0[npsf2-1,npsf1-2] = 0 +#Lens galaxy light profile +gal0 = pf.open('../Files/Galaxy.fits')[0].data +#Generation of lensed source +I2 = slit.Lens.source_to_image(newsource, nt1 ,nt2 , Fkappa) -xt,yt = np.where(np.zeros((nt1,nt2))==0) -I2 = Lens.source_to_image(newsource, nt1 ,nt2 , theta) -PSF0 = np.abs(PSF0/np.sum(PSF0)) -PSF = np.abs(PSF/np.sum(PSF)) - -##PSFconj = np.real(np.fft.ifft2(np.conjugate(np.fft.fft2(PSF)))) - - -I3 = np.copy(I2) -I2 = scp.fftconvolve(I2, PSF, mode = 'same') +#Noise levels SNR = 500 sigma = np.sqrt(np.sum(I2**2)/SNR/(nt1*nt2*size**2)) -gal = np.copy(gal0) -gal = scp.fftconvolve(gal, PSF, mode = 'same') - -Image = I2+gal+np.random.randn(nt1,nt2)*sigma - -plt.imshow(Image, cmap = cm.gist_stern, interpolation = 'nearest') -plt.title('Lens system') -plt.axis('off') -plt.colorbar() -plt.show() - -plt.imshow(newsource, cmap = cm.gist_stern, interpolation = 'nearest') -plt.title('Background source') -plt.axis('off') -plt.colorbar() -plt.show() +#Convolution by the PSF and generation of the final image +I2 = scp.fftconvolve(I2, PSF, mode = 'same') -plt.imshow(kappa, cmap = cm.gist_stern, interpolation = 'nearest') -plt.title('surface mass density') -plt.axis('off') -plt.colorbar() -plt.show() +#Convolution of the observed image +simu = scp.fftconvolve(gal0.astype(float)+I2, PSF.astype(float), mode = 'same') +#Sotring the convolved lens light profile: +gal = scp.fftconvolve(gal0.astype(float), PSF.astype(float), mode = 'same') +#Final simulated image +Image = simu+np.random.randn(nt1,nt2)*sigma +################################Running SLIT############################ +#Parameters kmax = 5 -niter =500 -kinter = 6 -ninter = 1 -pos_real = 1 -pos_wave = 0 -weight = 1 -#PSF=[0,0] -riter =200 - -FB = 0 -soft = 1 -repeat = 1 -PD = 0 +niter =100 +riter =50 +levels = [0] +#Comment the following to have the level estimation routine run (takes more time) +levels = pf.open('../Files/Noise_levels_SLIT_MCA.fits')[0].data + +#Start clock start = time.clock() -sourcesl, Imsl, G = slit.SLIT_old(Image, theta, kmax, niter,riter, kinter, ninter, size, - pos = pos_real, posit = pos_wave, decrease = 1, soft =soft, FB = FB, - reweighting = weight, Ginit =0, PSF = PSF, PSFconj = PSFconj, repeat = repeat, mrfilter = PD) +#Running SLIT_MCA +sourcesl, Imsl, G = slit.SLIT_MCA(Image, Fkappa, kmax, niter,riter, size,PSF, PSFconj, levels = levels) +#Stop clock elapsed = (time.clock()-start) - print('execution time:', elapsed, 'seconds') -source = sourcesl -Im=Imsl real_source = newsource @@ -172,105 +108,83 @@ def SIE(x0,y0,n1,n2,b,beta,q,xc,theta): /real_source[np.where(real_source!=0)]**2)/(np.size( np.where(real_source!=0))/2.) -GSF = np.copy(G) - - image_chi2 = np.std(Image-Imsl-G)**2/sigma**2 print('Residuals in source space', source_error) print('Residuals in image space',image_chi2) -hdus = pf.PrimaryHDU(sourcesl) -lists = pf.HDUList([hdus]) -lists.writeto('Souce.fits', clobber=True) - -hdus = pf.PrimaryHDU(G) -lists = pf.HDUList([hdus]) -lists.writeto('Galaxy_SLIT.fits', clobber=True) - - -hdus = pf.PrimaryHDU(real_source-source) -lists = pf.HDUList([hdus]) -lists.writeto('Diff.fits', clobber=True) - for i in [1]: ###Source - # plt.suptitle('FISTA: error per pixel on the source: '+str(source_error)+' image chi2:'+str(image_chi2)) -# plt.subplot(4,3,2) + plt.figure(0) plt.title('Source from SLIT') plt.imshow((sourcesl), vmin = np.min(real_source), vmax = np.max(real_source), cmap = cm.gist_stern, interpolation = 'nearest') plt.axis('off') plt.colorbar() - # plt.subplot(4,3,1) + plt.figure(1) plt.title('Original image of the source') plt.imshow(real_source, cmap = cm.gist_stern, interpolation = 'nearest') plt.axis('off') plt.colorbar() - # plt.subplot(4,3,3) + plt.figure(2) plt.title('relative difference') - diff = (real_source-sourcesl)#/real_source - # diff[np.where(real_source==0)] = 0 -# diff[np.where(diff>1)]= np.log(0.) + diff = (real_source-sourcesl) plt.imshow((np.abs(diff)), cmap = cm.gist_stern, interpolation = 'nearest') plt.axis('off') plt.colorbar() ####Lensed source - # plt.subplot(4,3,4) plt.figure(3) plt.title('Original lensed galaxy') plt.imshow(I2, cmap = cm.gist_stern, interpolation = 'nearest') plt.axis('off') plt.colorbar() - # plt.subplot(4,3,5) + plt.figure(4) plt.title('reconstructed lensed source') plt.imshow((Imsl), vmin = np.min(I2), vmax = np.max(I2), cmap = cm.gist_stern, interpolation = 'nearest') plt.axis('off') plt.colorbar() - # plt.subplot(4,3,6) + plt.figure(5) plt.title('error on the source in image plane') plt.imshow((I2-Imsl), cmap = cm.gist_stern, interpolation = 'nearest') plt.axis('off') plt.colorbar() ###Galaxy - # plt.subplot(4,3,7) plt.figure(6) plt.title('Original galaxy') plt.imshow((gal0), cmap = cm.gist_stern, interpolation = 'nearest') plt.axis('off') plt.colorbar() - # plt.subplot(4,3,8) + plt.figure(12) plt.title('Estimated Galaxy') plt.imshow((G), vmin = np.min(gal0), vmax = np.max(gal0), cmap = cm.gist_stern, interpolation = 'nearest') plt.axis('off') plt.colorbar() - # plt.subplot(4,3,9) + plt.figure(7) plt.title('Error on the galaxy') plt.imshow((gal-G), cmap = cm.gist_stern, interpolation = 'nearest') plt.axis('off') plt.colorbar() ###Image - # plt.subplot(4,3,10) plt.figure(8) plt.title('Image') plt.imshow(Image, cmap = cm.gist_stern, interpolation = 'nearest') plt.axis('off') plt.colorbar() - # plt.subplot(4,3,11) + plt.figure(9) plt.title('Reconstructed image') - plt.imshow(Imsl+GSF, cmap = cm.gist_stern, interpolation = 'nearest') + plt.imshow(Imsl+G, cmap = cm.gist_stern, interpolation = 'nearest') plt.axis('off') plt.colorbar() - # plt.subplot(4,3,12) + plt.figure(10) plt.title('difference with reconstructed image') - plt.imshow(Image-Imsl-GSF,cmap = cm.gist_stern, interpolation = 'nearest', vmin = -5*sigma, vmax = 5*sigma)#slit.fft_convolve(Im,PSF) + plt.imshow(Image-Imsl-G,cmap = cm.gist_stern, interpolation = 'nearest', vmin = -5*sigma, vmax = 5*sigma)#slit.fft_convolve(Im,PSF) plt.axis('off') plt.colorbar() diff --git a/Examples/Test_sparsity.py b/Examples/Test_sparsity.py new file mode 100644 index 0000000..7284ff6 --- /dev/null +++ b/Examples/Test_sparsity.py @@ -0,0 +1,104 @@ +from SLIT import Lens +import pyfits as pf +import matplotlib.pyplot as plt +import numpy as np +import matplotlib.cm as cm +from scipy import signal as scp +from SLIT import wave_transform as mw +import time +from scipy import signal as scp +import SLIT as slit +from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes +from mpl_toolkits.axes_grid1.inset_locator import mark_inset +import warnings +warnings.simplefilter("ignore") + + + +S = pf.open('../Files/source.fits')[0].data +G = pf.open('../Files/Galaxy.fits')[0].data +##Sizes in image and source planes +nt1= 100 +nt2 = 100 +size = 1 + +#Mass profile of the lens +kappa = pf.open('../Files/kappa.fits')[0].data +Fkappa = Lens.F(kappa, nt1,nt2, size,nt1/2.,nt2/2.) +lensed = slit.lens_one(Fkappa, nt1,nt2, size) + +#Starlet transforms of the lens and source in their respective planes +wG = mw.wave_transform(G, lvl = 6) +wS = mw.wave_transform(S, lvl = 6) +#Lensed source +FS = Lens.source_to_image(S, nt1, nt2,Fkappa) +#Unlensed lens +FG = Lens.image_to_source(G, size, Fkappa, lensed=lensed) +#Starlet transform of the unlensed lens +wFG = mw.wave_transform(FG, 6) +#Starlet transform of the lensed +wFS = mw.wave_transform(FS, 6) + +def mk_sort(X): + Y = np.sort(np.resize(np.abs(X), X.size)) + return Y[::-1] + +#Function that computes the reconstruction error from the p% highest coefficients +def error_rec_from(X, p, wave = 0): + Xcopy = np.copy(X) + Y = mk_sort(X) + + ymin = Y[p*Y.size/1000.] + Xcopy[np.abs(X)0: - PSF_norm = spectralNorm(nt1,nt2,20,0.00000001, PSF_deproj, PSF_proj) - mu = 1./(Star_norm_s*PSF_norm)**2 - else: - mu = 1./(Star_norm_s*F_norm)**2 + PSF_norm = spectralNorm(nb1,nb2,20,0.00000001, PSF_deproj, PSF_proj) + mu = 1./(Star_norm_s*PSF_norm*F_norm)**2 - - Res2 = [] + #Initialisation R=0 - - print(mu) S=np.random.randn(nb1,nb2)*sigma0 + FS = 0 alpha = np.zeros((lvl,nb1,nb2)) csi = np.copy(alpha) t =1. - - i=0 Res1 = [] - S=np.random.randn(nb1,nb2)*sigma0 - alpha = np.zeros((lvl,nb1,nb2)) - csi = np.copy(alpha) + i=0 while i < niter: - if np.sum(PSF) != 0: - R = mu*scp.fftconvolve((img-Im),(PSFconj),mode = 'same') - Rlens = Lens.image_to_source(R, size, beta, lensed=lensed) - else: - R = mu*(img-Im) - Rlens = Lens.image_to_source(R, size, beta, lensed=lensed) - + print(i) + #Gradient step + R = mu*scp.fftconvolve((img-FS),(PSFconj),mode = 'same') + Rlens = Lens.image_to_source(R, size, Fkappa, lensed=lensed) + + #Thresholding in starlet space alpha_new = csi+mw.wave_transform(Rlens, lvl, newwave = 1) - alpha_new= HT(alpha_new, kmax, levels, sigma0, soft = soft)*supp - - t_new = (1.+np.sqrt(1.+4.*t**2))/2. + alpha_new= ST(alpha_new, kmax, levels, sigma0)*supp + #Inertial step + t_new = (1.+np.sqrt(1.+4.*t**2))/2. csi = alpha_new+((t-1)/t_new)*(alpha_new-alpha) alpha = alpha_new + #Reconstruction of the source for next step S = mw.iuwt(csi) - Im = Lens.source_to_image(S, nt1, nt2,theta) - if np.sum(PSF)!=0: - Im = scp.fftconvolve(Im,PSF, mode = 'same') + FS = Lens.source_to_image(S, nt1, nt2,Fkappa) + FS = scp.fftconvolve(FS,PSF, mode = 'same')*mask t = np.copy(t_new) - Res1.append((np.std(img-Im)**2)/sigma0**2) + #Convergence condition + Res1.append((np.std(img-FS)**2)/sigma0**2) if i >20: if np.abs(Res1[i]-Res1[i-10])<0.001 and Res1[i]2: - lvlso = (scp.fftconvolve(noise[i,:,:]**2, wave_dirac[i,:,:]**2, - mode = 'same')) - else: - X = np.pad(noise,((np.int_((nn-nb1)/2),np.int_((nn-nb1)/2)),(np.int_((nn-nb2)/2),np.int_((nn-nb2)/2))), mode = 'edge') - Y = np.pad(wave_dirac[i],((np.int_((nn-nb1)/2),np.int_((nn-nb1)/2)),(np.int_((nn-nb2)/2),np.int_((nn-nb2)/2))), mode = 'edge') - lvlso = scp.fftconvolve(X**2,Y**2, - mode = 'same') - - lvlso = lvlso[(nn-nb1)/2.:(nn+nb1)/2.,(nn-nb2)/2.:(nn+nb2)/2.] - # plt.imshow((lvls)); plt.colorbar();plt.show() - - levels[i,:,:] = np.sqrt(np.abs(lvlso)) - return levels def spectralNorm(nx,ny,Niter,tol,f,finv): + ##DESCRIPTION: + ## Function that estimates the source light profile from an image of a lensed source given the mass density profile. + ## + ##INPUTS: + ## -nx,ny: shape of the input + ## -nz: number of decomposition scales (if the operator tis a multiscale decomposition for instance) + ## -Niter: number of iterations + ## -tol: tolerance error as a stopping criteria + ## -f: operator + ## -finv: inverse operator + ## + ##OUTPUTS: + ## -SspNorm: The spectral norm of the operator -## Inputs: -## nx: nombre de lignes l entree -## ny: nombre de colonnes de l entree -## nz: nombre d echelles (pour transformation en ondelette redondante) -## Niter: nombre d iterations -## tol: erreur de tolerance pour s arreter -## f: l operateur -## finv: l operateur inverse #### Initilize array with random numbers ### matA = np.random.randn(nx,ny) ### Normalize the input ### @@ -212,105 +400,60 @@ def spectralNorm(nx,ny,Niter,tol,f,finv): spNorm_new = LA.norm(matA) matA /= spNorm_new err = abs(spNorm_new - spNorm)/spNorm_new - print "Iteration:"+str(it)+",Error:"+str(err) spNorm = spNorm_new it += 1 return spNorm -def level_source(n1,n2,x,y,N1,N2,Fout,lvl,lensed, PSF=[0,0]): - - dirac = np.zeros((n1,n2)) - i1 = (np.max(x)+np.min(x))/2 - i2 = (np.max(y)+np.min(y))/2 - dirac[i1,i2] = 1. - wave_dirac = mw.wave_transform(dirac, lvl, newwave = 0) - n, w1, w2 = np.shape(wave_dirac) - if np.sum(PSF)!=0: - for j in np.linspace(0,n-1,n): - wave_dirac[j,:,:] = scp.fftconvolve(wave_dirac[j,:,:],PSF,mode = 'same') - - lensed_wave = np.zeros((n,N1,N2)) - for i in range(n): - lensed_wave[i,:,:] = Lens.image_to_source(wave_dirac[i,:,:],x,y,N1,N2,Fout,lensed=lensed)**2 - wave_sum = np.sum(np.sum(lensed_wave,1),1) - sig = np.sqrt(wave_sum) - return sig - -def lens_one(beta, nt1,nt2,size, ones = 0): - # n1,n2: size of the postage stamp in image plane - # N1,N2: size of the postage stamp in source plane +def lens_one(Fkappa, nt1,nt2,size): + ##DESCRIPTION: + ## Function that maps an all at one image to source plane. + ## + ##INPUTS: + ## -Fkappa: the mapping between source and image planes + ## -nt1,nt2: the shape of the image. + ## -size: the factor that scales the shape of the source relative to the shape of the image + ## + ##OUTPUTS: + ## -lensed: the projection to source plane of an all at aone image. dirac = np.ones((nt1,nt2)) - lensed = Lens.image_to_source(dirac, size,beta,lensed = [0])#multi = 1)# - if ones ==1: - lensed[np.where(lensed == 0)] =1 + lensed = Lens.image_to_source(dirac, size,Fkappa,lensed = [0]) return lensed -def unlens_one(F, xt,yt,nt1,nt2,nb1,nb2, ones = 0): - # n1,n2: size of the postage stamp in image plane - # N1,N2: size of the postage stamp in source plane - dirac = np.ones((nb1,nb2)) - unlensed = Lens.source_to_image(dirac,nt1, nt1,F) - if ones ==1: - unlensed[np.where(lensed == 0)] =1 - return unlensed - -def MOM(RG,RS,levelg,levels, lvlg,lvls, sigma): - nt1,nt2 = np.shape(RG) - ns1,ns2 = np.shape(RS) - levels[levels==0]=1 - wG = mw.wave_transform(RG,lvlg, newwave = 1)/levelg/sigma - wS = mw.wave_transform(RS,lvls, newwave = 1)/levels/sigma - - wGmax = np.max(wG) - wSmax = np.max(wS) - wmax = [wGmax, wSmax] - - k = np.min(wmax)+(max(wmax)-min(wmax))/100 - return k - - -def level(nt1,nt2,lvl = 6): - dirac = np.zeros((nt1,nt2)) - dirac[nt1/2,nt2/2] = 1 - wave_dirac = mw.wave_transform(dirac,lvl, newwave = 0) - - wave_sum = np.sqrt(np.sum(np.sum(wave_dirac**2,1),1)) - levels = np.multiply(np.ones((lvl,nt1,nt2)).T,wave_sum).T - - return levels - -def linorm(A,nit): - ns,nb = np.shape(A) - x0 = np.random.rand(nb) - x0 = x0/np.sqrt(np.sum(x0**2)) - - for i in np.linspace(0,nit-1,nit): - x = np.dot(A,x0) - xn = np.sqrt(np.sum(x**2)) - xp = x/xn - y = np.dot(A.T,xp) - yn = np.sqrt(np.sum(y**2)) - if yn < np.dot(y.T,x0) : - break - x0 = y/yn - - return xn - -def MAD(x,n=3,fil=1): - if fil == 1: - meda = med.median_filter(x,size = (n,n)) - else: - meda = np.median(x) - medfil = np.abs(x-meda) - sh = np.shape(x) - sigma = 1.48*np.median((medfil)) - return sigma - - -def HT(alpha, k, levels, sigma, pen = 0, soft = 0): +def MAD(x,n=3): + ##DESCRIPTION: + ## Estimates the noise standard deviation from Median Absolute Deviation + ## + ##INPUTS: + ## -x: a 2D image for which we look for the noise levels. + ## + ##OPTIONS: + ## -n: size of the median filter. Default is 3. + ## + ##OUTPUTS: + ## -S: the source light profile. + ## -FS: the lensed version of the estimated source light profile + meda = med.median_filter(x,size = (n,n)) + medfil = np.abs(x-meda) + sh = np.shape(x) + sigma = 1.48*np.median((medfil)) + return sigma + + +def ST(alpha, k, levels, sigma): + ##DESCRIPTION: + ## Soft thresholding operator. + ## + ##INPUTS: + ## -alpha: the starlet decomposition to be thresholded. + ## -k: the threshold in units of noise levels (usually 5). + ## -levels: the noise levels at each scale and location of the starlet decomposition. + ## -sigma: the noise standard deviation. + ## + ##OUTPUTS: + ## -alpha: The thresholded coefficients. lvl, n1,n2 = np.shape(alpha) th = np.ones((lvl,n1,n2))*k th[0,:,:] = th[0,:,:]+1 @@ -320,19 +463,30 @@ def HT(alpha, k, levels, sigma, pen = 0, soft = 0): alpha0 = np.copy(alpha) th = th*levels*sigma - if soft == 1: - alpha= np.sign(alpha0)*(np.abs(alpha0)-th) - alpha[np.where(np.abs(alpha)-th<0)]=0 - else: - alpha[np.where(np.abs(alpha)-th<0)] = 0 + alpha= np.sign(alpha0)*(np.abs(alpha0)-th) + alpha[np.where(np.abs(alpha)-th<0)]=0 - return alpha + return alpha -## -############################ SLIT MCA FISTA -def simulate_noise(nt1,nt2, size, beta, lensed, PSFconj): +def simulate_noise(nt1,nt2, size, Fkappa, lensed, PSFconj): + ##DESCRIPTION: + ## Simulates noise levels in source plane from lensing operator and convolution operator. + ## + ##INPUTS: + ## -nt1,nt2: the shape of the images for which to simulate noise maps. + ## -size: scaling factor for the shape of the source. + ## -Fkappa: Projection operator between lens and source plane. + ## -lensed: mapping of an all at one image to source plane. + ## -PSFconj: the conjugate of the PSF + ## + ##OPTIONS: + ## -n: size of the median filter. Default is 3. + ## + ##OUTPUTS: + ## -S: the source light profile. + ## -FS: the lensed version of the estimated source light profile n = 500. nb1,nb2 = nt1*size, nt2*size lvl = np.int(np.log2(nb1)) @@ -349,188 +503,12 @@ def simulate_noise(nt1,nt2, size, beta, lensed, PSFconj): for i in range(np.int(n)): noise = np.random.randn(nt1,nt2) noise = scp.fftconvolve(noise,PSFconj,mode = 'same') - noise_lens = (Lens.image_to_source(noise, size, beta, lensed =lensed)) + noise_lens = (Lens.image_to_source(noise, size, Fkappa, lensed =lensed)) noise_lens[noise_lens ==0] = 1 - storage[:,:,:,i] = mw.wave_transform(noise_lens, lvl = lvl) - w_levels = np.std(storage, axis = 3) - return w_levels -def SLIT_MCA(img, theta, kmax, niter, riter, kinter, ninter, size,decrease = 0, FB = 0, - pos = True, K=0, posit = False,reweighting = False, soft = 0, mask = [0,0], - PSF = [0], PSFconj = [0,0], gauss = 0, Ginit=0, repeat = 1, mrfilter = 0): - - beta = theta - #Initialisation of the lensed sourced - Im = np.copy(img)*0. - - nt1,nt2 = np.shape(img) - #Initialisation of the source - nb1= nt1*size - nb2 = nt2*size - - lvlg = np.int(np.log2(nt1)) - lvls = np.int(np.log2(nb1)) - #Masking if required - if np.sum(mask) == 0: - mask = np.ones((nt1,nt2)) - img = img*mask - M = [0] - #Noise in image plane - sigma0 = wine.MCA.MAD(img) - #Mapping of an all-at-one image - lensed = lens_one(beta, nt1,nt2, size) - bound = lensed/lensed - bound[lensed == 0] = 0 - #Noise level map with convolution by wavelet elements +weird normalisation - debound = mw.wave_transform(1-bound, lvls) - - levelg = level(nt1,nt1,lvl=lvlg) - - levels = simulate_noise(nt1,nt2, size, beta, lensed, PSFconj) - levels0 = lens_wnoise(nb1,nb2,size,beta, lvl = lvls) - - S =Lens.image_to_source(img, size, beta, lensed = lensed) - - #K thresholds in units of sigma - if np.sum(Ginit) == 0: - Ginit = np.random.randn(nt1,nt2)*sigma0 - alphag = mw.wave_transform(Ginit, lvlg, newwave = 1) - wFS = np.ones((lvlg,nt1,nt2)) - RG = np.copy(img) - RS = np.copy(S) - if np.sum(PSF) != 0: - RS = scp.fftconvolve(RS,(PSFconj),mode = 'same') - Res1 = [] - Res2 = [] - R=0 - i = 0 - S=np.random.randn(nb1,nb2)*sigma0 - G = Ginit - FS = 0 - k = 0 - alphas = np.zeros((lvls,nb1,nb2)) - csis = np.copy(alphas) - csig = np.copy(alphag) - K_s = np.zeros((niter)) - norm = np.zeros((3,niter)) - tg =1 - -## -##Compute spectral norms - def Finv_apply(I): - return Lens.image_to_source(I, size, beta, lensed = lensed) - def F_apply(Si): - return Lens.source_to_image(Si, nt1, nt2,theta) - def PSF_apply(i): - return scp.fftconvolve(i,PSF,mode = 'same') - def PSFT_apply(ii): - return scp.fftconvolve(ii,PSFconj,mode = 'same') - def star(x): - return mw.wave_transform(x, lvlg, newwave = 1) - def star_inv(x): - return mw.iuwt(x) - - F_norm = spectralNorm(nt1,nt2,20,1e-10,F_apply,Finv_apply) - Star_norm_im = spectralNorm(nt1,nt2,20,1e-10,star,star_inv) - Star_norm_s = spectralNorm(nb1,nb2,20,1e-10,star,star_inv) - if np.sum(PSF)>0: - PSF_norm = spectralNorm(nt1,nt2,20,0.0000001,PSF_apply,PSFT_apply) - muG = 1./(Star_norm_im**2) - muS = 1./(Star_norm_s*F_norm*PSF_norm)**2 - else: - muG = 1./(Star_norm_im**2) - muS = 1./(Star_norm_s*F_norm)**2 - - mu = 2./(Star_norm_s**2+F_norm**2) - - def HGT(X): - return X - def HG(X): - return X - def Phi(X): - return mw.iuwt(X) - def PhiT(X): - return mw.wave_transform(X, lvlg, newwave = 1) - - def HFT(X): - return Lens.image_to_source(X, size, beta, lensed=lensed) - def HF(X): - return Lens.source_to_image(X, nt1,nt2, theta) - - RS = Lens.image_to_source(img, size, beta, lensed = lensed) - RG = np.copy(img) - G = np.random.randn(nt1,nt2)*sigma0 - S = np.random.randn(nb1,nb2)*sigma0 - FS = np.random.randn(nt1,nt2)*sigma0 - i=0 - k = kmax - while i < niter: - print(i) - DS = img-G - j = 0 - ts = 1 - tr = np.zeros(riter) - - tauG = 1 - tauS = 1 - - while j < riter: - RS = muS*(DS-FS) - if np.sum(PSF)!=0: - RS = scp.fftconvolve(RS,(PSFconj),mode = 'same') - RS = Lens.image_to_source(RS, size, beta, lensed=lensed) - - alphas_new = csis+mw.wave_transform(RS, lvls, newwave = 1) - - alphas_new = HT(alphas_new, k, levels, sigma0, soft = soft)*bound - ts_new = (1.+np.sqrt(1.+4.*ts**2))/2. - csis = alphas_new+((ts-1)/ts_new)*(alphas_new-alphas) - alphas = alphas_new - ts = np.copy(ts_new) - S = mw.iuwt(csis) - - FS = Lens.source_to_image(S, nt1, nt2,theta) - - - if np.sum(PSF)!=0: - FS = scp.fftconvolve(FS,PSF,mode = 'same') - tr[j] = (np.std(img-FS-G)**2)/sigma0**2 - j+=1 - S = mw.iuwt(alphas) - FS = Lens.source_to_image(S, nt1, nt2,theta) - if np.sum(PSF)!=0: - FS = scp.fftconvolve(FS,PSF,mode = 'same') - - DG = img-FS - RG = muG*(DG-G) - - alphag_new = csig+mw.wave_transform(RG, lvlg, newwave = 1) - alphag_new = HT(alphag_new, kmax, levelg, sigma0, soft = soft) - tg_new = (1.+np.sqrt(1.+4.*tg**2))/2. - csig = alphag_new+((tg-1)/tg_new)*(alphag_new-alphag) - alphag = alphag_new - - tg = np.copy(tg_new) - G = mw.iuwt(csig) - - newres = (np.std(img-FS-G)**2)/sigma0**2 - K_s[i] = newres - - res = np.copy(newres) - - i = i+1 - - - S[np.where(S<0)] = 0 - FS = Lens.source_to_image(S, nt1, nt2,theta) - if np.sum(PSF)!=0: - FS = scp.fftconvolve(FS,PSF,mode = 'same') - - return S, FS,G - diff --git a/build/lib/SLIT/Lens.py b/build/lib/SLIT/Lens.py index 97ad3aa..74bcfae 100644 --- a/build/lib/SLIT/Lens.py +++ b/build/lib/SLIT/Lens.py @@ -5,61 +5,9 @@ import warnings warnings.simplefilter("ignore") -def SIS(n1, n2, x0,y0, Re): - # - #Generate a lens mass profile as a Singular Isothermic Spheroid with size n1xn2 - #Centered in (x0,y0) with Einstein Radii Re - x,y = np.where(np.zeros([n1,n1])==0) - kappa = np.zeros([n1,n2]) - kappa[x,y] = Re/(2*np.sqrt((x-x0)**2+(y-y0)**2)) +#Tool box for lensing - return kappa - -def NSIS(n1,n2,x0,y0,Re,Rc): - # - #Generate a lens mass profile as a NSingular Isothermic Spheroid with size n1xn2 - #Centered in (x0,y0) with Einstein Radii Re - - x,y = np.where(np.zeros([n1,n2])==0) - kappa = np.zeros([n1,n2]) - kappa[x,y] = Re/Rc*(1+((x-x0)**2+(y-y0)**2)/(2*Rc**2))*(1+((x-x0)**2+(y-y0)**2)/(Rc**2))**(-1.5) - return kappa - -def Dcomplex(x,y): - n1 = np.sqrt(x.size) - n2 = np.sqrt(y.size) - x = x-n1/2 - y = y-n2/2 - D1 = np.zeros((n1,n2)) - D2 = np.zeros((n1,n2))+0. - - for i in range(x.size): - rho = (x[i]**4+y[i]**4)+2*(x[i]*y[i])**2 - D1[x[i]+n1/2,y[i]+n2/2] = (y[i]**2-x[i]**2)/rho - D2[x[i]+n1/2,y[i]+n2/2] = (-2*x[i]*y[i])/rho - - return D1,D2 - -def gamma(kappa): - n1,n2 = np.shape(kappa) - g1 = np.zeros((n1,n2)) - g2 = np.zeros((n1,n2)) - x,y = np.where(g1 == 0) - D1,D2 = Dcomplex(x,y) - g1 = (1/np.pi)*scp.convolve2d(D1,kappa, mode = 'same', boundary = 'wrap') - g2 = (1/np.pi)*scp.convolve2d(D2,kappa, mode = 'same', boundary = 'wrap') - return g1,g2 - - -def make_big(array,l): - n1,n2 = np.shape(array) - for i in range(l): - array = np.append(np.reshape(array[0,:],(1,array[0,:].size)),array,axis = 0) - array = np.append(array,np.reshape(array[-1,:],(1,array[-1,:].size)),axis = 0) - array = np.append(np.reshape(array[:,0],(array[:,0].size,1)),array,axis = 1) - array = np.append(array,np.reshape(array[:,-1],(array[:,-1].size,1)),axis = 1) - return array def alpha_def(kappa, n1,n2,x0,y0,extra): @@ -71,8 +19,8 @@ def alpha_def(kappa, n1,n2,x0,y0,extra): #Coordonnees de la grille de l'espace image [x,y] = np.where(np.zeros([n1,n2])==0) - xc = np.reshape((x)-x0+1.,(n1,n2))#+extra/2. - yc = np.reshape((y)-y0+1.,(n1,n2))#+extra/2. + xc = np.reshape((x)-x0+1.,(n1,n2)) + yc = np.reshape((y)-y0+1.,(n1,n2)) r = np.reshape((xc**2+yc**2),(n1,n2)) @@ -83,9 +31,7 @@ def alpha_def(kappa, n1,n2,x0,y0,extra): taby[lx,ly]=0 l = 0 - # tabx = make_big(tabx,l) -# taby = make_big(taby,l) -# kappa = make_big(kappa,l) + kappa = kappa.astype(float) tabx = tabx.astype(float) @@ -114,18 +60,8 @@ def F(kappa, nt1,nt2, size, x0, y0, extra=100,local = False, alpha_file = 'none' if local == False: nk1,nk2 = np.shape(kappa) alpha = np.zeros((2,nt1,nt2)) - #for i in range(np.size(xk)): - # alpha[:,xk[i],yk[i]] = alpha_def(np.array([xk[i],yk[i]]),kappa) - #alpha_x = alpha[0,:,:] - #alpha_y = alpha[1,:,:] alpha_x,alpha_y = alpha_def(kappa,nt1,nt2,x0,y0,extra) - # alpha = alpha_def(np.array([xk,yk]),kappa)# -# alpha_x = alpha[0] -# alpha_y = alpha[1] - - # alpha_x = pf.open('Lenstools/dpl_y.fits')[0].data/(0.8) # - # alpha_y = pf.open('Lenstools/dpl_x.fits')[0].data/(0.8) alpha[0,:,:] = alpha_x alpha[1,:,:] = alpha_y @@ -151,11 +87,8 @@ def F(kappa, nt1,nt2, size, x0, y0, extra=100,local = False, alpha_file = 'none' #Scaling of the source grid #Scaling of the deflection grid - xa = xa*(np.float(nt1)/np.float(na1))-1#./2 - ya = ya*(np.float(nt2)/np.float(na2))-1#./2 - #Scaling of the deflection angles - # alpha_x = alpha_x -# alpha_y = alpha_y + xa = xa*(np.float(nt1)/np.float(na1))-1 + ya = ya*(np.float(nt2)/np.float(na2))-1 #Setting images coordinates in 2d xa2d = np.reshape(xa,(na1,na2)) ya2d = np.reshape(ya,(na1,na2)) @@ -180,12 +113,7 @@ def F(kappa, nt1,nt2, size, x0, y0, extra=100,local = False, alpha_file = 'none' F.append([0]) else: - F.append(np.int_(np.array(loc))) - hdus = pf.PrimaryHDU(alpha) - lists = pf.HDUList([hdus]) - lists.writeto('alphas.fits', clobber=True) - return F diff --git a/build/lib/SLIT/SLIT.py b/build/lib/SLIT/SLIT.py index 79bad78..315bd06 100644 --- a/build/lib/SLIT/SLIT.py +++ b/build/lib/SLIT/SLIT.py @@ -13,51 +13,71 @@ ##SLIT: Sparse Lens Inversion Technique -def SLIT(img, theta, kmax, niter, size, decrease = 0, mymy = 0, levels = [0], - pos = True, K=0, posit = False,weight = [0,0], soft = 0, mask = [0,0], PSF = [0], PSFconj = [0,0]): +def SLIT(img, Fkappa, kmax, niter, size, PSF, PSFconj, levels = [0], mask = [0]): + ##DESCRIPTION: + ## Function that estimates the source light profile from an image of a lensed source given the mass density profile. + ## + ##INPUTS: + ## -img: a 2-D image of a lensed source given as nt1xnt2 numpy array. + ## -Fkappa: an array giving the mapping between lens and source. This array is calculated from the lens mass density + ## using tools from SLIT.Lens + ## -kmax: the detection threshold in units of noise levels. We usualy set this value to 5 to get a 5 sigma + ## detection threshold. + ## -niter: maximal number of iterations of the algorithm. + ## -size: resoluution factor between lens and source grids such thathe size of the output source + ## will be nt1sizexnt2size + ## -PSF: the point spread function of the observation provided as a 2D array. + ## -PSFconj: The conjugate of the PSF. Usually computed via np.real(np.fft.ifft2(np.conjugate(np.fft.fft2(PSF0[:-1,:-1])))) + ## butthe user has to make sure that the conjugate is well centered. + ## + ##OPTIONS: + ## -levels: an array that contains the noise levels at each band of the wavelet decomposition of the source. + ## If not provided, the routine will compute the levels and save them in a fits file 'Noise_levels.fits' + ## so that they can be used at a later time. This option allows to save time when running the same + ## experiment several times. + ## -mask: an array of zeros and one with size nb1xnb2. The zeros will stand for masked data. + ## + ##OUTPUTS: + ## -S: the source light profile. + ## -FS: the lensed version of the estimated source light profile + ## + ##EXAMPLE: + ## S,FS = SLIT(img, Fkappa, 5, 100, 1, PSF, PSFconj) - beta = theta - #Initialisation of the lensed sourced - Im = np.copy(img)*0. nt1,nt2 = np.shape(img) - #Initialisation of the source - nb1= nt1*size - nb2 = nt2*size - + + #Size of the source + nb1,nb2 = nt1*size, nt2*size + #Number of starlet scales in source plane lvl = np.int(np.log2(nb2)) + #Masking if required if np.sum(mask) == 0: mask = np.ones((nt1,nt2)) img = img*mask - M = [0] + #Noise in image plane sigma0 = MAD(img) - #Mapping of an all-at-one image - lensed = lens_one(beta, nt1,nt2, size) - + #Mapping of an all-at-one image to source plane + lensed = lens_one(Fkappa, nt1,nt2, size) + #estimation of the frame of the image in source plane supp = np.zeros(lensed.shape) supp[lensed/lensed ==1] =1 - #Noise level map with convolution by wavelet elements +weird normalisation - + #Noise simulations to estimate noise levels in source plane if np.sum(levels)==0: - levels0 = lens_wnoise(nb1,nb2,lensed,lvl = lvl) - if np.sum(PSF)!=0: - levels = simulate_noise(nt1,nt2, size, beta, lensed, PSFconj) - else: - levels = levels0 - - levels0 = np.copy(levels) - S = Lens.image_to_source(img, size, beta, lensed = lensed) - k = kmax - + levels = simulate_noise(nt1,nt2, size, Fkappa, lensed, PSFconj) + #Saves levels + hdus = pf.PrimaryHDU(levels) + lists = pf.HDUList([hdus]) + lists.writeto('Noise_levels.fits', clobber=True) ##Compute spectral norms def Finv_apply(I): - return Lens.image_to_source(I, size, beta, lensed = lensed) + return Lens.image_to_source(I, size, Fkappa, lensed = lensed) def F_apply(Si): - return Lens.source_to_image(Si, nt1, nt2,theta) + return Lens.source_to_image(Si, nt1, nt2,Fkappa) def PSF_apply(i): return scp.fftconvolve(i,PSF,mode = 'same') def PSFT_apply(ii): @@ -70,73 +90,252 @@ def PSF_proj(x): return Finv_apply(PSFT_apply(x)) def PSF_deproj(x): return PSF_apply(F_apply(x)) - F_norm = spectralNorm(nb1,nb2,20,1e-10,F_apply,Finv_apply) Star_norm_s = spectralNorm(nb1,nb2,20,1e-10,star,star_inv) - if np.sum(PSF)>0: - PSF_norm = spectralNorm(nt1,nt2,20,0.00000001, PSF_deproj, PSF_proj) - mu = 1./(Star_norm_s*PSF_norm)**2 - else: - mu = 1./(Star_norm_s*F_norm)**2 + PSF_norm = spectralNorm(nb1,nb2,20,0.00000001, PSF_deproj, PSF_proj) + mu = 1./(Star_norm_s*PSF_norm*F_norm)**2 - - Res2 = [] + #Initialisation R=0 - - print(mu) S=np.random.randn(nb1,nb2)*sigma0 + FS = 0 alpha = np.zeros((lvl,nb1,nb2)) csi = np.copy(alpha) t =1. - - i=0 Res1 = [] - S=np.random.randn(nb1,nb2)*sigma0 - alpha = np.zeros((lvl,nb1,nb2)) - csi = np.copy(alpha) + i=0 while i < niter: - if np.sum(PSF) != 0: - R = mu*scp.fftconvolve((img-Im),(PSFconj),mode = 'same') - Rlens = Lens.image_to_source(R, size, beta, lensed=lensed) - else: - R = mu*(img-Im) - Rlens = Lens.image_to_source(R, size, beta, lensed=lensed) - + print(i) + #Gradient step + R = mu*scp.fftconvolve((img-FS),(PSFconj),mode = 'same') + Rlens = Lens.image_to_source(R, size, Fkappa, lensed=lensed) + + #Thresholding in starlet space alpha_new = csi+mw.wave_transform(Rlens, lvl, newwave = 1) - alpha_new= HT(alpha_new, kmax, levels, sigma0, soft = soft)*supp - - t_new = (1.+np.sqrt(1.+4.*t**2))/2. + alpha_new= ST(alpha_new, kmax, levels, sigma0)*supp + #Inertial step + t_new = (1.+np.sqrt(1.+4.*t**2))/2. csi = alpha_new+((t-1)/t_new)*(alpha_new-alpha) alpha = alpha_new + #Reconstruction of the source for next step S = mw.iuwt(csi) - Im = Lens.source_to_image(S, nt1, nt2,theta) - if np.sum(PSF)!=0: - Im = scp.fftconvolve(Im,PSF, mode = 'same') + FS = Lens.source_to_image(S, nt1, nt2,Fkappa) + FS = scp.fftconvolve(FS,PSF, mode = 'same')*mask t = np.copy(t_new) - Res1.append((np.std(img-Im)**2)/sigma0**2) + #Convergence condition + Res1.append((np.std(img-FS)**2)/sigma0**2) if i >20: if np.abs(Res1[i]-Res1[i-10])<0.001 and Res1[i]2: - lvlso = (scp.fftconvolve(noise[i,:,:]**2, wave_dirac[i,:,:]**2, - mode = 'same')) - else: - X = np.pad(noise,((np.int_((nn-nb1)/2),np.int_((nn-nb1)/2)),(np.int_((nn-nb2)/2),np.int_((nn-nb2)/2))), mode = 'edge') - Y = np.pad(wave_dirac[i],((np.int_((nn-nb1)/2),np.int_((nn-nb1)/2)),(np.int_((nn-nb2)/2),np.int_((nn-nb2)/2))), mode = 'edge') - lvlso = scp.fftconvolve(X**2,Y**2, - mode = 'same') - - lvlso = lvlso[(nn-nb1)/2.:(nn+nb1)/2.,(nn-nb2)/2.:(nn+nb2)/2.] - # plt.imshow((lvls)); plt.colorbar();plt.show() - - levels[i,:,:] = np.sqrt(np.abs(lvlso)) - return levels def spectralNorm(nx,ny,Niter,tol,f,finv): + ##DESCRIPTION: + ## Function that estimates the source light profile from an image of a lensed source given the mass density profile. + ## + ##INPUTS: + ## -nx,ny: shape of the input + ## -nz: number of decomposition scales (if the operator tis a multiscale decomposition for instance) + ## -Niter: number of iterations + ## -tol: tolerance error as a stopping criteria + ## -f: operator + ## -finv: inverse operator + ## + ##OUTPUTS: + ## -SspNorm: The spectral norm of the operator -## Inputs: -## nx: nombre de lignes l entree -## ny: nombre de colonnes de l entree -## nz: nombre d echelles (pour transformation en ondelette redondante) -## Niter: nombre d iterations -## tol: erreur de tolerance pour s arreter -## f: l operateur -## finv: l operateur inverse #### Initilize array with random numbers ### matA = np.random.randn(nx,ny) ### Normalize the input ### @@ -212,105 +400,60 @@ def spectralNorm(nx,ny,Niter,tol,f,finv): spNorm_new = LA.norm(matA) matA /= spNorm_new err = abs(spNorm_new - spNorm)/spNorm_new - print "Iteration:"+str(it)+",Error:"+str(err) spNorm = spNorm_new it += 1 return spNorm -def level_source(n1,n2,x,y,N1,N2,Fout,lvl,lensed, PSF=[0,0]): - - dirac = np.zeros((n1,n2)) - i1 = (np.max(x)+np.min(x))/2 - i2 = (np.max(y)+np.min(y))/2 - dirac[i1,i2] = 1. - wave_dirac = mw.wave_transform(dirac, lvl, newwave = 0) - n, w1, w2 = np.shape(wave_dirac) - if np.sum(PSF)!=0: - for j in np.linspace(0,n-1,n): - wave_dirac[j,:,:] = scp.fftconvolve(wave_dirac[j,:,:],PSF,mode = 'same') - - lensed_wave = np.zeros((n,N1,N2)) - for i in range(n): - lensed_wave[i,:,:] = Lens.image_to_source(wave_dirac[i,:,:],x,y,N1,N2,Fout,lensed=lensed)**2 - wave_sum = np.sum(np.sum(lensed_wave,1),1) - sig = np.sqrt(wave_sum) - return sig - -def lens_one(beta, nt1,nt2,size, ones = 0): - # n1,n2: size of the postage stamp in image plane - # N1,N2: size of the postage stamp in source plane +def lens_one(Fkappa, nt1,nt2,size): + ##DESCRIPTION: + ## Function that maps an all at one image to source plane. + ## + ##INPUTS: + ## -Fkappa: the mapping between source and image planes + ## -nt1,nt2: the shape of the image. + ## -size: the factor that scales the shape of the source relative to the shape of the image + ## + ##OUTPUTS: + ## -lensed: the projection to source plane of an all at aone image. dirac = np.ones((nt1,nt2)) - lensed = Lens.image_to_source(dirac, size,beta,lensed = [0])#multi = 1)# - if ones ==1: - lensed[np.where(lensed == 0)] =1 + lensed = Lens.image_to_source(dirac, size,Fkappa,lensed = [0]) return lensed -def unlens_one(F, xt,yt,nt1,nt2,nb1,nb2, ones = 0): - # n1,n2: size of the postage stamp in image plane - # N1,N2: size of the postage stamp in source plane - dirac = np.ones((nb1,nb2)) - unlensed = Lens.source_to_image(dirac,nt1, nt1,F) - if ones ==1: - unlensed[np.where(lensed == 0)] =1 - return unlensed - -def MOM(RG,RS,levelg,levels, lvlg,lvls, sigma): - nt1,nt2 = np.shape(RG) - ns1,ns2 = np.shape(RS) - levels[levels==0]=1 - wG = mw.wave_transform(RG,lvlg, newwave = 1)/levelg/sigma - wS = mw.wave_transform(RS,lvls, newwave = 1)/levels/sigma - - wGmax = np.max(wG) - wSmax = np.max(wS) - wmax = [wGmax, wSmax] - - k = np.min(wmax)+(max(wmax)-min(wmax))/100 - return k - - -def level(nt1,nt2,lvl = 6): - dirac = np.zeros((nt1,nt2)) - dirac[nt1/2,nt2/2] = 1 - wave_dirac = mw.wave_transform(dirac,lvl, newwave = 0) - - wave_sum = np.sqrt(np.sum(np.sum(wave_dirac**2,1),1)) - levels = np.multiply(np.ones((lvl,nt1,nt2)).T,wave_sum).T - - return levels - -def linorm(A,nit): - ns,nb = np.shape(A) - x0 = np.random.rand(nb) - x0 = x0/np.sqrt(np.sum(x0**2)) - - for i in np.linspace(0,nit-1,nit): - x = np.dot(A,x0) - xn = np.sqrt(np.sum(x**2)) - xp = x/xn - y = np.dot(A.T,xp) - yn = np.sqrt(np.sum(y**2)) - if yn < np.dot(y.T,x0) : - break - x0 = y/yn - - return xn - -def MAD(x,n=3,fil=1): - if fil == 1: - meda = med.median_filter(x,size = (n,n)) - else: - meda = np.median(x) - medfil = np.abs(x-meda) - sh = np.shape(x) - sigma = 1.48*np.median((medfil)) - return sigma - - -def HT(alpha, k, levels, sigma, pen = 0, soft = 0): +def MAD(x,n=3): + ##DESCRIPTION: + ## Estimates the noise standard deviation from Median Absolute Deviation + ## + ##INPUTS: + ## -x: a 2D image for which we look for the noise levels. + ## + ##OPTIONS: + ## -n: size of the median filter. Default is 3. + ## + ##OUTPUTS: + ## -S: the source light profile. + ## -FS: the lensed version of the estimated source light profile + meda = med.median_filter(x,size = (n,n)) + medfil = np.abs(x-meda) + sh = np.shape(x) + sigma = 1.48*np.median((medfil)) + return sigma + + +def ST(alpha, k, levels, sigma): + ##DESCRIPTION: + ## Soft thresholding operator. + ## + ##INPUTS: + ## -alpha: the starlet decomposition to be thresholded. + ## -k: the threshold in units of noise levels (usually 5). + ## -levels: the noise levels at each scale and location of the starlet decomposition. + ## -sigma: the noise standard deviation. + ## + ##OUTPUTS: + ## -alpha: The thresholded coefficients. lvl, n1,n2 = np.shape(alpha) th = np.ones((lvl,n1,n2))*k th[0,:,:] = th[0,:,:]+1 @@ -320,19 +463,30 @@ def HT(alpha, k, levels, sigma, pen = 0, soft = 0): alpha0 = np.copy(alpha) th = th*levels*sigma - if soft == 1: - alpha= np.sign(alpha0)*(np.abs(alpha0)-th) - alpha[np.where(np.abs(alpha)-th<0)]=0 - else: - alpha[np.where(np.abs(alpha)-th<0)] = 0 + alpha= np.sign(alpha0)*(np.abs(alpha0)-th) + alpha[np.where(np.abs(alpha)-th<0)]=0 - return alpha + return alpha -## -############################ SLIT MCA FISTA -def simulate_noise(nt1,nt2, size, beta, lensed, PSFconj): +def simulate_noise(nt1,nt2, size, Fkappa, lensed, PSFconj): + ##DESCRIPTION: + ## Simulates noise levels in source plane from lensing operator and convolution operator. + ## + ##INPUTS: + ## -nt1,nt2: the shape of the images for which to simulate noise maps. + ## -size: scaling factor for the shape of the source. + ## -Fkappa: Projection operator between lens and source plane. + ## -lensed: mapping of an all at one image to source plane. + ## -PSFconj: the conjugate of the PSF + ## + ##OPTIONS: + ## -n: size of the median filter. Default is 3. + ## + ##OUTPUTS: + ## -S: the source light profile. + ## -FS: the lensed version of the estimated source light profile n = 500. nb1,nb2 = nt1*size, nt2*size lvl = np.int(np.log2(nb1)) @@ -349,188 +503,12 @@ def simulate_noise(nt1,nt2, size, beta, lensed, PSFconj): for i in range(np.int(n)): noise = np.random.randn(nt1,nt2) noise = scp.fftconvolve(noise,PSFconj,mode = 'same') - noise_lens = (Lens.image_to_source(noise, size, beta, lensed =lensed)) + noise_lens = (Lens.image_to_source(noise, size, Fkappa, lensed =lensed)) noise_lens[noise_lens ==0] = 1 - storage[:,:,:,i] = mw.wave_transform(noise_lens, lvl = lvl) - w_levels = np.std(storage, axis = 3) - return w_levels -def SLIT_MCA(img, theta, kmax, niter, riter, kinter, ninter, size,decrease = 0, FB = 0, - pos = True, K=0, posit = False,reweighting = False, soft = 0, mask = [0,0], - PSF = [0], PSFconj = [0,0], gauss = 0, Ginit=0, repeat = 1, mrfilter = 0): - - beta = theta - #Initialisation of the lensed sourced - Im = np.copy(img)*0. - - nt1,nt2 = np.shape(img) - #Initialisation of the source - nb1= nt1*size - nb2 = nt2*size - - lvlg = np.int(np.log2(nt1)) - lvls = np.int(np.log2(nb1)) - #Masking if required - if np.sum(mask) == 0: - mask = np.ones((nt1,nt2)) - img = img*mask - M = [0] - #Noise in image plane - sigma0 = wine.MCA.MAD(img) - #Mapping of an all-at-one image - lensed = lens_one(beta, nt1,nt2, size) - bound = lensed/lensed - bound[lensed == 0] = 0 - #Noise level map with convolution by wavelet elements +weird normalisation - debound = mw.wave_transform(1-bound, lvls) - - levelg = level(nt1,nt1,lvl=lvlg) - - levels = simulate_noise(nt1,nt2, size, beta, lensed, PSFconj) - levels0 = lens_wnoise(nb1,nb2,size,beta, lvl = lvls) - - S =Lens.image_to_source(img, size, beta, lensed = lensed) - - #K thresholds in units of sigma - if np.sum(Ginit) == 0: - Ginit = np.random.randn(nt1,nt2)*sigma0 - alphag = mw.wave_transform(Ginit, lvlg, newwave = 1) - wFS = np.ones((lvlg,nt1,nt2)) - RG = np.copy(img) - RS = np.copy(S) - if np.sum(PSF) != 0: - RS = scp.fftconvolve(RS,(PSFconj),mode = 'same') - Res1 = [] - Res2 = [] - R=0 - i = 0 - S=np.random.randn(nb1,nb2)*sigma0 - G = Ginit - FS = 0 - k = 0 - alphas = np.zeros((lvls,nb1,nb2)) - csis = np.copy(alphas) - csig = np.copy(alphag) - K_s = np.zeros((niter)) - norm = np.zeros((3,niter)) - tg =1 - -## -##Compute spectral norms - def Finv_apply(I): - return Lens.image_to_source(I, size, beta, lensed = lensed) - def F_apply(Si): - return Lens.source_to_image(Si, nt1, nt2,theta) - def PSF_apply(i): - return scp.fftconvolve(i,PSF,mode = 'same') - def PSFT_apply(ii): - return scp.fftconvolve(ii,PSFconj,mode = 'same') - def star(x): - return mw.wave_transform(x, lvlg, newwave = 1) - def star_inv(x): - return mw.iuwt(x) - - F_norm = spectralNorm(nt1,nt2,20,1e-10,F_apply,Finv_apply) - Star_norm_im = spectralNorm(nt1,nt2,20,1e-10,star,star_inv) - Star_norm_s = spectralNorm(nb1,nb2,20,1e-10,star,star_inv) - if np.sum(PSF)>0: - PSF_norm = spectralNorm(nt1,nt2,20,0.0000001,PSF_apply,PSFT_apply) - muG = 1./(Star_norm_im**2) - muS = 1./(Star_norm_s*F_norm*PSF_norm)**2 - else: - muG = 1./(Star_norm_im**2) - muS = 1./(Star_norm_s*F_norm)**2 - - mu = 2./(Star_norm_s**2+F_norm**2) - - def HGT(X): - return X - def HG(X): - return X - def Phi(X): - return mw.iuwt(X) - def PhiT(X): - return mw.wave_transform(X, lvlg, newwave = 1) - - def HFT(X): - return Lens.image_to_source(X, size, beta, lensed=lensed) - def HF(X): - return Lens.source_to_image(X, nt1,nt2, theta) - - RS = Lens.image_to_source(img, size, beta, lensed = lensed) - RG = np.copy(img) - G = np.random.randn(nt1,nt2)*sigma0 - S = np.random.randn(nb1,nb2)*sigma0 - FS = np.random.randn(nt1,nt2)*sigma0 - i=0 - k = kmax - while i < niter: - print(i) - DS = img-G - j = 0 - ts = 1 - tr = np.zeros(riter) - - tauG = 1 - tauS = 1 - - while j < riter: - RS = muS*(DS-FS) - if np.sum(PSF)!=0: - RS = scp.fftconvolve(RS,(PSFconj),mode = 'same') - RS = Lens.image_to_source(RS, size, beta, lensed=lensed) - - alphas_new = csis+mw.wave_transform(RS, lvls, newwave = 1) - - alphas_new = HT(alphas_new, k, levels, sigma0, soft = soft)*bound - ts_new = (1.+np.sqrt(1.+4.*ts**2))/2. - csis = alphas_new+((ts-1)/ts_new)*(alphas_new-alphas) - alphas = alphas_new - ts = np.copy(ts_new) - S = mw.iuwt(csis) - - FS = Lens.source_to_image(S, nt1, nt2,theta) - - - if np.sum(PSF)!=0: - FS = scp.fftconvolve(FS,PSF,mode = 'same') - tr[j] = (np.std(img-FS-G)**2)/sigma0**2 - j+=1 - S = mw.iuwt(alphas) - FS = Lens.source_to_image(S, nt1, nt2,theta) - if np.sum(PSF)!=0: - FS = scp.fftconvolve(FS,PSF,mode = 'same') - - DG = img-FS - RG = muG*(DG-G) - - alphag_new = csig+mw.wave_transform(RG, lvlg, newwave = 1) - alphag_new = HT(alphag_new, kmax, levelg, sigma0, soft = soft) - tg_new = (1.+np.sqrt(1.+4.*tg**2))/2. - csig = alphag_new+((tg-1)/tg_new)*(alphag_new-alphag) - alphag = alphag_new - - tg = np.copy(tg_new) - G = mw.iuwt(csig) - - newres = (np.std(img-FS-G)**2)/sigma0**2 - K_s[i] = newres - - res = np.copy(newres) - - i = i+1 - - - S[np.where(S<0)] = 0 - FS = Lens.source_to_image(S, nt1, nt2,theta) - if np.sum(PSF)!=0: - FS = scp.fftconvolve(FS,PSF,mode = 'same') - - return S, FS,G - diff --git a/dist/SLIT-0.1-py2.7.egg b/dist/SLIT-0.1-py2.7.egg index 52b1e00..586f40b 100644 Binary files a/dist/SLIT-0.1-py2.7.egg and b/dist/SLIT-0.1-py2.7.egg differ