; $Id: mistral.pro,v 1.13 2006/06/01 13:35:45 mugnier Exp $ FUNCTION mistral, image, psf, OBJGUESS = objguess, VARNOISE = varnoise $ , REGOBJ = regobj, THRESHOLD = threshold $ , PSDOBJ = psdobj, AVGOBJ = avgobj, MYOPIC= myopic $ , PSDPSF = psdpsf $ , REGPSF = regpsf, AVGPSF = avgpsf, PSFES = psfes $ , NBIT = nbit , CONV = conv , PROJ = proj, NOPOS = nopos $ , VISU = visu, LOG = LOG, VERSION = version, HELP = help $ $ , COPYRIGHT = copyright, QUIET = quiet, VMLM = vmlm, SFIELD $ = sfield, info = info ;+ ;NAME : ; MISTRAL - edge preserving classical/myopic deconvolution ; ;CATEGORY : ; Signal Processing ; ;SYNTAX : ; objes = mistral(image, psf, [OBJGUESS, VARNOISE, REGOBJ, THRESHOLD, ; MYOPIC, PSDPSF, REGPSF, AVGPSF, PSFES, NBIT, CONV, /NOPOS,/LOG,/CRITERE_KEYBORD, ; /VISU, /VERSION , /HELP, /COPYRIGHT]) ; ;DESCRIPTION : ; ; mistral.pro is a shell that sets default values for various parameters ; before calling the deconvolution core "decomyro_gpa". Stopping rules : ; - conv criterion lower than conv (input parameters, default 1.e-5) ; - or iteration number equal itmax (input parameters, default 100000) ; - or Keybord value = Q ; ; ; MANDATORY ARGUMENTS : ; ; image : (input) image {dim x dim} to be deconvolved [photo-e-/pixel] ; flat-fielded and background-subtracted ; ; psf : (input) point-spread-function [PSF] {dim x dim} ; (starting point if myopic) ; ; OPTIONAL KEYWORDS : ; ; OBJGUESS : (input) starting point {dim x dim} for object estimation ; (default = image) ; ; VARNOISE : (input) noise variance {dim x dim} [e-^2/pixel] ; (default deduced from image) ; ; REGOBJ : (input) global object regularization factor {scalar} ; (default = 1.) ; ; THRESHOLD: (input) tuning {scalar} for threshold in object regularization ; (default = 1.) ; PSDOBJ : (input) PSD {dim x dim} for object regularization for a ; Wiener-like regularization ( ; (default PSD=0 => L2-L1 regularization) ; PSD obj must be given in photon^2 ; => max(psd) ~ (total(image))^2 ; AVGPOBJ : (input) mean obj {dim x dim} for DSPobj regularization ; (default = 0) ; (not used in L2-L1 regularization) ; ; /MYOPIC : (input) flag for myopic dec. (default = classical) ; ; PSDPSF : (input) PSD {dim x dim} for PSF regularization ; (default deduced from circular average of OTF^2) ; ; REGPSF : (input) global PSF regularization factor {scalar} ; (default = 1.) ; ; AVGPSF : (input) mean PSF {dim x dim} for regularization ; (default = psf starting point) ; ; PSFES : (output) estimated PSF {dim x dim} ; ; NBIT : (input) maximum number of iterations even if convergence criterion ; not satisfied (default = 100000) ; ; CONV : (input) stopping rule: convergence in the evolution of both the ; criterion and the parameters (default = 1.e-5) ; ; LOG : (input) log-scale visualisation. log gives the lower limit in ; comparison of the maximum of the object : (aff,alog10(obj > max(obj)/log)) ; ; QUIET : (input) desactivates prints at each iterations ; ; /NOPOS : (input) desactivates positivity constraint ; (default = positivity on object and PSF if myopic) ; ; /PROJ : (input) activates positivity using a projected gradient method ; (using Kuhn & Tucker condition) ; ; /VMLM : (input) use "OptimPack" (minimization package developped by E. ; Thiebaut in C) based on a VMLM approach instead of a modified ; conjugate-gradient method. The VMLM method is faster than ; conj.-grad. Warning, this needs to have the "OptimPack" ; package which is not included in Mistral package. ; ; /VISU : (input) activates displays ; (default = no display) ; ; /VERSION : (input) prints version number before running ; ; /HELP : (input) prints syntax and exits ; ; /COPYRIGHT: (input) prints information about copyright and exits ; ;ERROR DIAGNOSTIC : ; ; ;SEE ALSO : ; decomyro_gpa ; ;AUTHOR : ; $Author: mugnier $ ; ;HISTORY : ; $Log: mistral.pro,v $ ; Revision 1.13 2006/06/01 13:35:45 mugnier ; Nouveaux mots-cl�s /VMLM (minimisation by OptimPack) et /SFIELD (white prior, ; i.e., independent pixels, for star fields). ; ; Revision 1.12 2003/04/16 12:02:23 fusco ; New keyword : proj => positivity constraint using projection method. ; ; Revision 1.11 2003-02-18 19:42:28+01 fusco ; Changements mineurs : ; - possibilit� d'arret "propre" en tapant Q (shif q). Dans ce ; cas mistral s'arrete proprement a l'iteration suivante ; - mot cle log : possibilite de visu logaritmique : la dynamique est fixe par : ; max(image) et max(image)/log. ; ; Revision 1.10 2003-02-18 19:17:04+01 fusco ; *** empty log message *** ; ; Revision 1.9 2000-10-16 11:50:51+02 fusco ; correction d'un facteur de normailasation pour que mistral soit completement ; compatible avec wiener_est.pro et dsp_fit. ; ; Revision 1.8 2000-10-13 08:35:53+02 fusco ; *** empty log message *** ; ; Revision 1.7 2000-10-13 08:32:55+02 fusco ; *** empty log message *** ; ; Revision 1.6 2000-10-13 08:27:46+02 fusco ; changement de nom pour param2psf -> param2psf_mistral ; ; Revision 1.5 2000-10-13 08:22:56+02 fusco ; mise a jour de MISTRAL a partir de version de travail ThF : rajout des mots cle PSDOBJ (pour regul type Wiener) et AVGOBJ (idem) ; ; Revision 1.4 2000-02-16 16:38:41+01 conan ; ajout d'un call_procedure devant 'doc_library' (interet?) ; c'est la premiere version distribuee a Tatulli Roddier et Storrs ; ; Revision 1.3 2000-02-16 15:14:14+01 conan ; POS keywords replaced by NOPOS, COPYRIGHT keyword is added. ; ; Revision 1.2 2000-02-15 17:33:49+01 conan ; ajout de la doc, de copyright, THRESHOLD=1. par defaut (4 incorpore dans ; formule ad-hoc) ; ;- on_error,2 IF keyword_set(version) THEN $ message, '$Revision: 1.13 $, $Date: 2006/06/01 13:35:45 $', /INFO IF keyword_set(copyright) THEN BEGIN call_procedure, 'doc_library', 'usage_mistral' message, 'if you accept these conditions you can now run mistral.pro' ENDIF message,'* COPYRIGHT (C) Conan Mugnier Fusco - ONERA 1998-2000.', /info message,'This software, i.e., MISTRAL and the routines called therein ' + $ '(notably decomyro_gpa), are Copyright ONERA.', /info message,'The use of MISTRAL implies acceptance of the conditions described ' + $ 'in usage_mistral.txt', /info IF (n_params() NE 2) OR keyword_set(help) THEN BEGIN message, 'Syntax:', /INFO call_procedure, 'doc_library', 'mistral' message, 'Call MISTRAL with correct syntax.' ENDIF IF keyword_set(VMLM) THEN BEGIN print, '------------------------------------------------------------------------' print, "Use of OptimPack (minimization package developped by E. Thiebaut)" print, "based on a VMLM approach instead of a modified conjugate-gradient method" print, '------------------------------------------------------------------------' ENDIF ;Default parameters for decomyro : dimim = float((size(image))[1]) IF (n_elements(OBJGUESS) EQ 0L) THEN objguess = image > 0. IF (n_elements(VARNOISE) EQ 0L) THEN BEGIN minim = min(image) IF minim GE 0 THEN begin printf, -2 , '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' printf, -2 , '!!!!!! Warning: image strictly positive !!!!!!!!!!!!!!!!!!!!!' printf, -2 , '!!!!! Is the background correctly subtracted?? !!!!!' printf, -2 , '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' IF minim EQ 0 THEN BEGIN casdivzero = 1000. varnoise = image varnoise(where(varnoise LT max(varnoise)/casdivzero))= $ max(varnoise)/casdivzero ENDIF ELSE varnoise = image ENDIF ELSE $ ; varnoise = photon_noise + estimation of detector/background noise variance: varnoise = (image > 0 + (abs(mean(image(where(image lt $ 0))))*sqrt(!pi/2.))^2) ;with centered Gaussian statistics of standard deviation sigma, ;the mean on negative points gives: sigma/sqrt(pi/2) ENDIF ;"regobj" corresponds to \mu in Eq.3 [Conan-98]: IF (n_elements(PSDOBJ) EQ 0L) OR (n_elements(PSDOBJ) NE n_elements(image)) $ THEN BEGIN psdobj = 1B avgobj = 0B ENDIF ELSE BEGIN ;case of a Wiener-like regularization psdobj = psdobj/total(image)^2 ; PSD must be normalized (in term of photon ; number ENDELSE IF (n_elements(REGOBJ) EQ 0L) THEN BEGIN regobj = 1. ENDIF ELSE BEGIN ;cas d'une regularisation non uniforme IF n_elements(regobj) EQ dimim^2 THEN regnounif = 1B IF min(regobj) LE 0 THEN message, "WARNING: REGOBJ MUST BE POSITIVE" ENDELSE IF (n_elements(THRESHOLD) EQ 0L) THEN threshold = 1. ;input parameter THRESHOLD is only a factor applied to the following ad-hoc ;expression: seuil = threshold*mean([sqrt(variance(image-shift(image, 0, 1))), $ sqrt(variance(image-shift(image, 1, 0)))]) ;"seuil" is really the "threshold" \delta of Eq.3 [Conan-98] pos = 1 IF (n_elements(NBIT) EQ 0L) THEN nbit = 100000. IF (n_elements(CONV) EQ 0L) THEN conv = 1.e-5 IF keyword_set(NOPOS) THEN pos = 0B IF keyword_set(PROJ) THEN pos = 2 IF (n_elements(VISU) EQ 0L) THEN visu = 0B IF (n_elements(REGPSF) EQ 0L) THEN regpsf = 1 ELSE BEGIN IF (n_elements(REGPSF) GT 1L) THEN message, "WARNING : REGPSF MUST BE A " + $ "SCALAR" ENDELSE IF keyword_set(LOG) THEN BEGIN IF log LE 0 THEN message, 'Log must be a positive number' ENDIF IF (n_elements(MYOPIC) EQ 0L) THEN myopic = 0B ;parameter display: printf, -2, "" printf, -2, "INPUT PARAMETERS" printf, -2, "------------------------------" printf, -2, "" IF myopic NE 1B THEN printf, -2, "CLASSICAL DECONVOLUTION" ELSE BEGIN $ printf, -2, "MYOPIC DECONVOLUTION" printf, -2, "PSF reg: hyper = ", strcompress(string(regpsf), /remove_all) ENDELSE printf, -2, "% Positivite :", pos IF n_elements(psdobj) gt 1 THEN printf, -2, $ "% Object reg : Quadratic with PSD obj" ELSE BEGIN printf, -2, "% Object reg : Linear-Quadratic" printf, -2, "% Object reg : threshold = ", $ strcompress(string(threshold), /remove_all) ENDELSE printf, -2, "% Object reg : hyper = " , strcompress(string(regobj), /remove_all) printf, -2, "------------------------------" printf, -2, "" if keyword_set(info) then inf = 1B IF myopic NE 1B THEN BEGIN ;case of the classical deconvolution periode = round(float(dimim)^2/12.)+3 ;tuning parameter for minimization ;by Fmingpa decomyro_gpa, image,1./varnoise, param, OBJGUESS = objguess $ , PSFGUESS = psf, TERMREG_OBJ = psdobj, MOYOBJ = avgobj $ , IT = nbit, INFO = inf, hyper_obj=regobj $ , VISUIM = visu, LOG = log, CONV = conv, $ PERIODE = periode, POS = pos $ , /psfcste, seuil = seuil, VMLM = vmlm, $ $ CRITERE_KEYBORD = critere_keybord, QUIET = quiet, SFIELD = sfield,ERROR = valcritere objes = param2obj(param, pos = pos, dim= dimim,norm= inf(3)) info = inf ENDIF ELSE BEGIN ;case of the myopic deconvolution ;default value for PSDPSF [based on circular mean of OTF^2]: IF (n_elements(PSDPSF) EQ 0L) THEN BEGIN ftom2 = circmoy(abs(eclat(fft(psf, -1)))^2, /milieu) ftom2 = ftom2/max(ftom2) psdpsf = fltarr(dimim, dimim) selectfreq = eclat(distc(dimim)) FOR iii = 0, n_elements(ftom2)-2 DO BEGIN selectfreq1 = (selectfreq lt (iii+1)) and (selectfreq ge (iii)) psdpsf(where(selectfreq1 eq 1)) = ftom2(iii) ENDFOR psdpsf(where((psdpsf) LT 1.e-6))= 1.e-6 ;avoids division by 0 in PSF reg. ENDIF IF (n_elements(AVGPSF) EQ 0L) THEN avgpsf = psf ELSE BEGIN IF n_elements(avgpsf) ne float(dimim)^2 THEN message, "warning : " + $ "AVGPSF and IMAGE must have the same dimensions" ENDELSE periode = round(2.*float(dimim)^2/12.)+3 ;tuning parameter for ;minimization by Fmingpa decomyro_gpa, image,1./varnoise, param, OBJGUESS = objguess $ , PSFGUESS = psf, TERMREG_OBJ = psdobj, MOYOBJ = avgobj $ , IT = nbit, INFO = inf, hyper_obj=regobj $ , VISUIM = visu, LOG = log, CONV = conv, PERIODE = periode, POS = pos $ , /regpsf, seuil = seuil, AVGPSF = avgpsf, TERMREG_PSF = psdpsf, $ hyper_psf = regpsf, /normpsf , VMLM = vmlm, SFIELD = sfield, $ CRITERE_KEYBORD = critere_keybord, QUIET = quiet, VISUPSF= visu, error = error objes = param2obj(param, pos = pos, dim= dimim,norm= inf(3)) psfes = param2psf_mistral(param, 0) info = inf ENDELSE return, objes END