; $Id: decomyro_gpa.pro,v 1.5 2012/02/10 17:15:23 lmugnier Exp $ ;+ ;NOM : DECOMYRO_GPA.PRO ; ;SYNTAXE : ; DECOMYRO_GPA, images, wi_in, param, TERMREG_OBJ = termreg_obj, $ ; TERMREG_PSF = termreg_psf, VISUIM=visuim,$ ; VISUPSF=visupsf, OBJGUESS=objguess, PSFGUESS=psfguess,$ ; OBJREF=objref, POSITIVITE=positivite, $ ; POISSON=poisson,NORMOBJGUESS=normobjguess, $ ; CONTINUE=continue,REGPSF=regpsf, $ ; HYPER_OBJ = hyper_obj, HYPER_PSF = hyper_psf ; ITERATION=iteration,$ ; NORMPSF = normpsf, PSFCSTE=psfcste,SAVE1IT = save1it, SAVEIT = $ ; saveit, DIRECTORY= directory, NOMFICH = nomfich, RENORM = $ ; renorm, INFO = info, MOYOBJ = moyobj, AVGPSF = avgpsf, ; SEUIL=seuil, ; QUIET=quiet,CONV=conv,PERIODE=periode,LOG=log, SFIELD=sfield ; ; ;DESCRIPTION : Deconvolution myope a partir des donnees sur PSF (moyenne et ; DSP) utilisation de fmin_GPA ; - contrainte de positivite stricte sur l'objet ; - prise en compte d'un bruit a statistique gaussienne uniforme ; (bruit de detecteur) ou poissonienne (bruit de photon). ; - Hyper_obj : valeur de l'hyper parametre de regularisation ; objet ; - Hyper_psf : valeur de l'hyper parametre de regularisation ; PSF (1 par défaut) ; - regularisation de l'objet par ajout dans l'erreur de somme de ; |O-Om|^2/DSP_objet (ou carre de Laplacien), ou Lp, ou L1-L2 ; (parametre de regularisation: HYPER) ; - regularisation ou non de la psf par ajout dans l'erreur de ; somme de |h-hm|^2/DSP_psf (ou carre de Laplacien) avec meme ; parametre de regularisation (HYPER) que pour objet (a changer?) ; - TERMREG_OBJ permet de donner une régularisation sur l'objet : ; si TERMREG_OBJ est un tableau : régularisation quadratique, ; avec ce tableau comme estimation de la DSP de l'objet, et ; MOYOBJ comme objet moyen. ; si TERMREG_OBJ est un scalaire : ; TERMREG_OBJ = 1 : penalisation quadratique pour faibles ; ecarts et lineaire pour les grands écarts (cf SEUIL) ; TERMREG_OBJ = 2 : penalisation quadratique pour tous les écarts ; TERMREG_OBJ entre 1 et 2 : pénalisation L_p, p = termreg_obj ; - SEUIL est un seuil sur le gradient objet qui permet de régler ; le passage de la régularisation L2 à L1 (utilisé seulement ; pour TERMREG_OBJ = 1). ; - TERMREG_PSF permet de donner une estimation de la DSP de la PSF ; comme fonction regularisante ; - normalisation de la psf pour le calcul du gradient (NORMPSF) ; - PSFCSTE pour ne minimiser l'erreur que vis-a-vis des ; parametres objet ; - CONTINUE permet de poursuivre le calcul apres ; interruption => ne pas reinitialiser les parametres ; - SAVE1IT permet la sauvegarde des parametres à une itération ; donnée ; - SAVEIT permet la sauvegarde des parametres toutes les n ; itérations (SAVEIT = n) ; - DIRECTORY donne le directory où sauver les param de saveit ; - NOMFICH donne le nom generique de fichier où sauver les param ; - INFO sauve dans un tableau les valeurs de l'err a l'iteration ; d'arret, la valeur de l'hyper et la valeur du coef de ; renormalisation (à appliquer dans param2obj) ; - VISUIM et VISUPSF permettent un affichage graphique de ; l'objet et/ou la PSF pdt convergeance. ; - OBJREF mesure d'une distance a l'objet vrai si connu (=> en ; simulation) ; - QUIET evite les messages à chaque itération donc les .log ; énormes. ; -CONV critere d'arret ; -LOG affichage en echelle Log ; - critere L1L2 adapte au champs d'etoile : on ; s'interesse au pixel et pas a leur derivees ;; CRITERE_KEYBORD : (sortie) variable de sorie valant 0 SAUF si arret par 'Q' ; au clavier, dans ce cas la ca vaut 1 !!! Utile uniquement ; pour Mistral en Myope ;DIAGNOSTIC D'ERREUR : ; ; ATTENTION: le programme ne peut traiter que des images carrees !!!! ; ; ATTENTION: problème à fort flux lorsque l'on est en statistique ; Poissonienne (/poi). Il faut prendre une statistique Gaussienne ; (/pos) pour avoir une bonne convergence <- Probleme resulu par ; la double precision mais de toutes facons, une statitique ; Gaussienne non uniforme sera bientot implantee pour le bruit ; !!! ; ;VOIR AUSSI : ; decomyro.pro, decostars.pro, decosf.pro ; ;HISTORIQUE : ; $Log: decomyro_gpa.pro,v $ ; Revision 1.5 2012/02/10 17:15:23 lmugnier ; Correction du plantage avec /VMLM (errx pas defini). ; ; Revision 1.4 2012/02/10 17:06:28 lmugnier ; Mot-clé /VMLM (et modifs sur positivité par projection), ; check-in de la version par TF restée "checked-out" quelque temps... ; ; Revision 1.3 2003/04/16 12:05:17 fusco ; Rajout de la positivite par projection. Mot cle pos = 2 pour projection (pos=1 ; pour reparam) ; ; Revision 1.2 2003-02-18 19:15:37+01 fusco ; *** empty log message *** ; ; Revision 1.1 2000-10-13 08:28:20+02 fusco ; Initial revision ; ; Revision 1.14 1999-11-09 18:08:18+01 fusco ; *** empty log message *** ; ; Revision 1.13 1999-11-09 15:49:40+01 fusco ; *** empty log message *** ; ; Revision 1.12 1999-11-09 15:47:54+01 fusco ; *** empty log message *** ; ; Revision 1.11 1999-10-21 11:34:46+02 fusco ; Rajout de la regularisation on uniforme, pour le moment uniquement dans le ; cadre L1-L2 ou Lp ; ; Revision 1.10 1999-10-19 16:59:03+02 fusco ; *** empty log message *** ; ; Revision 1.9 1999-10-12 18:44:02+02 fusco ; *** empty log message *** ; ; Revision 1.8 1999-10-12 17:26:14+02 fusco ; *** empty log message *** ; ; Revision 1.7 1999-09-25 11:44:15+02 fusco ; *** empty log message *** ; ; Revision 1.6 1999-09-25 11:39:20+02 fusco ; *** empty log message *** ; ; Revision 1.5 1999-09-20 19:19:56+02 fusco ; rajout de la statistique de bruit gaussienne non uniforme. ; Activeé si wi = tableau et non entier ; ; Revision 1.4 1999-05-19 09:27:32+02 fusco ; *** empty log message *** ; ; Revision 1.3 1999-05-18 17:13:04+02 fusco ; *** empty log message *** ; ; Revision 1.2 1999-05-18 09:24:41+02 fusco ; rajout du mot cle periode ; ; Revision 1.1 1999-05-18 09:11:44+02 fusco ; Initial revision ; ;- FUNCTION convergence,old,new ;estimee indicative de l'evolution de la convergence conv = sqrt(total((new-old)^2) / total(new^2)) old = new return ,conv END ; FUNCTION param2psf_mistral ,param,jj ;passage du vecteur parametre a la psf numero jj ; COMMON JMCDEC, images_ou_tf, wi, status ; dim = status.dim ; dim2 = dim*dim ; deb = dim2 + jj*dim2 ; psf = dblarr(dim , dim) ; IF keyword_set(status.pos) THEN psf(*) = param(deb:deb + dim2-1)^2. $ ;cas positivite stricte ; ELSE psf(*) = param(deb:deb + dim2-1) ; ;normalisation tq tf_psf(0,0) = sqrt(fluxtot)/dim/status.norm ; IF (status.normpsf EQ 1B) THEN BEGIN ; cste = sqrt(status.flux)/float(dim)/status.norm*float(dim2) ; psf = psf/total(psf) * cste ; END ; return , psf ; END FUNCTION erreur_grad, param, grad, ERR2 = err2 ;calcul de l'erreur et du gradient (si necessaire) ;[renvois en option la 2ieme composante de l'erreur] ;[option valide que lors du calcul de l'erreur seule] COMMON JMCDEC, images_ou_tf, wi, status COMMON DECONV_REGUL, termregobj, termregpsf, mtfobj, mfto, seuil, hyperobj dim = status.dim dim2 = dim*dim ;constante de normalisation de la psf [cf. fonction param2psf_mistral] cste = sqrt(status.flux)/float(dim)/status.norm*float(dim2) dsp_ou_p = termregobj dsp_psf = termregpsf obj = param2obj(param,DIM=status.dim,POS=status.pos); domaine reel centre' tf_obj = fft( eclat( obj ), -1) ;############################ ; CALCUL ERREUR SEULE ;############################ IF n_params() EQ 1 THEN BEGIN err1 = 0. err2 = 0. err3 = 0. deb = dim2 FOR jj=0,status.nbimages-1 DO BEGIN tf_psf = fft( eclat( param2psf_mistral(param,jj) ), -1) ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;PREMIERE COMPOSANTE DE L'ERREUR [I-OS ou equivalent] ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;CAS STATISTIQUE NON UNIFORME (bruit gaussien non uniforme) ;IF (status.pois EQ 1B) THEN BEGIN IF (status.nonstat EQ 1B) THEN BEGIN gg = (real(fft(tf_obj*tf_psf,1))) err1 = err1 + total(wi*abs(images_ou_tf(*,*,jj)-gg)^2) ;CAS STATISTIQUE POISSONNIENE (bruit de photon) ;err1 = err1 + total((gg > 0.)- images_ou_tf(*,*,jj)*alog(gg > $ ; 0.)) ;=> gg > 0. positivite stricte ;Rq: images_ou_tf = images dans cas poisson END ELSE BEGIN ;CAS BRUIT GAUSSIEN UNIFORME (bruit de ;detecteur) : dans ce cas on calcule l'erreur ; directement en Fourier err1 = err1 + dim2*total( wi * abs(images_ou_tf(*,*,jj) - $ tf_obj*tf_psf)^2) ;le facteur dim2 provient de Parceval en discret et assure que: ; dim2*total( abs(images_ou_tf(*,*,jj) - tf_obj*tf_psf)^2) = total(images - gg) ; ou gg = obj convolue a psf = eclat(float(fft(tf_obj*tf_psf,1))) END ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;DEUXIEME COMPOSANTE DE L'ERREUR [contraintes sur ai OU regularisation psf] ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% IF NOT(status.psfcste) THEN BEGIN IF ((status.hyper_psf NE 0.)AND(status.regpsf NE 0)) THEN BEGIN ;regularisation FEP IF (total(where(dsp_psf NE 0))EQ -1.) THEN BEGIN err2 = status.hyper_psf *dim2/status.flux* total(abs((tf_psf-mfto) * distc(dim)^2 )^2 ) END ELSE err2 = status.hyper_psf *dim2/status.flux* total(abs(tf_psf-mfto)^2/dsp_psf) END END END ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;TROISIEME COMPOSANTE DE L'ERREUR ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;[regularisation par carre du laplacien soit ;ds domaine Fourier: ;tf_obj*freq_spatiale^2] IF (status.hyper_obj NE 0.) THEN BEGIN IF (n_elements(dsp_ou_p) GT 1) THEN BEGIN ; DSP err3 = hyperobj*dim2/status.flux* total(abs(tf_obj-mtfobj)^2/dsp_ou_p) END ELSE IF (dsp_ou_p EQ 0.) THEN $ ; Laplacien err3 = hyperobj*dim2/status.flux* total(abs((tf_obj-mtfobj) * distc(dim)^2 )^2 ) $ ELSE IF (dsp_ou_p GT 1.) THEN BEGIN ; norme Lp (pas utilisee... on ; pourrait la virer) IF status.sfield EQ 1B THEN BEGIN ;cas de champs d'etoiles : ;!!!! PAS ENCORE COMPLETEMENT DEBUGUE !!!! ;decorrelation complete entre pixels de l'obj => L1 L2 ;sur obj et non gradient => DSP = cste Dx = (obj)/seuil Dy = 0. ;pour rester coherent avec la suite et eviter de ;modifier les lignes suivantes ENDIF ELSE BEGIN Dx = obj - shift(obj, 1, 0) Dy = obj - shift(obj, 0, 1) ENDELSE err3 = dim2/status.flux * $ total(hyperobj*(Dx^2 + Dy^2)^(dsp_ou_p/2)) ENDIF ELSE BEGIN ; L1-L2 IF status.sfield EQ 1B THEN BEGIN ;cas de champs d'etoiles : ;!!!! PAS ENCORE COMPLETEMENT DEBUGUE !!!! ;decorrelation complete entre pixels de l'obj => L1 L2 ;sur obj et non gradient => DSP = cste Dx = (obj)/seuil Dy = 0. ;pour rester coherent avec la suite et eviter de ;modifier les lignes suivantes ENDIF ELSE BEGIN Dx = (obj - shift(obj, 1, 0))/seuil Dy = (obj - shift(obj, 0, 1))/seuil ENDELSE IF (dsp_ou_p EQ -1.) THEN BEGIN err3 = (dim2/status.flux*2.*seuil^2) * $ $ (total(hyperobj*(abs(Dx)-alog(abs(Dx)+1.)))+total(hyperobj*(abs(Dy)-alog(abs(Dy)+1.)))) END ELSE BEGIN err3 = (dim2/status.flux*2.*seuil^2) * $ $ (total(hyperobj*(sqrt(Dx^2+Dy^2)-alog(sqrt(Dx^2+Dy^2)+1.)))) ENDELSE ENDELSE END ;sommation des differentes composantes err = err1 + err2 + err3 END ELSE BEGIN ;############################ ; CALCUL ERREUR ET GRADIENT ;############################ err = 0. ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;PREMIERE & DEUXIEME COMPOSANTE => BOUCLE SUR IMAGES ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% deb = dim2 grad = param*0. grad_obj = dcomplexarr(dim,dim) t2 = dcomplexarr(1) t1 = dcomplexarr(1) FOR jj=0,status.nbimages-1 DO BEGIN psf = param2psf_mistral(param,jj) ;sera utile pour le calcul du gradient tf_psf = fft( eclat( psf ), -1) ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;PREMIERE COMPOSANTE DE L'ERREUR... ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;CAS STATISTIQUE NON UNIFORME (bruit de photon+detecteur) IF (status.nonstat EQ 1B) THEN BEGIN ;ERREUR (je rappelle que c'est seulement la premiere ;composante de l'erreur) gg = (real(fft(tf_obj*tf_psf,1))) err = err + total(wi*abs(images_ou_tf(*,*,jj)-gg)^2) ;CAS STATISTIQUE POISSONNIENE (bruit de photon) ;err = err + total( (gg > 0.)- images_ou_tf(*,*,jj)*alog(gg > 0.) ) ;=> gg > 0. positivite stricte, que faire si g(r) = 0. ???? ;Rq: images_ou_tf = images dans cas poisson ;...ET SON GRADIENT: dd = -2*wi*(images_ou_tf(*,*,jj)-gg) ;CAS STATISTIQUE POISSONNIENE (bruit de photon) ;dd = 1. - images_ou_tf(*,*,jj)/(gg > 0.) ;que faire si g(r) = 0. ???? tf_dd = fft(( dd ), -1) ;bufferisation pour calcul du GRADIENT/obj: derr/dobj grad_obj = grad_obj + tf_dd * conj(tf_psf) ;gradient/psf: derr/ds ;d(err)/d(psf) = derr/ds if not (status.psfcste) then grad_psf = tf_dd * conj(tf_obj) END ELSE BEGIN ;CAS BRUIT GAUSSIEN ADDITIF UNIFORME (bruit de detecteur) ;ERREUR (je rappelle que c'est seulement la premiere composante de l'erreur) err = err + dim2*total( wi * abs(images_ou_tf(*,*,jj) - tf_obj*tf_psf)^2) ;...ET SON GRADIENT: ;bufferisation pour calcul du GRADIENT/obj: derr/dobj grad_obj = grad_obj + 2. * wi * ( tf_obj*tf_psf - $ images_ou_tf(*,*,jj) ) * conj(tf_psf) ;gradient/psf: derr/ds (calcul en deux etapes: derr/dS -> derr/ds) ; en toute rigueur ces 2 etapes s'ecrivent: (...) ;d(err)/d(tf_psf) = derr/dS ;grad_psf = 2. * dim2 * wi * (tf_psf*tf_obj - images_ou_tf(*,*,jj)) * conj(tf_obj) ;d(err)/d(psf) = derr/ds ;grad_psf = 1./float(dim2) * eclat( float( fft(grad_psf, 1) ) ) ; (...) mais en shintant le *dim2 suivi du /dim2 on ecrit: if not (status.psfcste) then grad_psf = 2. * wi * $ ( tf_obj*tf_psf - images_ou_tf(*,*,jj)) * conj(tf_obj) END ;FIN DES CALCULS LIES A LA PREMIERE PARTIE DE L'ERREUR ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;DEUXIEME PARTIE DE L'ERREUR ... [regularisation psf] ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% IF NOT(status.psfcste) THEN BEGIN IF (status.hyper_psf NE 0.)AND(status.regpsf NE 0) THEN BEGIN ;ERREUR ... [cas regularisation psf] IF (total(where(dsp_psf NE 0))EQ -1.) THEN BEGIN err = err + status.hyper_psf *dim2/status.flux* $ $ total( abs((tf_psf-mfto) * distc(dim)^2 )^2 ) ;... ET SON GRADIENT (ajouter au gradient du premier ;terme precedemment calcule) grad_psf = grad_psf + 2. * status.hyper_psf*dim2/status.flux * (tf_psf-mfto) * $ distc(dim)^4/dim2 ENDIF ELSE BEGIN err = err + status.hyper_psf*dim2/status.flux * total( abs(tf_psf-mfto)^2/dsp_psf ) ;... ET SON GRADIENT (ajouter au gradient du premier ;terme precedemment calcule) grad_psf = grad_psf + 2. * status.hyper_psf*dim2/status.flux * $ (tf_psf-mfto)/(dsp_psf*dim2) ENDELSE ENDIF ;GRADIENT/PARAMETRES [cf. fonction param2psf_mistral] et Paragraphe 4.B du Thiebaut/Conan 1995 ;+ STOCKAGE DU RESULTAT DANS LE VECTEUR GRADIENT grad_psf = eclat( float( fft(grad_psf, 1) ) ) ;vecteur parametre courant paramjj = param(deb : deb + dim2-1) IF (status.pos EQ 1B) THEN BEGIN IF (status.normpsf EQ 1B) THEN $ grad(deb : deb + dim2-1) = 2. *paramjj/ total(paramjj^2) * $ (cste*grad_psf(*)-total(psf(*)*grad_psf(*)))$ ELSE grad(deb : deb + dim2-1) = 2.*paramjj*grad_psf(*) END ELSE BEGIN IF (status.normpsf EQ 1B) THEN $ grad(deb : deb + dim2-1) = (cste*grad_psf(*) - total(psf(*)*grad_psf(*))) / total(paramjj) $ ELSE grad(deb : deb + dim2-1) = grad_psf(*) END deb = deb + dim2 END END ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;TROISIEME COMPOSANTE DE L'ERREUR : REGULARISATION OBJET ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% IF (status.hyper_obj NE 0.) THEN BEGIN IF (n_elements(dsp_ou_p) GT 1) THEN BEGIN ; DSP err = err + hyperobj*dim2/status.flux * total(abs(tf_obj-mtfobj)^2/dsp_ou_p ) ;[regularisation par DSP] ;grad_obj_ap = 2. * hyperobj*dim2/status.flux * (tf_obj-mtfobj)/(dsp_ou_p*dim2) ;grad_obj = grad_obj + grad_obj_ap grad_obj = grad_obj + 2. * hyperobj*dim2/status.flux * (tf_obj-mtfobj)/(dsp_ou_p*dim2) grad_obj = eclat(float( fft(grad_obj, 1) ) ) ;retour domaine reel ;centre' BOFF !!! il y a un eclat qui sert a ;rien mais de toutes façon j'utilise plus la DSP ENDIF ELSE IF (dsp_ou_p EQ 0.) THEN BEGIN ; Laplacien err = err + hyperobj*dim2/status.flux * total(abs((tf_obj-mtfobj) * distc(dim)^2 )^2 ) ;[regularisation par carre du laplacien] grad_obj = grad_obj + 2. * hyperobj*dim2/status.flux * (tf_obj-mtfobj) * $ distc(dim)^4/dim2 grad_obj = eclat( float( fft(grad_obj, 1) ) );retour domaine reel centre' ENDIF ELSE IF (dsp_ou_p GT 1.) THEN BEGIN ; [norme Lp] grad_obj = eclat( float( fft(grad_obj, 1) ) );retour domaine reel centre' IF status.sfield EQ 1B THEN BEGIN ;cas de champs d'etoiles : ;decorrelation complete entre pixels de l'obj => L1 L2 ;sur obj et non gradient => DSP = cste Dx = (obj)/seuil Dy = 0. ;pour rester coherent avec la suite et eviter de ;modifier les lignes suivantes ENDIF ELSE BEGIN Dx = obj - shift(obj, 1, 0) Dy = obj - shift(obj, 0, 1) ENDELSE ; print, 'erreur donnees', err err = err +dim2/status.flux * total(hyperobj*(Dx^2 + Dy^2)^(dsp_ou_p/2)) ; print, 'erreur totale', err epsilon = 1e-5 temp = Dx / ((Dx^2 + Dy^2) > epsilon)^(1.-dsp_ou_p/2.) ; grad_obj = grad_obj + status.hyper_obj*dim2/status.flux * dsp_ou_p*(temp - shift(temp, -1, 0)) grad_obj_ap = hyperobj*dim2/status.flux * dsp_ou_p*(temp - shift(temp, -1, 0)) temp = Dy / ((Dx^2 + Dy^2) > epsilon)^(1.-dsp_ou_p/2.) ; grad_obj = grad_obj + status.hyper_obj*dim2/status.flux * dsp_ou_p*(temp - shift(temp, 0, -1)) grad_obj_ap = grad_obj_ap + hyperobj*dim2/status.flux * dsp_ou_p*(temp - shift(temp, 0, -1)) grad_obj = grad_obj + grad_obj_ap ; o en sqrt(N) ENDIF ELSE BEGIN ; [norme L2-L1] grad_obj = eclat( float( fft(grad_obj, 1) ) );retour domaine reel centre' IF status.sfield EQ 1B THEN BEGIN ;cas de champs d'etoiles : ;decorrelation complete entre pixels de l'obj => L1 L2 ;sur obj et non gradient => DSP = cste Dx = (obj)/seuil Dy = 0. ;pour rester coherent avec la suite et eviter de ;modifier les lignes suivantes ENDIF ELSE BEGIN Dx = (obj - shift(obj, 1, 0))/seuil Dy = (obj - shift(obj, 0, 1))/seuil ENDELSE IF (dsp_ou_p EQ -1.) THEN BEGIN IF status.sfield EQ 1B THEN BEGIN ;cas de champs d'etoiles : ;decorrelation complete entre pixels de l'obj => L1 L2 ;sur obj et non gradient => DSP = cste Dx = (obj)/seuil Dy = 0. ;pour rester coherent avec la suite et eviter de ;modifier les lignes suivantes ENDIF ELSE BEGIN Dx = (obj - shift(obj, 1, 0))/seuil Dy = (obj - shift(obj, 0, 1))/seuil ENDELSE err = err + (dim2/status.flux*2.*seuil^2) * $ (total(hyperobj*(abs(Dx)-alog(abs(Dx)+1.))) + $ total(hyperobj*(abs(Dy)-alog(abs(Dy)+1.)))) temp = Dx / (abs(Dx) + 1.) grad_obj_ap = (hyperobj*dim2/status.flux*2.*seuil) * $ (temp-shift(temp, -1, 0)) temp = Dy / (abs(Dy) + 1.) grad_obj_ap = grad_obj_ap+(hyperobj*dim2/status.flux*2.*seuil) * $ (temp-shift(temp, 0, -1)) END ELSE BEGIN err = err + (dim2/status.flux*2.*seuil^2) * $ $ (total(hyperobj*(sqrt(Dx^2+Dy^2)-alog(sqrt(Dx^2+Dy^2)+1.)))) temp = Dx/(sqrt(Dx^2+Dy^2) + 1.) ;ouah ca s'arrange ;vachement bien cette derivée, c'est pas beau la science ? grad_obj_ap = (hyperobj*dim2/status.flux*2.*seuil) * $ (temp-shift(temp, -1, 0)) temp = Dy/(sqrt(Dx^2+Dy^2)+ 1.) grad_obj_ap = grad_obj_ap + $ (hyperobj*dim2/status.flux*2.*seuil) * $ (temp-shift(temp, 0, -1)) END grad_obj = grad_obj + grad_obj_ap ; o en sqrt(N) ENDELSE ENDIF ; sur status.hyper NE 0 ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;ON RANGE LE GRADIENT OBJET DANS LE VECTEUR GRADIENT/PARAMETRES_OBJET: ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;cas positivite stricte [particulier dans cas poisson] : IF (status.pos EQ 1B) THEN grad(0:dim2-1) = 2. * param(0:dim2-1) * grad_obj $ ELSE grad(0:dim2-1) = grad_obj ENDELSE ; fin du calcul erreur et gradient return , err END PRO decomyro_gpa, images, wi_in, param, TERMREG_OBJ = termreg_obj, $ TERMREG_PSF = termreg_psf, VISUIM=visuim,$ VISUPSF=visupsf, OBJGUESS=objguess, PSFGUESS=psfguess,$ OBJREF=objref, POSITIVITE=positivite, $ POISSON=poisson,NORMOBJGUESS=normobjguess, $ CONTINUE=continue,REGPSF=regpsf, $ HYPER_OBJ = hyper_obj,HYPER_PSF = hyper_psf,$ ITERATION=iteration,$ NORMPSF = normpsf, PSFCSTE=psfcste,SAVE1IT = save1it, SAVEIT = $ saveit, DIRECTORY= directory, NOMFICH = nomfich, RENORM = $ renorm, INFO = info, MOYOBJ = moyobj, AVGPSF = avgpsf, $ SEUIL=seuil, QUIET=quiet, CONV = conv, PERIODE = periode, $ LOG = log, SFIELD = sfield, REGNOUNIF = regnounif, $ ERROR = error, VERROR = verror, $ CRITERE_KEYBORD = critere_keybord, VMLM = vmlm COMMON JMCDEC, images_ou_tf, wi, status COMMON DECONV_REGUL, termregobj,termregpsf, mtfobj, mfto, leseuil, hyperobj on_error,2 true = 1B false = 0B status = {struct_status, dim:0L, nbimages:0L, flux:0., norm:0.,$ hyper_psf:1., pos:0, $ pois:false, normpsf:false, regpsf:false, psfcste:false, $ nonstat:false, sfield:false, hyper_obj:0., regnounif:false} leseuil = 0. IF keyword_set(version) THEN $ message, "$Revision: 1.5 $, $Date: 2012/02/10 17:15:23 $", /INFO IF ((n_params() LT 3) OR (n_params() GT 5)) THEN $ message, $ "SYNTAXE : taper doc_library, 'decomyro'" if keyword_set(psfcste) then status.psfcste = true IF NOT(keyword_set(conv)) then conv = 1.e-10 IF (N_ELEMENTS(iteration) EQ 0) THEN BEGIN itmax = 1000 erreur = dblarr(itmax) END ELSE itmax = iteration IF NOT(keyword_set(PERIODE)) THEN periode = itmax critere_keybord = 0B sz = size(images) IF (sz(1) NE sz(2)) THEN message,'WARNING: decomyro ne fonctionne que sur des images carrees !!' status.dim = sz(1) IF (sz(0) EQ 2) THEN status.nbimages = 1 ELSE status.nbimages = sz(3) dim = status.dim dim2 = dim*dim IF (keyword_set(termreg_obj)) THEN BEGIN IF (n_elements(termreg_obj) GT 1) THEN $ termregobj = eclat(termreg_obj) $ ;regularisation par DSP ELSE BEGIN termreg_obj = float(termreg_obj) termregobj = termreg_obj IF ((termreg_obj EQ 1.0) OR (termreg_obj EQ -1.0)) THEN BEGIN IF NOT keyword_set(seuil) THEN $ message, "En régularisation L1-L2 il faut passer SEUIL" leseuil = seuil print, "regularisation L2-L1, seuil=", leseuil ENDIF ELSE print, "regularisation Lp, p=", termreg_obj ENDELSE ENDIF ELSE termregobj = 0 termregpsf = dblarr(dim, dim) IF (keyword_set(termreg_psf)) THEN termregpsf = eclat(termreg_psf) IF NOT keyword_set(objguess) THEN message,"mot cle OBJGUESS doit etre initialise" ;on pourrait s'en passer avec /CONTINUE mais vu que objguess donne status.flux ... ;rend compatible les normalisation de l'objet et l'image IF keyword_set(normobjguess) THEN BEGIN status.flux = total(images)/status.nbimages objguess = objguess/total(objguess)*status.flux END ELSE status.flux = total(objguess) ;suppose que objguess est normalise correctement au flux total IF keyword_set(regpsf) THEN status.regpsf = true IF keyword_set(SFIELD) THEN status.sfield = true IF keyword_set(normpsf)THEN status.normpsf = true IF positivite eq 0 THEN BEGIN status.pos = 0B active_set = 0B proj = 0 ENDIF IF positivite eq 1 THEN BEGIN status.pos = 1B active_set = 0B proj = 0 ENDIF IF positivite eq 2 THEN BEGIN status.pos = 0B proj = 1 ENDIF IF N_elements(hyper_psf) THEN status.hyper_psf = hyper_psf IF keyword_set(HYPER_OBJ) THEN BEGIN status.hyper_obj = 1B hyperobj = hyper_obj IF keyword_set(REGNOUNIF) THEN BEGIN status.regnounif = 1B IF n_elements(hyper_obj) NE n_elements(images) THEN message, $ 'Attention hyper obj doit etre un tableau de la taille de l''image' IF n_elements(dsp) GT 1 THEN message, $ $ 'Attention : regularisation adapative uniquement dans le cas L1-L2 ou Lp' ENDIF ENDIF IF keyword_set(poisson) THEN print, $ ;debug "Option POISSON choisie, on imposera donc aussi la POSITIVITE stricte sur l'objet" err_offset = 0. IF keyword_set(poisson) THEN BEGIN status.pois = true status.pos = 1B ;positivite obligatoire dans le cas d'un bruit de photon err_offset = total((.5*alog(2.*!pi) + (images+.5)*alog(images>1.)-(images>1.))*(images gt 1.)) ;estimation de somme/i[ alog(factorielle de images(*,*,i)) ] en utilisant Stirling ;terme constant de l'erreur n'intervenant pas dans la minimisation mais permettant d'obtenir ;a l'affichage une erreur positive END ;CAS D'UNE PSF DECRITE PIXEL A PIXEL IF NOT keyword_set(continue) THEN BEGIN IF (NOT keyword_set(psfguess)) THEN message,'PSFGUESS non initialise' $ ELSE BEGIN szzz = size(psfguess) IF (szzz(0) NE 2) AND (szzz(0) NE 3) THEN message,'PSFGUESS doit etre une image ou un cube d''images' IF (szzz(0) EQ 3) AND (szzz(3) NE status.nbimages) THEN message,'nb d''images dans PSFGUESS et IMAGES incompatible' IF (status.nbimages GT 1) AND (szzz(0) EQ 2) THEN print,'L''unique image PSFGUESS servira pour toutes les psf' IF (szzz(0) EQ 2) THEN nbpsfg = 1 ELSE nbpsfg = szzz(3) END END ;INITIALISATION des variables du COMMON JMCDEC ;IMAGES_OU_TF ;Rq: images_ou_tf = images dans cas poisson ; et = fft des images dans le cas bruit gaussien additif IF n_elements(wi_in) GT 1 THEN images_ou_tf = eclat(images) $ ELSE BEGIN images_ou_tf=dcomplexarr(dim,dim,status.nbimages) FOR jj=0,status.nbimages-1 DO images_ou_tf(*,*,jj) = fft( eclat( images(*,*,jj) ), -1) END ;WI (utile a priori que pour le bruit gaussien mais bon...) IF n_elements(wi_in) GT 1 THEN BEGIN status.nonstat = 1B wi_in = eclat(wi_in) end wi = wi_in ;QQUES CONSTANTES DE NORMALISATION PAS FORCEMENT ;INTELLEGIBLE MAIS IL FAUT ME CROIRE ;-) status.norm = 1. ;le facteur status.norm est necessaire pour eviter des pb numeriques de convergence ds fmin_conjgrad ;C'etait utile por la dec par aso. Je l'ai laisse au cas ou mais mis a 1. kkk = float(dim)/sqrt(status.flux) ;DEBUG - DEBUG kkk = kkk * status.norm ;tel que tf_obj(0) = sqrt(fluxtot)/dim*status.norm a la premiere iteration ;(l'energie de l'objet peut ensuite evoluer) ;l'energie de la psf, elle, est fixee une fois pour toute a sqrt(fluxtot)/dim/status.norm ;ce qui est compatible avec les images puisque: tf_images(0) ~ 1./dim2 * fluxtot ; ~ tf_obj(0)*tf_psf(0) leseuil = kkk * leseuil IF keyword_set(moyobj) THEN BEGIN moyobj = moyobj/total(moyobj)*status.flux*kkk mtfobj = fft(eclat(moyobj), -1) END ELSE mtfobj = dblarr(dim, dim) IF keyword_set(avgpsf) THEN BEGIN avgpsf = avgpsf/total(avgpsf)*float(dim)*sqrt(status.flux)/status.norm mfto = fft(eclat(avgpsf), -1) END ELSE mfto = dblarr(dim, dim) ;INITIALISATION DES PARAMETRES IF NOT keyword_set(continue) THEN BEGIN param = dblarr(dim2 + dim2*status.nbimages) ;PARAMETRES LIES A LA PSF FOR jj=0,status.nbimages-1 DO BEGIN kkkpsf = float(dim)*sqrt(status.flux)/status.norm IF (nbpsfg EQ 1) THEN guess = psfguess/total(psfguess)*kkkpsf ELSE $ guess = psfguess(*,*,jj)/total(psfguess(*,*,jj))*kkkpsf IF (status.pos EQ 0B) THEN BEGIN IF proj EQ 1 THEN param(dim2 + jj*dim2: dim2 + $ jj*dim2 + dim2-1) = guess > 0 $ ELSE param(dim2 + jj*dim2: dim2 + $ jj*dim2 + dim2-1) = guess ENDIF IF (status.pos EQ 1B) THEN param(dim2 + jj*dim2: dim2 + $ jj*dim2 + dim2-1) = sqrt(guess>0.) ENDFOR ;PARAMETRES LIES A L'OBJET ;avec ce qui suit deconv suppose que ;objguess est compatible avec images: ;exemple: si images en photons => objguess en photons guess = objguess*kkk IF (status.pos EQ 1B) THEN param[0:dim2-1] = sqrt(guess>0.) $ ELSE BEGIN IF (proj EQ 1) THEN param[0:dim2-1] = guess > 0. ELSE param[0:dim2-1] = guess ENDELSE objold = guess objnew = objold ENDIF ELSE BEGIN ;PARAMETRES LIES A L'OBJET objold = param2obj(param,DIM=status.dim,POS=status.pos) objnew = objold END ;JMC - DEBUG psfold = dblarr(dim,dim,status.nbimages) FOR jj=0, status.nbimages-1 DO psfold(*,*,jj) = param2psf_mistral(param,jj) psfnew = psfold ;initialisation de parametres pour le calcul d'une distance a l'objet vrai ;initialisation des affichages IF keyword_set(visuim) THEN window,3,xs=512,ys=256, $ xpos = 0, ypos = 770 ; BEGIN ; IF (status.dim GE 256) THEN window,3,xs=2*status.dim,ys=status.dim, $ ; xpos = 0, ypos = 770 ELSE window,3,xs=512,ys=256, $ ; xpos = 0, ypos = 770 ; END ;pour l'instant affichage de la psf uniquement si 1 image !!! IF keyword_set(visupsf) THEN BEGIN IF (status.nbimages EQ 1) THEN BEGIN window,4,xs=512,ys=256, xpos = 680, ypos = 770 ; IF (status.dim lt 256) THEN window,4,xs=2*status.dim,$ ; ys=status.dim $ ; ELSE window,4,xs=512,ys=256, xpos = 680, ypos = 770 END ELSE BEGIN nbc = fix(sqrt(status.nbimages))+1 window,4,xs=dim*nbc,ys=dim*nbc END END print,"### conditions initiales ###" ; print, 'erreur_offset', err_offset err = erreur_grad(param) + err_offset print,"erreur: ",err ;pas panique pour err_offset ca ne joue pas dans la minimisation ;c'etait juste pour l'histoire du choix de l'hyper-parametre de regul errbis = err IF keyword_set(visuim) THEN BEGIN wset,3 IF keyword_set(LOG) THEN aff,alog10(objnew > max(objnew)/log),0,/sa ELSE $ aff,objnew,0,/sa aff,objnew-objold,1,/sa END IF keyword_set(visupsf) THEN BEGIN IF (status.nbimages EQ 1) THEN BEGIN wset,4 aff,psfnew,0,/sa aff,psfnew-psfold,1,/sa END ELSE BEGIN wset,4 for ii=0,status.nbimages-1 do tvscl,psfnew(*,*,ii),ii END END paraminit = param info = dblarr(4) ;deconvolution i = 0L IF keyword_set(OBJREF) THEN error = fltarr(itmax) ;plot des erreurs ;initialisation des garde-fous continue = 1B ; ws = 1B ; gradient = 1B IF proj EQ 1 THEN active_set = param*0. + 1B ; = active_mask = pixels sur lesquels on ; minimise ;stop IF (itmax NE 0) THEN BEGIN REPEAT BEGIN IF keyword_set(VMLM) THEN BEGIN errx = 0. ; pas de errx dans FMIN_OP. fmin_op, param, err, errf = errf, ws = ws, GRADIENT $ = gradient, $ FUNC = "erreur_grad", INIT = (i EQ 0), $ ACTIVE_SET = active_set, CONV_THRESHOLD= conv ENDIF ELSE BEGIN fmin_gpa, param, err, errx = errx, errf = errf,$ FUNC = "erreur_grad", INIT = $ (i EQ 0), REINIT = (i MOD periode EQ 0), $ ACTIVE_SET = active_set ENDELSE ;fmin_conjgrad, param, err, conv, FUNC = ;"erreur_grad", /INITIALIZE ,/USE_DERIV ;DEBUG ;noramilasation de param pour ;eviter les cstes multiplicatives qui ;n'influencerais pas la psf mais son gradient IF (status.normpsf EQ 1B) THEN $ param(dim2 :*) = param(dim2 :*)/total(param(dim2:*))*total(paraminit(dim2 :*)) err1 = err+err_offset err = erreur_grad(param,ERR2=err2) IF ((NOT keyword_set(quiet)) OR (i EQ (itmax-1))) THEN BEGIN print,"" print,"### iteration",i," ###" print,"erreur : ", strcompress(string(err), $ /remove_all) ; print,"err1: ",err - err2," err2:",err2 ; print, "" print,"convergence erreur : ", strcompress(string(errf), $ /remove_all) print,"convergence parametres : ", strcompress(string(errx), /remove_all) print,"" ENDIF objnew = param2obj(param,DIM=status.dim,POS=status.pos) renorm = kkk IF keyword_set(visuim) THEN BEGIN wset,3 IF keyword_set(LOG) THEN aff,alog10(objnew > $ max(objnew)/log),0,/sa ELSE $ aff,objnew,0,/sa diff = objnew-objold diff(0:status.dim/10.,0:status.dim/10.) = 0. aff,diff,1,/sa xyouts,256+5,5,'0',/dev END ;test de convergence et de valeur de l'erreur errconv = convergence(objold,objnew) if not (status.psfcste) then begin ;JMC - DEDUG FOR jj=0, status.nbimages-1 DO psfnew(*,*,jj) = param2psf_mistral(param,jj) IF keyword_set(visupsf) THEN BEGIN IF (status.nbimages EQ 1) THEN BEGIN wset,4 aff,psfnew,0,/sa aff,psfnew-psfold,1,/sa END ELSE BEGIN wset,4 for ii=0,status.nbimages-1 do tvscl,psfnew(*,*,ii),ii END END END errbis = err1 ; IF keyword_set(save1it) THEN BEGIN ; IF (save1it EQ i) THEN BEGIN ; IF keyword_set(directory) THEN cd, directory ; IF keyword_set(nomfich) THEN BEGIN ; writefits,strcompress(nomfich+string(i)+'it.fits',/remove_all), param ; END ELSE BEGIN ; writefits,strcompress('param_'+string(i)+'it.fits',/remove_all), param ; ENDELSE ; print, " Sauvegarde du vecteur paramètre à la", i, " itération" ; ENDIF ; ENDIF ; IF keyword_set(saveit) THEN BEGIN ; IF ((i MOD saveit) EQ 0) THEN BEGIN ; IF keyword_set(directory) THEN cd, directory ; IF keyword_set(nomfich) THEN BEGIN ; writefits,strcompress(nomfich+string(i)+'it.fits',/remove_all), param ; END ELSE BEGIN ; writefits,strcompress('param_'+string(i)+'it.fits',/remove_all), param ; END ; print, " Sauvegarde du vecteur paramètre à la", i, " itération" ; END ; ENDIF info(0) = err info(1) = status.hyper_obj info(2) = status.hyper_psf info(3) = kkk IF keyword_set(OBJREF) THEN BEGIN objest = param2obj(param, pos = status.pos $ , dim= status.dim,norm= info(3)) error[i] = eqm(objest(where(objref NE 0)), objref(where(objref NE 0))) IF keyword_set(VERROR) THEN BEGIN IF i EQ 0 THEN BEGIN window, 0, xsize = 500, ysize = 400 ENDIF ELSE BEGIN wset, 0 IF i EQ 1 THEN plot_oi, findgen(i+1) $ , error(0:i), xr = [1, itmax], yr = [1, error(1)]ELSE $ IF i GT 1 THEN oplot, findgen(i+1), error(0:i) ENDELSE ENDIF ENDIF IF (getenv('DISPLAY') NE '') THEN BEGIN;get_kbrd doit avoir un tty en input IF (get_kbrd(0) EQ 'Q') THEN BEGIN; interruption propre des itérations continue = 0B message, 'interruption des itérations.', /INFO critere_keybord = 1B ENDIF ENDIF i = i+1L END UNTIL ((errx LT conv) AND (errf LT conv)) OR (i GE itmax) OR $ (NOT continue) END END