PRO FILTR_ECH, nx, ie, iemax, a_in, a_ex, $ ; entree a_red_fil, $ ; entree ax_e0, filtre ; sortie ; But: produire un filtre annulaire dans l'espace de Fourier, de facon a selec- ; tionner un domaie d'echelles spatiales. ; Le filtre est un anneau a bord adoucis, avec un rapport environ 2 entre ; le bord externe et le bord interne. On prend un anneau elliptique dont ; les axes EW et NS sont definis plus bas. Le bord externe est 2 fois ; moins raide que le bord interne, de facon que la somme de filtres ; consecutifs soit plate. ; Pour ie = iemax (petites echelles, la limite exterieurs du filtre est pri- ; se, par simplicite), comme l'ellipse englobant les barres de E2 avec les ; NS et de NS45 avec les EW, et la limite exterieure est abrupte. ; Les limites interieures sont progressive, sauf pour ie = 0 ou le domaine ; est un disque englobant l'origine. ; Rem : le filtre en sortie est centr'e au coin inferieur gauche. ; Modifications : ; 01 jan 10 : choix possible de iemax dans WCLEAN. ; jan 11 : suivant la remarque de Kerdraon on prends des zones rigoureuse- ; ment homothetiques et un peu plus petites. En fait qui semble ; le plus important est que la somme des filtres fasse 1. ; 12 : ie = 0 correspond aux plus grandes echelles et iemax aux ; petites. ; 22 : passage en argument d'un facteur de reduction des zones de fil- ; trage. ; 23 : transition externe progressive (et non plus brutale) pour l'an- ; neau de filtrage le plus externe (echelles les plus petites). ; 02 avr 17 : ajout de ax_e0 a la liste de sortie pour le lissage final dans ; RH_CLEAN. ; 02 mai 24 : la couronne externe doit etre contrainte par le bord du pave. Les ; essais du 24 mai 02 montrent que les resultats sur la TF de ; l'image propre non filtree sont ainsi bien meilleurs, surtout ; avec une source compacte. L'effet de a_red_fil=0.7 etait simple ; ment de reduire assez la couronne externe pour la contraindre. ; 06 fev 23: amlioration generale et simplification: ; - les limites des zones de filtrage sont exactement circulaires, ; ce qui m'a evit'e de me fatiguer pour m'assurer que la somme de ; filtres elliptiques est exactement egale a 1. ; - les limites des zones de filtrages successives sont exactement ; dans un facteur 2 (ce qui n'etait pas exactement le cas, alors ; que les largeurs des transitions l'etaient). Il en resulte que ; la somme des filtres est maintenant exactement egale a l'unite. ; - les transitions sont adoucies (non plus lineaires, mais en 1/2 ; periode de cosinus) et elargies au maximum, de facon que les ; filtres n'aient plus de plateau strict (il est indiqu'e dans ; rh_dpnew que a_in et a_ex doivent etre 2/3 et 4/3). ; - il est conseill'e dans rh_dpnew de figer aussi larg_2 = [18, 25] ; => on arrive a des artefacts ~ 1-2%, dus a la couronne externe, et ; sans doutes dus a la nature de la couverture uv au dela du ; pave central. ; Notations: ; a_in a_ex definissent les transitions progressives internes et externes ; des zones annulaires de filtrage. ; Choisis dans RH_DPNEW (usuel 0.85 et 1.15). ; On appelle limite interne (resp externe) d'une transition la ; plus proche (resp la plus loin) de l'origine. Cela ne fait ; pas reference a l'interieur ou l'exterieur de l'anneau du ; filtre. ; ie numero du filtrage (0 pour le filtrage le plus interne). ; iemax numero du filtrage le plus externe (0 pour n'avoir aucun ; filtrage d'echelle). ; nx nbre de pts sur une dimension de l'image (en entree). if (iemax eq 0) then begin ; pas de filtrage d'echelle) filtre = replicate(1, nx, nx) goto, a1 endif ; Production de matrices nx*nx mx et my contenant respectivement les ; abscisses et les ordonnees des points du plan (u,v) pr au centre ; (nx/2, nx/2). vx = findgen(nx) - nx/2 vy = findgen(nx) - nx/2 ; vx et vy sont les vecteurs des colonnes et lignes mx = vx # replicate (1, n_elements(vy)) ; nl lignes identiques. my = replicate (1, n_elements(vx)) # vy ; nc colonnes identiques. ; Remplissage bidon des elements (64,64) pour eviter des ennuis avec atan. mx(nx/2, nx/2) = 0.001 my(nx/2, nx/2) = 0.001 ; Controle imp = 0 if (imp eq 1) then begin help, Mx, My shade_surf, congrid(mx, 512, 512, cubic=-0.5) window, 1 & wset, 1 shade_surf, congrid(my, 512, 512, cubic=-0.5) stop endif ; Conclusion: OK. ; Fin de controle ; Axes des iemax+1 ellipses delimitant les bandes passantes des filtres. axes = fltarr(iemax+1, 2) ; 1er indice de "axe" : numero de l'ellipse, 2eme indice : EW ou NS. ; idees guide : ; . environ un facteur 2 pour chaque zone ; . le NS n'est pas contraint au dela de 45, inutile de remplir au ; dela. ; . l'EW va a 64, donc on peut remplir un peu au dela de la barre de ; E2 a 48 ; . la couronne externe doit etre contrainte par le bord du pave. Les ; essais du 24 mai 02 montrent que les resultats sur la TF de ; l'image propre non filtree sont ainsi bien meilleurs, surtout ; avec une source compacte. L'effet de a_red_fil=0.7 etait simple ; ment de reduire assez la couronne externe pour la contraindre. if (iemax eq 4) then begin axes(0, *) = [ 3, 4] axes(1, *) = [ 6, 8] axes(2, *) = [10, 12] axes(3, *) = [15, 22] axes(4, *) = [60, 48] ; garder les barres de E2 et NS45 endif ; if (iemax eq 4) then begin ; cadavre avec la couronne externe non ; axes(0, *) = [ 3, 4] ; contrainte par le pave. ; axes(1, *) = [ 6, 8] ; axes(2, *) = [12, 14] ; axes(3, *) = [23, 27] ; axes(4, *) = [60, 48] ; garder les barres de E2 et NS45 ; endif ; Sauvegarde du 23 fev 2006 pour iemax = 3 ; if (iemax eq 3) then begin ; en general choisis dans RH_DPNEW. ; axes(0, *) = [ 4, 4] ; axes(1, *) = [ 8, 10] ; axes(2, *) = [15, 22] ; axes(3, *) = [60, 48] ; garder les barres de E2 et NS45 ; endif ; Rem du 23 fev 2006 : il faut que la somme des filtres soit egale ; a l'unite, et le rapport des pentes internes et externes doit ; etre exactement egal au rapport des limtes internes et externes ; pour chaque filtre, sauf pour la couronne externe qui peut etre ; libre. if (iemax eq 3) then begin ; d'apres essais du 23 fev 2006. axes(0, *) = [ 4, 4] axes(1, *) = [ 8, 8] axes(2, *) = [16, 16] axes(3, *) = [64, 64] endif if (iemax eq 2) then begin axes(0, *) = [10, 10] axes(1, *) = [15, 22] axes(2, *) = [52, 48] ; garder les barres de E2 et NS45 endif if (iemax eq 1) then begin axes(0, *) = [12, 12] axes(1, *) = [52, 48] ; garder les barres de E2 et NS45 endif ; Reduction des zones de filtrage axes = a_red_fil * axes ; Sauvegarde des zones definies avant le 10 jan 01 ; axes(0, *) = [ 3, 4] ; axes(1, *) = [ 5, 8] ; axes(2, *) = [11, 16] ; axes(3, *) = [30, 30] ; axes(4, *) = [64, 50] ; axe EW, axes NS ; Calcul des axes des 4 ellipses definissant pour chaque filtre les limites in- ; ternes et externes des transitions internes et externes ; Rappel : on appelle limite interne (resp externe) d'une transition ; celle qui est la plus proche (resp la plus loin) de l'origine. Cela ; ne fait pas reference a l'interieur ou l'exterieur de l'anneau du ; filtre. ; Adoucissement de la transition externe de la couronne pour iemax da = 1 - a_in ; 1/2 largeur d'une transition interne. On pourrait ; aussi bien la definir a partir de a_ex. fac_adou = 1.0 ; facteur d'adoucissement , < 2.5 pour eviter secon- ; daires visibles sur la croix des secondaires. ; fac_adou < 1 durcit et fac_adou > 1 adoucit. ; Apres le 23 fev 2006, les filtres sont constitues de ; deux transitions en 1/ periodes de cosinus, sans ; plateau constant. Il faut donc prendre fac_adou = 1 ; pour que la somme des filtres fasse 1. ; b_in = 1 - fac_adou * da ; avant 23 fev 2006 ; b_ex = 1 + fac_adou * da ; avant 23 fev 2006 ;; b_in = 1 - fac_adou * da ; apres 23 fev 2006 b_in = 1 - 2.0 * da ; apres 23 fev 2006 car la couronne ; externe est tres large b_ex = 1 ; apres 23 fev 2006 ; on fait en sorte que le filtre tombe a zero a la limite du plan uv ; explore. La formulation d'avant le 26 fev 2006 faisait que le filtre ; etait encore a 0.5 a cette limite, et a cette limite le filtre tom- ; bait evidemment brutalement a zero => fortes oscillations sur image. ; Rappel: "axes" : 1er indice numero de l'ellipse, 2eme indice EW ou NS ; notations : suffixe "0" pour la transition du cote ou elle ; est nulle. ; En pratique c'est isotrope depuis le 23 fev 2006 (et iemax est ; fig'e a 3). ; Notation pour les transitions ci-dessous : ; Transition interne, limite interne if (ie eq 0) then ax_i0 = [0.0, 0.0] if (ie eq 1) then ax_i0 = [0.0, 0.0] if (ie gt 1) then ax_i0 = reform (a_in * axes (ie - 1, *)) ; Transition interne, limite externe if (ie eq 0) then ax_i1 = [0.0, 0.0] if (ie gt 0) then ax_i1 = reform (a_ex * axes (ie - 1, *)) ; Transition externe, limite interne if (ie eq 0) then ax_e1 = [0.0, 0.0] if ( (ie gt 0) and $ (ie lt iemax) ) then ax_e1 = reform (a_in * axes (ie , *)) if (ie eq iemax) then ax_e1 = reform (b_in * axes (ie , *)) ; Transition externe, limite externe if (ie lt iemax) then ax_e0 = reform(a_ex * axes (ie , *)) if (ie eq iemax) then begin ax_e0 = reform(b_ex * axes(ie , *)) ax_e0_2 = [55, 52] ax_e0 = ax_e0 > ax_e0_2 ; Pour englober de toutes facons les barres de E2 et NS45 (modif du ; 08 mar 02) ; La derniere ligne peut se mettre en commentaire pour ne pas en- ; glober les barres E2 et NS45. endif ; Controle des ellipses d'un filtre icon = 0 ; if ((icon eq 1) and (ie eq iemax)) then begin if (icon eq 1) then begin print, ' ' ch_format ="('FILTR_ECH : axes des ellipses pour filtrage', i2," + $ "' EW NS')" print, format=ch_format, ie print, ' limite interne de la transition interne:', ax_i0 print, ' limite externe ------------------------:', ax_i1 print, ' limite interne de la transition externe:', ax_e1 print, ' limite externe ------------------------:', ax_e0 stop endif ; Conclusion: OK. ; Fin de controle ; Intersections des directions portant les rayons vecteurs definis par mx et my ; avec les ellipses. ; Repres d'une ellipse par : x = a cos(phi), y = b sin(phi) (1) ; ou x et y sont les coord d'un point sur l'ellipse. ; L'angle polaire ap du pt (x, y) est telle que: ; tan(ap) = y/x = b/a tan(phi) => phi = atan(a/b * y/x) (2) ; L'eq (2) associe une valeur phi a tous les points situes sur la droi- ; te D joignant l'origine au point choisi sur l'ellipse. ; Les coord xI et yI de l'intersection de l'ellipse et de la droite D ; sont alors donnees par les relations (1). ; En IDL: ATAN(Y, X) = ATAN (Y/X) [atan(sin, cos) comme en fortran]. rv = sqrt (mx^2 + my^2) ; long du rayon vecteur. ; Parametres phi correspondant aux intersections avec les 4 ellipses. ; rappel : les ax_* contiennent 2 valeurs : l'axe EW et l'axe NS. phi_e0 = atan (ax_e0(0)*my, ax_e0(1)*mx) phi_e1 = atan (ax_e1(0)*my, ax_e1(1)*mx) if (ie gt 0) then begin ; transition interne seulement si ie>0. phi_i1 = atan (ax_i1(0)*my, ax_i1(1)*mx) phi_i0 = atan (ax_i0(0)*my, ax_i0(1)*mx) endif ; Longueurs des rayons vecteurs des points d'intersection avec les ellipses. r_e0 = sqrt ( (ax_e0(0)*cos(phi_e0))^2 + (ax_e0(1)*sin(phi_e0))^2 ) r_e1 = sqrt ( (ax_e1(0)*cos(phi_e1))^2 + (ax_e1(1)*sin(phi_e1))^2 ) if (ie gt 0) then begin ; pas de transition interne pour ie=0. r_i1 = sqrt ( (ax_i1(0)*cos(phi_i1))^2 + (ax_i1(1)*sin(phi_i1))^2 ) r_i0 = sqrt ( (ax_i0(0)*cos(phi_i0))^2 + (ax_i0(1)*sin(phi_i0))^2 ) endif ; Controle r_e0, r_e1, r_i1, r_i0 icon = 0 ; du 8 mars 02. if (icon eq 1) then begin print, 'Filtrage : ', ie if (ie gt 0) then begin print, ' min(r_i0), max(r_i1) =', min(r_i0), max(r_i1) print, ' min(r_i1), max(r_i1) =', min(r_i1), max(r_i1) endif print, ' min(r_e1), max(r_e1) =', min(r_e1), max(r_e1) print, ' min(r_e0), max(r_e0) =', min(r_e0), max(r_e0) stop endif ; Conclusion : OK. icon = 0 ; vereux if (icon eq 1) then begin print, 'numero de filtrage d''echelle =', ie window, 0 & wset, 0 ; print, 'r_e0 =', r_e0 ; print, 'r_e1 =', r_e1 ; print, 'r_i1 =', r_i1 ; print, 'r_i0 =', r_i0 ; stop shade_surf, congrid(r_e0, 512, 512, cubic=-0.5) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'longueur rayon vecteur transition externe, limite externe' window, 1 & wset, 1 tvscl , congrid(r_e0, 512, 512, cubic=-0.5) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'long ray vec transition externe, limite externe' stop window, 0 & wset, 0 shade_surf, congrid(r_e1, 512, 512, cubic=-0.5) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'long ray vec transition externe, limite interne' window, 1 & wset, 1 tvscl , congrid(r_e1, 512, 512, cubic=-0.5) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'long ray vec transition externe, limite interne' stop window, 0 & wset, 0 shade_surf, congrid(r_i1, 512, 512, cubic=-0.5) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'long ray vec transition interne, limite externe' window, 1 & wset, 1 tvscl , congrid(r_i1, 512, 512, cubic=-0.5) stop window, 0 & wset, 0 shade_surf, congrid(r_i0, 512, 512, cubic=-0.5) window, 1 & wset, 1 tvscl , congrid(r_i0, 512, 512, cubic=-0.5) stop endif ; Conclusion: OK. ; Fin de controle ; Definition du filtre filtre = replicate(1., nx, nx) ; 1 partout. filtre_1 = filtre ; ne sera pas modifie. ; Zone interne (seulement pour ie > 0) ; if (ie gt 0) then begin ; -> 2 juin 2006 if (ie gt 1) then begin filtre (where(rv le r_i0)) = 0.0 endif ; Zone de transition interne (seulement pour ie > 0) ; if (ie gt 0) then begin ; -> 2 juin 2006 if (ie ge 1) then begin frac = (rv - r_i0) / (r_i1 - r_i0) zone = where( (frac gt 0) and (frac lt 1) ) ; filtre (zone) = frac (zone) ; -> 23 fev 2006. ; filtre (zone) = frac ne marche pas car fitre(zone) et frac n'ont ; pas les memes dimensions (fitre(zone) a une seule dim). ; frac croit de 0 a 1 le long du rayon vecteur, et on veut engender ; une arche de sinus qui varie entre ses extrema de 0 a 1. ; On prend dans la zone de transition interne : ; filtre = 1/2 ( 1 - cos(pi * frac)) filtre (zone) = 0.5 * (1 - cos(!pi * frac(zone))) endif ; Zone de transition externe frac = (rv - r_e1) / (r_e0 - r_e1) zone = where( (frac gt 0) and (frac lt 1) ) ; filtre(zone) = filtre_1 - frac(zone) ; -> 23 fev 2006. ; rappel : canular avec "filtre" ; frac croit de 0 a 1 le long du rayon vecteur, et on veut engender ; une arche de sinus qui varie entre ses extrema de 1 a 0. ; On pose : dr = rv - r_i0 et Dr = r_i1 - r_i0 ; On prend dans la zone de transition externe : ; filtre = 1/2 (cos(pi * frac)) ; filtre(zone) = filtre_1 - frac(zone) ; avant fev 2006. filtre(zone) = 0.5 * (1 + cos( !pi * frac(zone) ) ) ; apres 23 fev 2006 ; Zone externe filtre (where(rv ge r_e0)) = 0.0 a1: J_M = 0 ; fin du saut dans le cas iemax=0 (pas de filtrage d'echelle). ; Controle final du filtre icon = 0 ; if ((icon eq 1) and (ie eq iemax)) then begin if (icon eq 1) then begin print, 'FILTR_ECH : filtrage', ie window, 0 shade_surf, congrid(filtre, 512, 512, cubic=-0.5) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'filtre centre en milieu de tableau' window, 1 tvscl , congrid(filtre, 512, 512, cubic=-0.5) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'filtre centre en milieu de tableau' window, 2 umax = fix(1.1 * (ax_e0(0) > ax_e0(1))) ; if (umax gt nx-1) then umax = nx - 1 ; -> 2 juin 2006. if (umax gt nx/2 - 1) then umax = nx/2 - 1 absc = -umax + indgen(2 * umax + 1) plot, absc, filtre(nx/2-umax:nx/2+umax, nx/2), yrange=[0, 1.1], $ xstyle=1 oplot, absc, filtre(nx/2, nx/2-umax:nx/2+umax) / 2, linestyle=1 ; /2 pour separer. xyouts, 0.5, 0.97, /normal, alignment=05, $ 'Coupes du filtre en u (___) et v(- - -) stop endif ; Conclusion: OK. ; Fin de controle ; Centrage dans les coins (cas de filtrage d'echelle effectif) if (iemax gt 0) then filtre = shift(filtre, -nx/2, -nx/2) end ; fin fe FILTR_ECH ;______________________________________________________________________________ ; appel de CLEAN_ECHELLE dans RH_CLEAN : ; ( verifie le 13 jun 02) ; CLEAN_ECHELLE, $ ; ; entree: ; mode_clean, i_ech_min, i_ech_max, i_ech, $ ; im_ech, l_th_ech, i_fil, dom_c_ech, $ ; dom_a, $ ; itermax, jmax, crit_arret, gain, $ ; i_seuil(i_ech), fac_crois_res, seuil_amp, $ ; i_seuil_max_rms, fac_seuil, seuil_flux, $ ; i_test_arret, $ ; i_seuil_max_rms et fac_seuil pour affichage. ; ; sortie: ; im_ech_s, im_ech_r ; ; Rappel 1 : im_ech, l_th_ech et i_fil sont ici definis sur le ; ; nbre de pts reduit. ; ; Rappel 2 : im_ech_s est "l'image synthetique", cad la somme ; ; des composantes clean habillees avec i_fil, sans addi- ; ; tion d'image residuelle. PRO CLEAN_ECHELLE, $ ; entree mode_clean, i_ech_min, i_ech_max, i_ech, $ im, l_th, i_fil, dom_c, $ dom_a, $ itermax, jmax, crit_arret, gain, $ i_seuil, fac_crois_res, seuil_amp, $ i_seuil_max_rms, fac_seuil, seuil_flux, $ i_test_arret, $ ; i_seuil_max_rms et fac_seuil pour affichage ; sortie im_s, im_r ; Rappel : i_seuil = 1 pour toutes les echelles apres le 15 mai 02. ; But : Nettoyage d'une image, filtree d'echelle lors de l'appel par RH_CLEAN. ; Contraintes a satisfaire : ; - le nettoyage doit conserver le flux : ; . en l'absence de filtrage d'echelle la fonction de remplacement doit avoir ; le meme flux que la fonction de soustraction. Si la fonction de sous- ; traction est l_th (de flux 1) la fontion de remplacement doit etre un ; dirac unite (le lissage de l'image propre peut se faire en aval). ; . avec filtrage d'echelle, l'image nettoyee finale est la somme des images ; nettoyees filtrees d'echelle. Supposons la procedure suivante : pour ; chaque filtrage d'echelle on definit : ; . l_th_ech qui est le resultat du filtrage de l_th (de total 1) par le ; filtre d'echelle (d'amplitude 1 dans sa bande passante) ; . i_fil qui est la TF du filtre (ou "image" du filtre) ; On adopte l_th_ech comme fonction de soustraction et i_fil comme fonc- ; tion de remplacement. ; im_s = Sigma(im_ech_s) (_s comme "synthetique") ; Flux(im_s) = Somme(im_s) = Somme(Sigma(im_ech_s)) ; Il est evident que le flux est conserv'e si l'image initiale im est ; constituee d'une ou quelques sources ponctuelles bien separees, puis- ; que la somme des flux des l_th_ech (flux tous nuls sauf un) vaut 1, de ; meme que la somme des flux des i_fil. Le resultat est vrai si les sour- ; ces sont reconnaissables et en nombre identique a toutes les echelles. ; Le 13 dec 00 j'ai tres envie que le resultat soit vrai en general pour ; une image quelconque, mais je ne l'ai pas demontre. ; - le nombre d'iterations ne doit pas etre trop eleve: il est plus rapide de ; retirer une composante compacte "a la louche" qu'a la "petite cuiller". ; Or max(l_th)<<1 =>> max(l_th_ech)<<1 et on ne retire a chaque opera- ; tion qu'une petite fraction du maximum de chaque image filtree d'echelle. ; On peut accelerer le processus et reduire le nombre des iterations en ; remplacant : ; . l_th_ech par l_th_ech / max(l_th_ech) ; . i_fil par i_fil / max(l_th_ech) ; Cela devrait reduire le nombre d'iterations d'un facteur 1/max(l_th_ech) ; Modifications : ; 00 dec 4 passage de dom_c a CLEAN_ECHELLE dans la liste d'arguments. ; 7 nouveau critere d'arret (en plus des autres) : debut d'appari- ; tions erratiques de composantes clean de signe oppos'e. ; 13 definition des fonctions de soustraction et de remplacement ; (voir commentaire plus haut) a partir de l_th_ech (obtenu ; par filtrage (gain 1 dans l'anneau passant) dans RH_CLEAN a ; partir de l_th calcule dans MALC_IM_2D avec un total unite. ; Rappel : le l_th connu dans CLEAN_ECHELLE est DEJA filtre ; d'echelle. ; 20 supression de toutes les delouiseries. niter disparait. ; 01 jan 25 modification du clean selectif destine a : ; . ameliorer la lisibilite des details faibles (ex arche de ; CME) en presence d'un centre intense, ; . nettoyer un centre intense de ses propres secondaires ; seulement pour gagner du temps ; Description detaillee dans dpnew. ; fev 14 tous les parametres de description du clean sont reportes dans ; dpnew. ; oct 4 ajout de fac_seuil dans la liste d'entree pour affichage. ; 02 mar 6 on etend le critere d'arret sur le changement de signe des ; composantes clean au cas de toutes les echelles. Cela donne ; de bons resultats pour supprimer les gros secondaires nega- ; tifs dus a une mauvaise calibration, si on n'ajoute pas ; l'image residuelle a l'image synthetique (mode_clean=2). ; annulation de la soustraction des composantes clean de l'image ; residuelle et de leur ajout a l'image synthetique pour les ; dernieres iterations (ayant provoque l'arret). ; 02 mai 24 la couronne externe doit etre contrainte par le bord du pave. Les ; essais du 24 mai 02 montrent que les resultats sur la TF de ; l'image propre non filtree sont ainsi bien meilleurs, surtout ; avec une source compacte. L'effet de a_red_fil=0.7 etait simple ; ment de reduire assez la couronne externe pour la contraindre. ; (modif faite dans FILTRE_ECHELLE). ; 02 mai 29 la moyenne des maxima des jmax images residuelles consecutives ; est faite sur les valeurs absolues et non les valeurs signees. ; Si seul le critere d'arret sur l'amplitude etait actif, cela ; pouvait arreter l'iteration trop tot car la moyenne pouvait ; etre << ses termes. En effet ca arrange beaucoup en simulation, ; il n'y a plus besoin de prendre coeff_seuil de l'ordre de 10^-6 ; pour depasser 50 iterations ! ; 02 oct 17 modifications pour diminuer le temps de calcul par iteration ; (8.6 msec pour la version jusqu'au 17 oct 02). Voir a "temps ; d'execution" plus bas. ; Notations: ; entree: ; im image a traiter, avec eventuellement un filtage d'echelle, ; i_ech indice du filtrage d'echelle. Ne sert qu'a un controle. ; i_fil TF du filtre d'echelle (image du filtre) centree en milieu de ; tableau, comme l_th. ; l_th lobe theorique utilise pour les soustactions CLEAN, centre en ; milieu de tableau, comme i_fil. ; seuil_amp valeur du max (en valeur absolue) de l'image residuelle pour ; laquelle on arrete l'iteration CLEAN. ; ; sortie: ; im_s image "synthetique". Somme des composantes clean habillees ; avec i_fil), sans addition de l'image residuelle. ; im_r image residuelle = (image en entree) - (composantes CLEAN) tab = size (im) nx = tab (1) ; dimension x de l'image. ny = tab (2) ; Controle: comparaison des lobes (filtres d'echelle) sale (soustrait) et pro- ; pre (ajout'e). icon = 0 if (icon eq 1) then begin print, 'CLEAN_ECHELLE : i_ech =', i_ech print, ' Total(l_th) =', total(l_th), $ ' Total(i_fil) =', total(i_fil) print, ' Max (l_th) =', max (l_th), $ ' max (i_fil) =', max(i_fil) print, ' ' ch_format = "(i2)" ch_id = string(format=ch_format, i_ech) window, 0 & wset, 0 shade_surf, congrid(l_th , 512, 512, cubic=-0.5) xyouts, 0.5, 0.97, /normal, alignment=0.5, 'l_th_ech i_ech=' + ch_id window, 1 & wset, 1 shade_surf, congrid(i_fil, 512, 512, cubic=-0.5) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'i_fil i_ech=' + ch_id window, 2 & wset, 2 shade_surf, congrid(im, 512, 512, cubic=-0.5) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'image filtree d echelle initiale i_ech=' + ch_id S_l_th = total(abs(l_th)) S_i_fil = total(abs(i_fil)) ; print, 'CLEAN_ECHELLE : i_ech =', i_ech ; help, nx, dom_c ; print, ' Totaux (abs(lobes sales et propre))', S_l_th, S_i_fil stop endif ; Conclusion: i_fil est seulement un peu plus propre que l_th pour i_ech=0 ; (ce qui est attendu du fait de la bonne couverture (u,v) ). La diffe- ; rence augmente evidemment avec i_ech croissant. ; Determination du maximum en valeur absolue de l'image avec son signe if((mode_clean eq 1) or (mode_clean eq 2) or (mode_clean eq 3)) then begin ; recherche des composantes sur toute l'image. max_im = max(abs(im), ip) endif if ( (mode_clean eq 4) or (mode_clean eq 5) ) then begin ; recherche des composantes dans dom_c seulement : ; . mode_clean = 4 on ne retire qu'elles, l'image finale est ; sale et residuelle. ; . ------------ 5 -- n'ajoute -------- l'image finale est ; partiellement propre. ; Rem : precaution pour que le calcul de la position du max soit ; le meme que dans le cas mode_clean=1 : ; max_im = max(abs(im(dom_c)), ip) ; donne un resultat totalement fantaisiste (vu le 8 fev 01). im_sous = im im_sous(dom_c) = 0 im_2 = im - im_sous ; nulle a l'ext de dom_c max_im = max(abs(im_2), ip) endif k_max = (ip mod nx) ; ind de colonne (orig coin inf gauche) l_max = (ip / nx) ; ------ ligne im_max = im(k_max, l_max) ; maximum avec son signe. ; Initialisation de l'mage synthetique (somme des composantes clean habillees ; avec i_fil) et de l'image residuelle im_s = 0. * float(im) im_r = im im_s_1 = 0. * float(im) ; idem du groupe d'iteration precedent. im_r_1 = im ; ------------------------------------ ; Initialisation des maxima sign'es de jmax images residuelles successives et ; de leurs positions et des valeurs quadratiques moyennes des images ; residuelles. max_r = fltarr(jmax) ; max d'image residuelle. max_x = fltarr(jmax) ; pour controles seulement. max_y = fltarr(jmax) ; ------------------------ quad_r = fltarr(jmax) ; val qud moy d'image residuelle sur dom_c. j_cum = 0L ; initialis compteur d'iterations ("cumul"). t = systime(1) ; Initialisation max image residuelle (avec signe) et drapeau d'iteration. if (i_seuil eq 1) then begin amp_preced = im_max ; en vue arret par croissance residu. endif if (i_seuil eq 2) then begin ; rappel : i_seuil = 1 pour toutes les echelles apres le 15 mai 02. if((mode_clean eq 1) or (mode_clean eq 2) or (mode_clean eq 3)) $ then mom = moment(im) if((mode_clean eq 4) or (mode_clean eq 5)) then mom = moment(im(dom_c)) amp_preced = sqrt(mom(1) + mom(0)^2) ; en vue arret par croissance residu. endif drap_iter = 1 ; passage oblige a l'entree dans CLEAN_ECHELLE. im_r est l'image resi- ; duelle, au debut identique a l'image. ; drap_iter peut etre mis a zero plus bas et arreter le while. ; Calcul flux de l'image initiale filtree d'echelle (sert seulement a un con- ; trole pour i_ech = 0) if (i_ech eq 0) then flux_initial = total(im) ; Fabrication de la fonction "voute" ; But: eviter de trouver le max au bord du champ quand on cherche le ; max d'une image obtenue sans aucun des harm de E0 (le champ in- ; terfero reel est moitie du champ adopte) le max d'une source ; au centre du champ (le Cygne par exemple) se trouve repet'e ; exactement sur les bords, et la fonction IDL "max" ne retient ; que le premier des maxima egaux trouves. ; Pour ecarter le cas on ajoute a l'image une voute basse culmi- ; nant a 1/100 du max d'image au centre du champ et diminuant a ; zero au bord du champ. ar = shift( dist(nx) / (nx/2), nx/2, nx/2) ; ray vect a partir du centre. voute = 1 - ar^2 ; voute unite ; Boucle sur la recherche des composantes CLEAN, par groupes de jmax : ; . calcul de l'image residuelle (par soustraction des composantes ; CLEAN de l'image initiale), ; . construction de l'image propre (composantes CLEAN seules) ; . tests d'arret t0 = systime(1) ; Definitions des fonctions de soustraction et de remplacement dans l'image ; residuelle (sorties de la boucle le 17 oct 02 pour gagner du temps) ; Ca gagne 1.4 msec par iteration (8.6 -> 7.2 msec) cmul = 1 ; pour essais seulement. f_sous = cmul * l_th / max(l_th) ; fonction de soustraction. f_remp = cmul * i_fil / max(l_th) ; fonction de remplacement. ; Rappel : total(l_th) = total(i_fil) = 1 while ((drap_iter eq 1) and (j_cum lt itermax)) do begin ; j_cum est le compteur d'iteration ("cumul"). ; Controle de l'image residuelle apres des goupes de jmax iterations icon = 0 ; if (icon eq 1) then begin if ((icon eq 1) and (i_ech eq 0)) then begin print, 'CLEAN_ECHELLE : i_ech =', i_ech shade_surf, im_r ch_nbre = string(format="(i4)", j_cum) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Image residuelle apres ' + ch_nbre + ' iterations' if (i_ech eq 0) then begin ; seul cas de flux nin nul. if (j_cum eq 0) then print, ' flux_initial =', flux_initial if (j_cum gt 0) then begin rap_flux = flux / flux_initial print, ' flux / flux_initial =', rap_flux endif endif stop endif ; Boucle de recherche de jmax composantes clean sans test d'arret. jc1 = j_cum ; pour essai for j_par = 0, jmax-1 do begin ; indice "j partiel". jc1 = jc1 + 1 ; pour essai ; Detection d'underflows icon = 0 if (icon eq 1) then begin if ( (i_ech eq 1) and (j_cum gt -1) ) then begin print, ' i_ech =', i_ech print, ' Debut boucle jmax : j_cum =', j_cum, $ ' j partiel =', j_par stop endif endif ; fin de recherche d'underflows. ; Controle de l'image residuelle apres chaque iteration icon = 0 ; if (icon eq 1) then begin if ((icon eq 1) and (i_ech eq 0)) then begin print, 'CLEAN_ECHELLE : i_ech =', i_ech ; on n'ouvre pas de fentre pour n'avoir pas a cliquer pour ; faire .c . nz = 8 ; zoomer de 64 - nz a 64 + nz. shade_surf, im_r(64-nz:64+nz, 60-nz:60+nz) ch_nbre = string(format="(i4)", j_cum + j_par) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Image residuelle apres ' + ch_nbre + ' iterations' if (i_ech eq 0) then begin ; seul cas de flux nin nul. if (j_cum eq 0) then print, ' flux_initial =', $ flux_initial if (j_cum gt 0) then begin rap_flux = flux / flux_initial print, ' flux / flux_initial =', rap_flux endif endif stop endif ; <- niveau de la boucle de recherche de jmax composantes sans ; tests d'arret ; Normalisation d'une voute faible ajoutee a l'image residuelle pour ; discriminer la position du centre d'orage de son alias (on ; suppose que le max de l'image est plus proche du centre du ; champ que son alias et on recherchera le max de la somme de ; l'image residuelle et d'un voute centree faible. ; Le but etant de distinguer entre un centre et son alias, on ; normalise la voute au max de l'image (en general un centre ; ou un sursaut), qu'il soit ou non dans dom_c. ; Rem : il peut arriver, notamment en haute frequence, que ; l'alias soit plus proche du centre du champ EW que la sour- ; ce. Alors la voute conduit a conserver l'alias (cas signale ; par K le 17 avr 02), probablement parce qu'elle contreba- ; lance l'effet "desaliasant" de E0. Pour cette raison on ; met un coeff de voute presque nul : 0.0001. ; Rem (17 oct 02) : compte tenu de la remarque precedente et du ; fait que l'effet desaliasant de E0 est surement > 0.0001, ; on supprime totalement le calcul de la voute et son ajout. ; toto = max(abs(im_r), ip) ; k = (ip mod nx) ; pour origine en bas a gauche. ; l = (ip / nx) ; ---------------------------- ; I MOD J is equal to the remainder when I is divided by J. ; coeff_voute = 0.01 * im_r(k, l) ; 0.01 * max d'image avec son ; signe. ; On ne cherche pas ici a dis- ; tinguer le centre de son ; alias. C'est suffisant pour ; normaliser la voute. ; Recherche du max en valeur absolue (reaffecte de son signe) sur le ; domaine dom_c de l'image residuelle et de sa position en ; vue de determiner une composante "sale" ; Rem : Pour mode_clean = 4 ou 5 il faut rechercher le max ; dans dom_c, mais il faut sa position dans le champ total. ; Or l'instruction max(image(dom_c), ip) ; fournit un indice ip qui ne traduit pas simplement la posi- ; tion dans le champ total. Cette erreur etait faite avant ; le 8 fev 01. La position 2D calculee etait fausse et les ; operations de soustraction et de remplacement fantaisistes. ; C'est pourquoi on annule im_r_2 en dehors de dom_c et on ; cherche le max dans le champ total. Le calcul de la posi- ; tion a 2D a partir de ip est alors le meme pour toutes les ; valeurs de mode_clean. ; ; <- niveau de la boucle de recherche de jmax composantes sans ; tests d'arret ; im_r_2 = im_r ; im_r_2 = im_r + coeff_voute * voute ; leger soulevement au centre. if ( (mode_clean eq 4) or (mode_clean eq 5) ) then begin ; on recherche le maximum dans dom_c . im_sous = im_r ; = im_r_2 -> 17 oct 02. im_sous(dom_c) = 0 im_r_3 = im_r - im_sous ; = im_r_2 -> 17 oct 02. icon = 0 if (icon eq 1) then begin window, 0 tvscl, congrid(im_r, 512, 512, cubic=-0.5) window, 1 tvscl, congrid(im_sous, 512, 512, cubic=-0.5) window, 2 tvscl, congrid(im_r_3, 512, 512, cubic=-0.5) stop endif im_r = im_r_3 ; im_r_2 = im_r_3 -> 17 oct 02. endif toto = max(abs(im_r), ip) ; im_r_2 -> 17 oct 02. k = (ip mod nx) ; pour origine en bas a gauche. l = (ip / nx) ; ---------------------------- x = k - nx/2 ; pour origine au centre. y = l - ny/2 ; ---------------------- ; I MOD J is equal to the remainder when I is divided by J. im_r_kl = im_r(k, l) ; max d'image avec son signe, ; (perdu -> 30 sept 99). indice = ip ; pour comparer avec dom_c. max_x(j_par) = x ; pour controles seulement. max_y(j_par) = y ; ------------------------ ; <- niveau de la boucle de recherche de jmax composantes sans ; tests d'arret ; Soustraction d'une composante de l'image brute ; Les composantes sont retirees de l'image sale sous la forme ; de l_th_ech, avec une couverture (u,v) pauvre au moins pour ; l'anneau de filtrage externe), et seront ajoutees a l'image ; propre sous la forme de i_fil avec une couverture (u,v) con- ; tinue. ; L'interet de l'operation est clair pour les anneaux de fil- ; trage externes. Il l'est moins pour les anneaux internes, si- ; tues dans le pave, ou le clean ne doit pas faire beaucoup ; plus qu'une simple apodisation. ; Rappel : i_fil et l_th sont centres en milieu de tableau. ; ; Suivant les valeurs de mode_clean : ; 1 2 ou 3 : soustraction des composantes sales quelles que ; soient leur position sur toute l'image. ; 4 ou 5 : soustraction des composantes sales interieures a ; dom_c ; Comme le max de l'image a ete determine plus haut sur toute ; l'image ou sur dom_c selon la valeur de mode_clean, il est ; inutile de distinguer ici entre les differents cas selon ; la valeur de mode_clean et on fait toujours la soustraction ; (ce qui ne sera pas le cas pour l'addition des composantes ; synthetiques. im_r = im_r - gain * im_r_kl * shift(f_sous, x, y) ; Detection d'underflows icon = 0 if (icon eq 1) then begin if ( (i_ech eq 1) and (j_cum gt -1) ) then begin print, ' i_ech =', i_ech print, ' plus loin boucle jmax : j_cum =', j_cum, $ ' j partiel =', j_par stop endif endif ; fin de recherche d'underflows. ; <- niveau de la boucle de recherche de jmax composantes sans ; tests d'arret ; Construction de l'image synthetique. ; Suivant les valeurs de mode_clean : ; 1 ou 2 (clean usuel) : addition de composantes synthetiques ; quelle que soit leur position. ; 3 : addition des composantes synthetiques si leur position ; est exterieure a dom_c. ; 4 : pas d'addition des composantes synthetiques (qui sont ; toutes interieures a dom_c). ; 5 : addition des composantes synthetiques (comme pour 1, mais ; elles sont interieures a dom_c). if((mode_clean eq 1) or (mode_clean eq 2)) then begin im_s = im_s + gain * im_r_kl * shift(f_remp, x, y) ; jusqu'au 30 sept 99 on soustrayait abs(im_r_kl) au lieu ; de im_r_kl), ce qui n'avait pas forcement le bon signe. endif if(mode_clean eq 3) then begin ; add comp synth ext a dom_c. toto = where( (dom_c - indice) eq 0, ic) ; si "indice" coincide avec un element de dom_c, ic=1. ; dans ce cas on n'ajoute pas la composante synthetique. if(ic eq 0) then begin im_s = im_s + gain * im_r_kl * shift(f_remp, x, y) endif endif if(mode_clean eq 4) then begin ; on ne fait rien : l'image finale est une image residuelle. endif if(mode_clean eq 5) then begin ; add comp synth int a dom_c. toto = where( (dom_c - indice) eq 0, ic) if(ic gt 0) then begin im_s = im_s + gain * im_r_kl * shift(f_remp, x, y) endif endif ; if ((i_ech eq 3) and(jc1 ge 105)) then begin ; print, 'i_ech =', i_ech, ' jc1 =', jc1, $ ; ' 1-1' ; stop ; endif ; <- niveau de la boucle de recherche de jmax composantes sans ; tests d'arret ; Calcul du maximum de l'image residuelle et de sa valeur quadratique ; moyenne, tous deux sur le domaine ou on recherche les com- ; posantes "sales", en vue des tests d'arret. ; Ce domaine est : ; . tout le champ pour mode_clean = 1 2 ou 3 (bien que ; dans ce dernier cas dom_c soit plus petit). ; . dom_c pour mode_clean = 4 ou 5 . max_r(j_par) = im_r_kl ; max_r est maximum sign'e d'image residuelle sur dom_c. ; Pour i=0 c'est le maximum de l'image initiale. ; Il est inutile de distinguer suivant les valeurs de ; mode_clean puisque la recherche du maximuma deja ete ; limitee a dom_c . if((mode_clean eq 1) or (mode_clean eq 2) or (mode_clean eq 3)) $ then begin ; mom = moment(im_r) ; Peut produire un underflow pour mom(0). Pas grave. ; Abandonne le 17 oct 02 au profit du calcul direct sui- ; vant. Le temps moyen par iteration passe alors de ; 7.2 msec a 2.8 msec !! le calcul direct suivant prend ; environ 0.25 msec sur mesopz. ; quad_r(j_par) = total(im_r)^2 / (nx * ny) ; Abandon (17 oct 02) du calcul direct qui ne sert que dans ; le cas iseuil=2 (plus jamais le cas depuis le 24 mai 02. endif if ((mode_clean eq 4) or (mode_clean eq 5)) then begin ; mom = moment(im_r(dom_c)) ; quad_r(j_par) = total(im_r(dom_c)^2) / n_elements(dom_c) ; abandom calcul direct 17 oct 02 (idem 5 lignes plus haut) endif ; <- niveau de la boucle de recherche de jmax composantes sans ; tests d'arret ; Controle d'essai ; if ((i_ech eq 2) and (j_cum gt 320) and (j_par ge 2)) then begin ; print, ' 1er dans boucle jmax : j_cum =', $ ; j_cum, ' j_par =', j_par ; shade_surf, im_r ; stop ; endif ; if ( (i_ech eq 2) and (j_cum gt 320) and (j_par ge 2) ) then begin ; print, ' 2eme dans boucle jmax : j_cum =', $ ; j_cum, ' j_par =', j_par ; print, ' mom =', mom ; stop ; endif ; quad_r(j_par) = sqrt(mom(1) + mom(0)^2) ; remplace 17 oct 02 par calcul direct pour gagner du temps. endfor ; fin de recherche de jmax composantes "sales" (indice j_par). ; Detection d'underflows icon = 0 if (icon eq 1) then begin if ( (i_ech eq 1) ) then begin print, ' i_ech =', i_ech print, ' Fin de boucle jmax : j_cum =', j_cum stop endif endif ; fin de recherche d'underflows. ; Controle de la position des composantes clean icon = 0 if ((icon eq 1) and (i_ech eq 0)) then begin print, 'CLEAN_ECHELLE : i_ech =', i_ech window, 0 plot, max_x xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Absc pos des max dans une serie de jmax' window, 1 plot, max_y xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Ord pos des max dans une serie de jmax' window, 2 plot, max_r xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Max d image residuelle dans une serie de jmax' print, 'Absc des max d image residuelle' print, max_x print, 'Ord des max d image residuelle' print, max_y print, 'Amplitude des max d image residuelle' print, max_r stop endif ; TESTS D'ARRET ; Arret quand la moyenne (sur jmax images residuelles) des maxima en valeur ; absolue ou des valeurs rms de l'image residuelle sur le domaine ; dom_c devient < seuil_amp (fraction du max de valeur absolue ; ou de la valeur rms de l'image initiale). ; Ce critere devrait arreter l'iteration pour les plus grandes ; echelles avant un trop grand nombre d'iterations. En effet dans ; le cas des plus grandes echelles l'iteration conduit a une ima- ; ge residuelle de type "plateau" qui baisse lentement et il faut ; un tres grand nombre d'iterations pour que les composantes ; CLEAN changent de signe, alors que l'image residuelle filtree ; d'echelle est devenue tres faible. Ce critere permet alors ; d'arreter avant un grand nombre d'iterations inutiles. ; Calcul preliminaire du maximum d'amplitude de l'image residuelle. if (i_seuil eq 1) then begin ; seul cas depuis 24 mai 02 ; rappel : max_r est maximum sign'e d'image residuelle ; sur dom_c. if (jmax eq 1) then begin amp_moy = abs(max_r(0)) ; max_r a un seul element. endif if (jmax ge 2) then begin ; mom_max_abs = moment(abs(max_r)) ; vecteur des moments des jmax maxima en val absolue ; d'images residuelles sur dom_c. ; moy_max_abs_res = mom_max_abs(0) moy_max_abs_res = total(abs(max_r)) / jmax ; moy_max_abs_res est la moyenne des jmax maxima ; en val absolue de l'image residuelle sur dom_c. ; abandon de "moment" pour ganer du temps (17 oct 02) amp_moy = moy_max_abs_res ; apres 29 mai 02. ; mom_max_sign = moment(max_r) ; jusqu'au 29 mai 02 ; vecteur des moments des jmax maxima sign'es d'images ; residuelles sur dom_c. ; moy_max_sign_res = mom_max_sign(0) ; moy_max_sign_res est la moyenne des jmax maxima ; sign'es (maxima en val absolue mais affectes de ; leur signe) de l'image residuelle sur le domaine ; dom_c. ; amp_moy = moy_max_sign_res ; jusqu'au 29 mai 02. pouvait ; etre beaucoup plus petit ; que ses elements et arreter ; l'iteration trop tot. endif endif if (i_seuil eq 2) then begin ; jamais depuis 24 mai 02. mom_quad = moment(quad_r) ; vecteur des moment des jmax valeurs quadratiques moyen- ; nes d'image residuelles sur dom_c. moy_quad_res = mom_quad(0) ; moy_quad_res est la moyenne des jmax valeurs quadrati- ; ques moyennes sur le domaine dom_c des jmax images resi- ; duelles . amp_moy = moy_quad_res endif if (crit_arret(0) eq 1) then begin if (abs(amp_moy) le seuil_amp) then begin ; rem : amp_moy est necessairement >0 depuis le 29 mai 02. drap_iter = 0 ; => arret du while j_tot = j_cum + jmax ; j_tot (utilise pour l'impression locale) est egal a j_cum ; apres son increment juste avant le endwhile. if (i_test_arret eq 1) then begin ch_format = "(' Filtrage ', i2, i6, ' iterations:" + $ " passage max image residuelle sous seuil " + $ "ampl.')" print, format=ch_format, i_ech, j_tot if (i_seuil_max_rms eq 1) then ch_format = $ "(20x, 'd apres max, fact mult du seuil " + $ "d arret =', f5.2)" if (i_seuil_max_rms eq 2) then ch_format = $ "(20x, 'd apres rms, fact mult du seuil " + $ "d arret =', f5.2)" print, format=ch_format, fac_seuil endif endif endif ; Arret si le flux de l'image residuelle passe sous une fraction du flux ; de l'image initiale ; Ce critere d'arret est limit'e au cas i_ech = 0 , le seul ; ou les images ont un flux non nul. Le flux des images polar ; peut etre negatif. if ((crit_arret(1) eq 1) and (i_ech eq 0)) then begin flux = total(im_r) if (abs(flux) lt seuil_flux) then begin drap_iter = 0 j_tot = j_cum + jmax if (i_test_arret eq 1) then begin ch_format = "(' Filtrage ', i2, i6, ' iterations:" + $ " flux residuel < ', f6.4, ' flux initial.')" frac = abs(seuil_flux/ flux_initial) print, format=ch_format, i_ech, j_tot, frac endif endif ; fin de la decision d'arret. endif ; fin du test special pour i_ech=i_ech_max. ; Rafraichissement en vue test sur croissance du max d'image residuelle ; (ce test est fait ci dessus) a l'iteration suivante. if (i_seuil eq 1) then amp_preced = moy_max_abs_res ; moy_max_sign_res jusqu'au 29 mai 02. ; Arret si le max des images residuelles (sur dom_c) ou leur valeur quadra- ; tique moyenne, suivant la valeur de i_seuil, commence a croitre ; nettement en valeur absolue. Rappel : ; . moy_max_sign_res est la moyenne des jmax maxima sign'es (maxi ; ma en val absolue mais affectes de leur signe) de l'image ; residuelle sur le domaine dom_c. if (crit_arret(2) eq 1) then begin ; if (i_seuil eq 1) then amp_moy = moy_max_sign_res ; -> 7 jun 02. if (i_seuil eq 1) then amp_moy = moy_max_abs_res if (i_seuil eq 2) then amp_moy = moy_quad_res ; rappel iseuil est toujours 1 depuis le 24 mai 02 if (abs(amp_moy) gt fac_crois_res*abs(amp_preced)) then begin drap_iter = 0 j_tot = j_cum + jmax if (i_test_arret eq 1) then begin ch_format = "(' Filtrage ', i2, i6, ' iterations:" + $ " croissance de l image residuelle.')" print, format=ch_format, i_ech, j_tot endif endif endif ; Arret si, dans une serie de jmax iterations, la proportion des maxima ; d'image residuelles exhibant un changement de signe excede ; une valeur donnee (0.1 a 0.4). L'idee est que la soustraction ; d'une composante clean acheve de faire disparaitre le dernier ; pic dominant. La composante clean suivante concerne alors un ; secondaire negatif oubli'e ou du bruit. ; Remarque : ce critere n'est clair que si l'image a un signe domi- ; nant, c'est a dire qu'elle a un flux non nul. Ce n'est le cas ; que pour les plus grandes echelles (i_ech=0). Cependant on ; l'utilise quand meme car cela donne de bons resultats quand ; la calibration est mauvaise avec un gros secondaire negatif ; genre anti-lobe : cela interdit de nettoyer ce secondaire, et ; si en outre on ne rajoute pas le residu a l'image finale (avec ; mode_clean=2) le gros seondaire negatif n'existe plus sur l'ima ; ge finale (essai du 6 mars 02 sur le soleil du 01 jan 02 a ; 09:08:40). ; Remarque : ce test parait meilleur que celui fait plus bas et qui ; utilise des valeurs moyennees sur une serie de jmax iterations. ; Dans de dernier cas la valeur moyenne peut etre beaucoup plus ; petite que les max_r a cause des chgts de signe. ; Rappel : max_r est le max sign'e de jmax images residuelles suc- ; cessives . max_r(0) = max d'image initiale du paquet, ; cad de l'image initiale s'il s'agit du 1er paquet. ; Rappel : On traite d'abord les grandes echelles (i_ech=0) depuis ; le 12 jan 01. if (crit_arret(3) eq 1) then begin rap_stop = 0.1; 0.1 ; de quand jmax etait grand. toto = where(im_max*max_r gt 0, n_s_id) ; pour "signe idem". toto = where(im_max*max_r le 0, n_s_op) ; pour "signe oppose". if(n_s_op ne 0) then begin if (n_s_id eq 0) then goto, a1 ; tous signes chang'es. if (float(n_s_op)/float(n_s_id) lt rap_stop) then goto, a2 a1: drap_iter = 0 j_tot = j_cum + jmax if (i_test_arret eq 1) then begin ch_format = "(' Filtrage ', i2, i6, ' iterations: " + $ "signe des composantes : ', i2,' idem', i4, ' oppose')" print, format=ch_format, i_ech, j_tot, n_s_id, n_s_op endif endif a2: j_M = 0 ; a2 non admis devant endif, endif ; Arret si le signe du (max sign'e d'image residuelle, moyenn'e sur ; jmax iterations) devient oppose a celui de l'image initiale. ; MIS EN COMMENTAIRE le 7 dec 00 apres le test sur le debut de ; l'apparition des chgts de signe dans une serie max_r. if ( (i_seuil eq 1) and (im_max*amp_moy lt 0) ) then begin ; j_tot = j_cum + jmax ; drap_iter = 0 if (i_test_arret eq 1) then begin ; ch_format = "(' Filtrage ', i2, i6, ' iterations:" + $ ; " chgt de signe du max d image residuelle.')" ; print, format=ch_format, i_ech, j_tot endif ; stop endif ; Arret si l'image devient plus negative que positive, c'est a dire si ; abs(min) > max sur l'image. Critere douteux pour la polar. ; De toutes facons c'est fait en mieux avec le critere sur le ; chgt de signe du max absolu. LAISSE EN COMMENTAIRE. ; if (max(im_r) lt abs( min (im_r))) then begin ; max_im_r = 0. ; max_res_jmax = 0.8 * seuil_amp ; endif ; FIN DES TESTS D'ARRET. ; Actualisation de l'indice d'iteration pour nouvelle serie de jmax. j_cum = j_cum + jmax ; -> j = 0L, jmax, 2*jmax, etc. ; Sauvegarde des images synthetique et residuelle avant nouvelle iteration if (drap_iter eq 1) then begin ; print, 'RH_CLEAN : filtrage', i_ech, ' sauvegarde' im_s_1 = im_s im_r_1 = im_r endif endwhile ; fin d'iteration de soustraction de composantes CLEAN (indice ; j_cum). ; Fin de recherche de composantes clean. ; Temps d'execution t4 = systime(1) ; temps en sec depuis le 1er jan 70. ch_1 = "(' CLEAN_ECHELLE : temps moyen par iteration = ', " ch_2 = "' msec')" ch_format = ch_1 + " f5.2," + ch_2 ; print, format=ch_format, 1000 * (t4 - t0) / j_cum ; Resultats du 17 oct 02 sur mesopz : ; A) programme dans sa version -> 17 oct 02 : ; . temps moyen pour une iteration avec jmax=100 = 8.3 msec ; . ----------------------------------- jmax= 1 = 8.7 ---- ; On en deduit que les tests d'arret prennent 0.4 msec. ; B) apres transfert hors de la boucle des fonctions de soustrac- ; tion et de substitution (gardant jmax=2) : gain de 1.4 msec ; (8.6 -> 7.2 msec) ; C) abandon de la fonction IDL "moment" qui calcule des moments ; d'ordre 3 et 4 qu'on n'utilise pas. Calcul direct de la ; valeur quadratique moyenne. Gain tres important : ; 7.2 msec -> 2.7 msec ; D) abandon complet de la voute (l'alias peut etre plus proche du ; centre que l'orage de bruit). Le temps par iteration passe ; de 2.7 msec a 1.90 msec . ; E) abandon du calcul direct de la valeur quadratique moyenne de ; l'image residuelle, qui ne sert plus depuis que i_seuil est ; toujours 1. Le temps par iteration passe de 1.90 msec a ; 1.60 msec . ; F) abandon de "moment" dans le calcul preliminaire du maximum ; d'amplitude de l'image residuelle. Le temps par iteration ; passe de 1.60 msec a 1.5 msec ; Le gain cumul'e est 8.6 / 1.5 = 5.7 ; Controle de l'arret par changement de signe des composanates clean OBSOLETE icon = 0 if (icon eq 1) then begin print, ' Filtrage ', i_ech, ' test arret par chgt signe' ch_format = "(' nbre iterations =', i4, ' n_s_id =', i2, " + $ "' n_s_op =', i2)" print, format=ch_format, jc1, n_s_id, n_s_op ; stop endif ;print,j_cum ; Affichage d'arret par depassement du nbre d'iterations permises, if (j_cum ge itermax) then begin ch_format = "(' Filtrage ', i2, i6, ' iterations" + $ ": depassement du nbre d iterations permises.')" print, format=ch_format, i_ech, j_cum endif ; if (i_ech eq 2) then stop ; recherche d'underflows. ; Restauration des images synthetiques et residuelles avant les iterations ; ayant ayant provoque l'arret ; On considere que ces composantes clean ne sont pas significatives ; (par exemple quand les composantes clean changent de signe du fait ; d'une mauvaise calibration). On revient donc aux images precedant ; ces iterations. On peut toutefois decider de sauter cette restaura ; tion. i_saut = 0 ; 0 POUR RESTAURER if (i_saut eq 0) then begin im_s = im_s_1 im_r = im_r_1 endif ; Controle: comparaison images sale et propre apres iterations icon = 0 if (icon eq 1) then begin ch_i_ech = string(format="('iech = ', i1)", i_ech) ; window, 0 & wset, 0 ; print, 'Nombre d iterations =', j ; shade_surf, congrid(im , 512, 512, cubic=-0.5) ; xyouts, 0.5, 0.97, /normal, alignment=0.5, $ ; 'Image filtree d echelle initiale' ; window, 1 & wset, 1 ; shade_surf, congrid(im_r, 512, 512, cubic=-0.5) ; xyouts, 0.5, 0.97, /normal, alignment=0.5, $ ; 'Image filtree d echelle residuelle' ; window, 2 & wset, 2 ; shade_surf, congrid(im - im_r, 512, 512, cubic=-0.5) ; xyouts, 0.5, 0.97, /normal, alignment=0.5, $ ; 'Diff images filtrees d echelle initiale - residuelle' ; window, 3 & wset, 3 ; shade_surf, congrid(im_s, 512, 512, cubic=-0.5) ; xyouts, 0.5, 0.97, /normal, alignment=0.5, $ ; 'Image filtree d echelle propre' ; stop print, 'CLEAN_ECHELLE : resultats avant lissage (dans RH_CLEAN)' loadct, 0 window, 0 & wset, 0 shade_surf, congrid(im , 512, 512, cubic=-0.5) xyouts, 0.5, 0.97, /normal, alignment=0.5, charsize=1.3, $ ch_i_ech + ' : image filtree d echelle initiale' window, 1 & wset, 1 shade_surf, congrid(im_r, 512, 512, cubic=-0.5) xyouts, 0.5, 0.97, /normal, alignment=0.5, charsize=1.3, $ ch_i_ech + ' : image residuelle' window, 2 & wset, 2 shade_surf, congrid(im_s + im_r, 512, 512, cubic=-0.5) xyouts, 0.5, 0.97, /normal, alignment=0.5, charsize=1.3, $ ch_i_ech + $ ' : somme images filtrees d echelle propre et residuelle' stop endif ; Fin de controle im_s = float(im_s) ; precaution. end ; fin de CLEAN_ECHELLE. ; Musee delouisien: ; print, "p_clean 2d : the clean algorythm (sic) start (re-sic) ", $ ; "diverging at iter =", j ;____________________________________________________________________________ PRO MAX_RMS, i_ech, im, $ ; entree max_im, rms_im ; sortie ; Creation 1 mai 2002. Test'e avec TEST_MAX_RMS ; But : calculer les valeurs maximum et rms d'une image interferometrique ; filtree d'echelle sur un domaine raisonnablement du meme ordre que ; l'image solaire et non nettement plus grand comme c'est le cas pour le ; champ interferometrique, specialement a basse frequence. ; En vue de calculer le seuil d'arret de l'iteration clean. ; Role de i_ech : ; -1 pas de filtrage d'echelle. Il s'agit d'un essai sur un objet ; simul'e. ; 0, 1, 2, 3 l'image en entree est filtree d'echelle. taille = size(im) npi = taille(1) ; Definition d'un domaine entourant le centre du champ, suffisamment reduit ; pour ne pas deborder du soleil calme, meme aux frequences basses ; et assez grand pour n'etre pas recouvert par une source compacte ; (pouvant etre intense). ; Si l'image n'est pas filtree d'echelle, s'il s'agit d'un soleil calme ; simul'e, il peut y avoit une grande partie du champ ou elle est ; quasi nulle et sans signification. ; Si l'image est filtree d'echelle, elle est nettement plus delocali- ; see. Les simulations montrent que pour i_ech=0 et 1 l'image oscille ; partout, avec eventuellement un pic dominant. Pour i_ech=2 et 3 ; elle occupe encore toute la largeur en EW du fait de l'aliasing ; mais elle est resreinte a une bande EW plus limitee en NS, dont la ; hauteur depend evidemment de la frequence d'observation. ; Conclusion sur la definition du domaine : ; . i_ech=-1 : on prend un rectangle au centre. ; . i_ech=0 ou 1 : -------- tout le champ. ; . i_ech=2 ou 3 : -------- une bande centrale limitee en NS. if (i_ech eq -1) then begin ; image non filtree d'echelle. l2_u = 25 ; 1/2 largeur en u l2_v = 15 ; -------------- v, peut-etre un peu grande en hiver ; au meridien. endif if ((i_ech eq 0) or (i_ech eq 1)) then begin l2_u = 63 l2_v = 63 endif if ((i_ech eq 2) or (i_ech eq 3) or (i_ech eq 4)) then begin l2_u = 63 l2_v = 32 ; imparfait car ca devrait dependre de la frequence. endif im_red = im(npi/2 - l2_u : npi/2 + l2_u, npi/2 - l2_v : npi/2 + l2_v) n_pts = (2 * l2_u + 1) * (2 * l2_v + 1) rms_im = sqrt( total(im_red^2) / n_pts) rms_0 = rms_im ; sauvegarde. max_im = max(im) ; sur tout le champ. ; Detection du cas ou il existe une ou plusieurs sources compactes nettement ; plus intenses que le niveau moyen ; Le nbre de pts retenu sur l'image est environ 50 * 30 = 1500. Une ; source compacte solaire a une surface typique de 10 - 20 pixels ; soit 20/1500 = 1/75 de l'aire definie par l2_u et l2_v. Une telle ; source compacte commence a dominer rms_im si son rapport au niveau ; moyen est > sqrt(75) ~ 9. Pour calculer une valeur rms typique du ; fond, il faut exciser cette source et recommencer le calcul de ; l'amplitude rms. ; Rem du 16 mai 02 : tout ca n'est vrai que pour une image non filtree ; d'echelle. Si elle est filtree d'echelle le pic n'est jamais forte- ; ment dominant (meme pour une source ponctuelle) a cause des ; enormes secondaires. Pour des image filtrees la 2eme valeur de la ; moyenne rms n'est pas tres differente (abaissement de 25% au max) ; de la 1ere valeur, aussi bien pour une d'une source ponctuelle que ; pour un soleil complexe. max_rms = max_im / rms_im rap_crit = 9.0 ; on n'excise que si un pic domine vraiment. if (max_rms ge rap_crit) then begin ; 1) il faut retirer la source la ou elle excede le fond moyen ; multipli'e par ~ rap_crit/2. A ce stade le fond moyen, eva- ; lu'e par rms_im, peut etre biais'e mais on s'en contente. ; 2) pour des images filtrees d'echelle la condition ne sera pas ; remplie, en particulier pour i_ech=0 (le "lobe" est tres ; large) et i_ech=1 (gros secondaires). dom = where(im_red le (rap_crit/2)*rms_im) nbre_min = (2*l2_u * 2*l2_v) / 4 ; on admets qu'il faut faire la moyenne sur un nbre suffisant de ; points, ici le quart du domaine initial. if (n_elements(dom) ge nbre_min) then begin im_red_fond = im_red(dom) n_red_fond = n_elements(im_red_fond) rms_im = sqrt(total(im_red_fond^2) / n_red_fond) ; help, amp_rel, rms_0, (rap_crit/2)*rms_0, dom, n_pts, $ ; n_red_fond, rms_im endif endif ; Verification limitable aux images non filtrees d'echelle) icon = 0 ; if (icon eq 1) then begin if ((icon eq 1) and (i_ech eq -1)) then begin print, 'MAX_RMS :' print, ' i_ech=', i_ech print, ' nbre d''elements du domaine :', $ n_elements(im_red) print, ' nbre d''elements du domaine excise du max :', $ n_elements(dom) print, ' max image =', max_im, ' 1ere val rms =', rms_0 print, ' 2eme val rms =', rms_im stop endif end ; fin de MAX_RMS. ;______________________________________________________________________________ PRO RH_CLEAN, $ ; entree: ; donnees image freq, ch_pol, $ i_simul_soleil, i_pave_plein, $ tf_objet_simule, objet_simule, $ tf_im, im, tf_l_th, l_th, $ ; clean mode_clean, dom_c, dom_a, $ i_ech_min, i_ech_max , a_in, a_ex, $ a_red_fil, itermax, gain_iter, jmax, crit_arret, $ i_seuil, coeff_seuil, frac_flux, fac_crois_res, $ i_test_arret, i_lissage, poids_c, larg_2, $ ; sortie tf_lobe_synt, lobe_synt, im_p, $ objet_simule_lisse ; Indications sommaires sur les variables ; . tf_im et tf_l_th sont centrees en milieu de tableau. ; . affichage et controle seulement : ; freq, ch_pol, i_pave_plein, ; . simulation (autrement varaibles bidon) ; i_simul_soleil, tf_objet_simule, objet_simule ; MODIFICATIONS (y compris celles de l'ancien WCLEAN, incorpore dans RH_CLEAN : ; 00 dec 5 introduction de : ; - en entree dom_c, domaine autour d'un centre d'orage dont ; veut debarasser l'image, et auquel est limite la recher- ; che des composantes clean ; - en sortie production de im_orage, image propre de l'orage ; sans ajout de residu. ; 01 sep 21 ajout de objet_simule dans la liste d'appel pour transmission ; a RH_CLEAN (alors WCLEAN) ; oct 2 ajout de i_pave_plein dans la liste d'appel pour transmission ; a RH_CLEAN (alors WCLEAN) ; oct 18 regroupement de REGCLEAN et WCLEAN pour faire RH_CLEAN. Abandon ; de toutes les traces de delouiseries. ; Modifications de l'ancien WCLEAN : ; 99 jun 2: TF directe (drapeau -1) pour passer des visibilites aux images ; (meme convention que dans MALC_IM_2D). ; 99 jun 24: normalisation du filtre d'echelle dans l'espace de Fourier. ; 99 sep 27: Definition (dans FILTRE_ECH) de 4 zones de filtrage annulaires ; reliees au remplissage du plan (u,v) et non plus au lobe gaus- ; sien. Le domaine radial est toujours environ 2 mais l'anneau ; exterieur a pour axes 64 en EW et 50 en NS (environ 2 fois plus ; que dans le cas J-M. ; 99 oct 3: Le nbre des zones de filtrage passe a 5 pour eviter un trop large ; domaine d'echelles dans l'anneau exterieur. ; 99 oct 6: Realisation d'un leger lissage apres le CLEAN. Dans le cas d'une ; image 128*128 ce lissage doit etre leger sinon on perd la ; resolution apportee par les extensions. pour faire un lissage ; gaussien il faut ume image au moins 256*256. ; 00 nov 30: Transformation de function en procedure. ; Ajout de dom_c en argument pour limiter le domaine des positions ; des composantes clean recherchees. ; 00 dec 5: Dans le cas du retrait de la contribution d'un centre d'orage ; localis'e, fourniture en sortie de l'image propre sans residu ; de ce centre d'orage (dernier argument dans la liste im_orage) ; 8: Reduction du nombre de points sur les images filtrees d'echelle ; qui sont tres sur-echantillonnes pour les grandes echelles ; (parametre i_ech =3 ou 4). Ceci dans l'espoir de reduire le ; nombre d'iterations. ; dec 12: Clarification de la normalisation de l_th (normalise a un maximum ; de 1 et de i_fil (TF de "filtre"*total(l_th)). ; Normalisation de l_th_ech a un max de 1, correspondant a son usa- ; ge dans CLEAN_ECHELLE (soustraction de composantes clean), et ; normalisation correspondante de i_fil. ; 00 dec 20: grand nettoyage d'hiver : suppression de toutes les delouiseries, ; widgeteries etc. ; 01 jan 4: visualisation du contour du domaine dom_c. ; 9: simulation possible d'un soleil calme bien merdique pour tester ; le clean. ; 11: i_ech=0 correspond aux grandes echelles (et plus aux petites). ; 01 sep 25: calcul de seuil_amp inclut le facteur rms_im_nf/rms_im_ech pour ; limiter le nombre des iterations pour les petites echelles dans ; le cas des images assez lisses. ; oct 2: ajout de i_pave_plein dans la liste d'appel pour affichage de ; controle. ; oct 18: eradication de tous les vestiges delouisiens, en particulier im_f ; (accumulateurs de filtres d'echelles) et de tout ce qui est ; li'e au lobe gaussien lg. ; 02 mar 6: option de ne pas rajouter l'image residuelle a l'image synthetique ; (mode_clean=2), ; on etend (dans CLEAN_ECHELLE) le critere d'arret sur le changement ; de signe des composantes clean au cas de toutes les echelles. ; Cela donne de bons resultats pour supprimer les gros secondai- ; res negatifs dus a une mauvaise calibration, si on n'ajoute pas ; l'image residuelle a l'image synthetique (mode_clean=2). ; annulation (dans CLEAN_ECHELLE) de la soustraction des composantes ; clean de l'image residuelle et de leur ajout a l'image synthe- ; tique pour les dernieres iterations (ayant provoque l'arret). ; 02 avr 17: utilisation de axes (ajout'e a la liste de sortie de FILTR_ECH) ; pour definir le lissage final : le filtre de fourier tombe a ; zero au bord. On realise le 18 avril que c'est probablement ; nuisible puisqu'il y a un filtrage final gaussien (commentaires ; vers ligne 1880). ; 02 mai 15: calcul du seuil d'arret : le plus grand des deux seuils : ; . fraction (petite) du max de l'image ; . -------- (plus grande) d'une valeur typique du fond du SC, ; sources compactes intenses exclues ; Le critere d'arret porte toujours sur le max de l'image resi- ; duelle et non plus sur une valeur rms comme c'etait la cas pour ; les echelles moyennes et petites (dans le cas d'une source ; compacte le seuil rms arrete l'iteration trop tot). On conserve ; i_seuil bien qu'il soit maintenant inutile. ; 02 mai 24: la couronne externe doit etre contrainte par le bord du pave. Les ; essais du 24 mai 02 montrent que les resultats sur la TF de ; l'image propre non filtree sont ainsi bien meilleurs, surtout ; avec une source compacte. L'effet de a_red_fil=0.7 etait simple ; ment de reduire assez la couronne externe pour la contraindre. ; (modif faite dans FILTRE_ECHELLE). ; 02 jun 7: ajout de objet_simule_lisse a la liste de sortie pour comparaison ; dans RH_CALC_IMAGE_HELIO. ; 02 mai 14 - jun 17 : ; La dynamique sur une image nettoyee depend certainement de ; l'aspect de l'image. Il est intuitif que cette dynamique est ; plus importante quand l'image est dominee par une source com- ; pacte que quand c'est un magma complique. ; - dans le 1er cas on admet une fraction du max, mettons 1%, ; - dans le 2eme cas --------------------- de la valeur rms sur ; un domaine ou la source est importante (domaine plus petit ; que le soleil calme) mettons 5%. ; Dans le cas general on calcule la valeur rms de la source en ; excluant les sources compactes puissantes (ce que fair MAX_RMS) ; et on adopte comme seuil d'arret la plus grande de ces deux ; valeurs. ; Ca parait plus logique, pour les sources compliquees, que de ; comparer la moyenne rms du residu a celle de l'image initiale ; puisque la soustraction des composantes clean est essentielle- ; ment locale. ; ; En fait ce raisonnement intuitif ne vaut que pour des images ; non filtrees d'echelle. Une image non filtree peut consister en ; un fond de soleil calme bossele complexe plus ou moins domin'e ; par une ou plusieurs sources compactes localisees (beaucoup ; moins etendues que le SC). ; Dans les images filtrees d'echelle la base et les sources ; compactes ont des echelles moins separees (c'est d'ailleurs le ; but du filtrage d'echelle) : ; . pour i_ech=0, qui ne comporte que des tres bas harmoniques) ; l'image est tres large et on ne distingue pas entre base ; et source compacte. ; . pour i_ech=1 ou >1 il y a de forts secondaires et une source ; compacte est moins localisee que sur l'image non filtree. ; La procedure MAX_RMS (dont le but est de calculer les valeurs ; maximum et rms d'une image interferometrique filtree d'echelle) ; excise bien les sources compactes pour des images non filtrees ; et on obtient alors le max et une valeur rms typique du fond ; (source compacte exclue). ; Pour des images filtrees d'echelle MAX_RMS perd son efficacite ; Les essais du 14 juin 02 montrent que les valeurs max et rms ; sont essentiellement proportionnelles pour des sources consti- ; tuees d'un SC complexe (toujours la meme) et d'une source com- ; pacte dont on fait varier la brillance de 1 a 100 fois celle ; du SC. ; => pour des images filtrees d'echelle on ne peut prendre comme ; seuil d'arret qu'une fraction du max et non de la valeur ; moyenne rms (conclusion du 14 juin 02). ; Toutefois on prendra (automatiquement) un seuil 3 fois plus ; faible s'il existe une source compacte tres dominante. ; 02 oct 17 modifications dans CLEAN_ECHELLE pour diminuer le temps de calcul ; par iteration (8.6 msec pour la version jusqu'au 17 oct 02). ; Voir a "temps d'execution" dans CLEAN_ECHELLE. ; 02 oct 22: - ajout de freq et ch_pol a la liste d'entree pour affichage. ; - facteur freq/164 appliqu'e a seuil_amp pour tenir compte du fait ; que les hautes frequences sont plus bruiteuses et eviter un ; nombre d'iterations trop elev'e. Ce n'est qu'une rustine. ; 03 mar 20: - ajout de tf_im et tf_l_th a la liste d'entree. ; 06 fev 23: amlioration generale et simplification: ; - les limites des zones de filtrage sont exactement circulaires, ; ce qui m'a evit'e de me fatiguer pour m'assurer que la somme de ; filtres elliptiques est exactement egale a 1. ; - les limites des zones de filtrages successives sont exactement ; dans un facteur 2 (ce qui n'etait pas exactement le cas, alors ; que les largeurs des transitions l'etaient). Il en resulte que ; la somme des filtres est maintenant exactement egale a l'unite. ; - les transitions sont adoucies (non plus lineaires, mais en 1/2 ; periode de cosinus) et elargies au maximum, de facon que les ; filtres n'aient plus de plateau strict (il est indiqu'e dans ; rh_dpnew que a_in et a_ex doivent etre 2/3 et 4/3). ; - il est conseill'e dans rh_dpnew de figer aussi larg_2 = [18, 25] ; => on arrive a des artefacts ~ 1-2%, dus a la couronne externe, et ; sans doutes dus a la nature de la couverture uv au dela du ; pave central. ; 06 fev 24: - calcul du lobe synthetique propre lobe_synt et incrustation ; sur l'image propre. ; Notations ; entree: ; a_in a_ex definissent les transitions progressives internes et externes ; des zones annulaires de filtrage. Usuel 0.85 et 1.15 . ; Choisis dans RH_DPNEW. ; dom_c domaine de recherche des composantes clean. Retrait d'orage ; si dom_c est plus petit que l'image. ; dom_a domaine du contour de dom_c pour visualisation. ; l_th lobe theorique (calcul'e comme la TF de "1" puis normalis'e ; `a un max de 1 dans MALC_IM_2D. Cette normalisation est uti- ; le pour l'usage qui en est fait dans CLEAN_ECHELLE (soustrac- ; tion de composantes clean). Elle est differente de celle de ; im depuis le 11 dec 00. ; im image a nettoyer en entree , centr'ee en milieu de tableau. Non ; modifiee dans WCLEAN. ; l_th lobe theorique, centre milieu de tableau, sans apodisation, ; normalis'e dans MALC_IM_2D a un total de 1. ; reg pour compenser couverture (u, v) irreguliere. ; mode_clean choix du mode de clean : ; . 0 pas de clean ; . 1 pour clean normal sur champ entier. L'image finale est la ; somme de l'image synthetique et de l'image residuelle. ; . 2 idem mais on n'ajoute pas l'image residuelle a l'image ; synthetique. ; . 3 pour chercher partout les composantes clean mais ne pas ; restituer celles interieures au domaine dom_c defini ; par x_dom, y_dom et dr_dom. ; . 4 pour retirer les composantes clean de dom_c et obtenir ; une image residuelle non nettoyee du reste. ; . 5 pour rechercher les composantes sales dans dom_c et les ; restituer sous forme propre (nettoyage rapide d'un ; centre intense ; i_ech_min, definissent les zones de filtrage d'echelle : ; i_ech_max, . nuls tous les deux : aucun filtrage d'echelle ; . i_ech_max > i_ech_min filtrage d'echelle 1 disque central ; et 3 zones annulaires., ; Le 18 oct 01 on abandonne le dernier pan de la saga deloui- ; sienne avec i_filtr_ech = 2 (filtrage 4 zones et criteres ; d'arret obscurs). ; interne: ; l_th_ech lobe theorique filtr'e d'echelle, calcule et normalise a un ; max de par multiplication de la TF de l_th par "filtre". ; Son total est nul sauf pour les echelles les plus grandes pour ; lesquelles "filtre" inclut la composante continue. ; La "fonction de nettoyage" de l'image sale, pour en soustraire ; les sources compactes, doit avoir un maximum de l'ordre de 1 ; pour retirer ces composantes "a la louche" et non "a la petite ; cuiller", de facon a eviter un trop grand nombre d'iterations. ; On pourra ainsi utiliser l_th_ech/max(l_th_ech) dans l'equa- ; tion de Thompson utilisee dans CLEAN_ECHELLE. ; filtre Filtre d'echelle : anneau dans le plan (u, v), avec un gain max de ; 1 et des bords adoucis tels que le recouvrement des bords de 2 ; anneaux successifs soit egal a 1. La somme des filtres sur tout ; le plan (u,v) est donc constante et egale a 1. ; i_fil TF(filtre), image du filtre d'echelle (centre en milieu de ta- ; bleau comme l_th). ; seuil_amp valeur (en max absolu ou valeur rms, suivant la valeur de ; i_seuil)) de l'image residuelle pour laquelle on arrete l'ite- ; ration CLEAN. ; seuil_flux valeur du max (en valeur absolue) du flux d'image residuelle ; pour laquelle on arrete l'iteration CLEAN. ; l_th_ech resultat du filtrage d'echelle sur le lobe theorique l_th (lui- ; meme normalise a un maximum de 1 depuis le 11 dec 00). On ; rappelle que l'amplitude du filtre est 1 dans sa bande. La som- ; des points de l_th_ech est donc nulle sauf pour i_ech=4, ou le ; filtre passe la composante continue a l'origine. ; l_th_ech est centre en milieu de tableau comme i_fil et l_th. ; sortie ; im_p image nettoyee en sortie. im = float(im) l_th = float(l_th) tab = size (im) nx = tab (1) ny = tab (2) ; tab est inutilise plus bas. if (nx ne ny) then begin print, 'RH_CLEAN : attention image en entree non carree' print, ' dimensions :', nx, ny endif np = nx ; Controle icon = 0 if (icon eq 1) then begin print, 'entree dans RH_CLEAN: ' print, max(float(im)), min(float(im)) print, max(float(l_th)) , min(float(l_th)) shade_surf, im help, champ stop endif ; Fin de controle ; Parametres a choisir dans RH_DPNEW : ; . a_red_fil (0.7 soleil calme bossele ou 1.0 pour types III forts) ; . itermax = 8000 ou moins ; . gain_iter = 0.2 (gain d'iteration). ; . frac_flux = 0.001 a 0.02 (fraction de flux d'image initiale sous ; lequel on arrete l'iteration/ ; . larg_2 = [15, 20] largeurs a mi puissance du filtre gaussien ; final. ; usuel [150, 20] pour soleil calme bossele , [24, 36] au max ; pour une source genre type III complexe intense. ; Rappel sur le choix des parametres de clean choisis dans dpnew depuis le 14 ; fev 01 ; Commentaire refait le 23 jan 01 apres avoir rendue progressive la ; transition externe de la zone externe. ; Les valeurs proposees dans dpnew ont ete choisies apres les essais ; suivants : ; 1) Soleil (simule par SIMUL_VISIB_2, contenu dans RH_SIMULATION) : ; base, couloir sombre de filament, 2 bacilles, 15 sources epar- ; ses, le tout fait avec des sources elementaires gaussiennes de ; 1/2 largeur a mi-puissance larg_c. ; Essais avec larg_c=3, 2, 0.2 . Dans tous les cas le meilleur ; resultat est obtenu avec a_red_fil=0.7 et larg_2=24. ; Pour a_red_fil > 0.7 a 0.75 il y a trop de details inventes. ; Pour larg_2 > 24 il y a des anneaux visibles autour des sources. ; 2) Source genre type III constituee de 5 sources avec larg_c=0.2 et ; disposees en 5 de domino. ; Meilleur resultat pour : ; . a_red_fil=1 (le taux d'anneaux augmente quand on reduit les ; couronnes de filtrage avec a_red_fil), ; . larg_2 = 24 ou au max 32 . Au dela la resolution de l'image ; propre n'augmente plus et lea anneaux sont de plus en plus ; visibles. ; Conclusion : ; . le seuil rms permet d'aller plus loin dans les iterations (et ; avec de meilleurs resultats) que le seuil sur la valeur max, ; pour toutes les domaines de filtrage autres que le domaine cen- ; tral. On prend donc i_seuil = [1, 2, 2, 2]. ; correctif du 1 mai 02 : dans le cas d'une source compacte les ; images filtrees d'echelle ont un max bien net et le seuil rms ; arrete l'iteration trop tot. On choisit donc plus bas le mode ; de calcul su seuil d'arret selon le caractere de l'image. ; . i_ech_max = 3, avec une couronne centrale pas trop petite, permet ; un bon nettoyage pour les grandes echelles avec une iteration ; limitee (env 300). ; . soleil calme bossele : a_red_fil = 0.7 larg_2 = 24 ; . genre type III : 1.0 24 ou 32 ; Initialisation de l'image finale "propre" et de l'accu des images des filtres im_p = fltarr(nx, ny) ; image propre finale. "somme" de J-M. im_r = fltarr(nx, ny) ; image residuelle finale. ; Remarque: fft(fft) = identique a 2E-7 pres ; (essai du 1 jul 99). ; Evaluation de l'importance des sources compactes par rapport au soleil calme ; (sur l'image non filtree d'echelle) i_ech_test = -1 ; pour absence de filtrage d'echelle. MAX_RMS, i_ech_test, im, $ ; entree. max_im, rms_im ; sortie. ; max_im est le max de la valeur ; absolue de l'image filtree ; d'echelle. ; rms_im est calculee sur un domaine ; important du soleil calme (sour- ; ces compactes exclues). rap_max_rms = max_im / rms_im rap_1 = 3 ; correspond a fac_mul = fac_max rap_2 = 30 ; -------------------- = 1 ; (pour une source compacte intense). ; Raccord du coeff de seuil entre les sources compactes intenses et faibles fac_max = 3 ; variation max cu coeff de seuil. if (rap_max_rms gt rap_2) then fac_mul = 1 ; sources compactes intenses if ((rap_max_rms lt rap_2) and (rap_max_rms gt rap_1)) then begin ; Dans le plan (log(fac_mul), log(rap_max_rms)) on raccorde les ; points (rap_1, fac_max) et (rap_2, 1) par une droite, ce qui ; donne une relation de puissance entre fac_mul et rap_max_rms. ; Rem : dans le diagramme log log, l'ordonnee du 2eme point ; (fac_mul=1) est nulle. aa = (-alog(fac_max)) / (alog(rap_2) - alog(rap_1)) bb = alog(rap_1) * (alog(fac_max)) / $ (alog(rap_2) - alog(rap_1)) + alog(fac_max) fac_mul = exp(bb) * rap_max_rms^aa endif if (rap_max_rms lt rap_1) then fac_mul = fac_max ; pour sources compactes faibles. ; print, 'RH_CLEAN : facteur multiplicatif du seuil :', fac_mul icon = 0 if (icon eq 1) then begin print, 'RH_CLEAN, fact mult du seuil d ampl' shade_surf, im xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Image non filtree d''echelle' print, ' rap_max_rms =', rap_max_rms print, ' aa =', aa, ' bb =', bb print, ' fac_mul =', fac_mul stop endif ; FILTRAGE AVEC DES DOMAINES D'ECHELLES SUCCESSIFS (un domaine = un facteur 2). ; Initialisation des indicateurs de niveau des images filtrees d'echelle rms_im_ech = fltarr(6) max_im_ech = fltarr(6) rms_im_nf = sqrt(total(im^2) / nx*ny) ; rms_im_nf = sqrt(total(im^2)) / nx*ny ; jusqu'au 1 mai 02. Erreur ; sans consequence car utili- ; s'e seulement pour impres- ; sion. ; Debut boucle sur domaines d'echelles. sigma_filtre = fltarr(nx, nx) ; accumulateur des filtres. for i_ech=i_ech_min, i_ech_max do begin ; for i_ech=i_ech_min, i_ech_min + 2 do begin ; On traite d'abord les grandes echelles spatiales (i_ech=0). ; Calcul du filtre d'echelle. taille = size(im) npew = taille(1) & npns = taille(2) FILTR_ECH, npew, i_ech, i_ech_max, a_in, a_ex, $ ; entree a_red_fil, $ ; entree ax_e0, filtre ; sortie ; Rem : "filtre" est centr'e dans les angles du tableau. max_f = max(abs(filtre)) sigma_filtre = sigma_filtre + shift(filtre, nx/2, nx/2); pour controle. ; if (i_ech eq i_ech_max) then begin ; dom_ext = where(sigma_filtre lt 0.95) ; pour des controles. ; ; on definit ainsi le domaine (en retrecissant un peu) ou le ; ; filtre n'est pas nul. ; endif ; Calcul de l'image du filtre, centree en milieu de tableau i_fil = float (fft (filtre, -1)) ; cbeam de J-M. i_fil = i_fil * total(l_th) ; normalisation du 12 dec 00, ; adaptee a l'usage de i_fil ; dans CLEAN_ECHELLE. i_fil = shift (i_fil, nx/2, ny/2) ; centrage. ; i_fil est l'image, centree en milieu de tableau dans l'espace ; reel, de "filtre". La difference avec l_th_ech (defini plus ; bas) est que dans ce dernier cas on applique "filtre" a un ; plan (u, v) tres incompletement rempli. ; Rappel: FFT directe (drapeau -1) se fait avec fact 1/N ; FFT inverse -------- 1 ------- sans -------- ; Rappel : FFT(FFT) = identique a 2E-7 pres (essai 1 jul 99). ; Controle du filtre et son image icon = 0 if (icon eq 1) then begin ch_i_ech = string(format="('iech = ', i1)", i_ech) help, filtre filtre_c = shift(filtre, nx/2, nx/2) window, 0, xsize=512, ysize=512 loadct, 13 tvscl, congrid(filtre_c, 512, 512, cubic=-0.5) loadct, 0 xyouts, 0.5, 0.95, /normal, alignment=0.5, charsize=1.3, $ ch_i_ech + ' filtre' window, 1 loadct, 0 shade_surf, congrid(i_fil, 512, 512, cubic=-0.5) xyouts, 0.5, 0.95, /normal, alignment=0.5, charsize=1.3, $ ch_i_ech + ' image du filtre' print, 'max de i_fil =', max(i_fil), $ ' total(i_fil) =', total(i_fil) stop endif ; Conclusion: tableau du max de i_fil ; i_ech 0 1 2 3 ; facon J-M 0.056 0.014 0.0035 0.0012 ; facon CM 0.440 0.140 0.026 0.008 ; total(i_fil) = 1 si i_ech=3, 0 pour i_ech=0,1,2. ; Fin de controle du filtre et son image ; Filtrage d'echelle sur l'image et le lobe theorique. ; Rappel : "filtre" (espace TF) en sortie de FILTR_ECH, est centre ; au coin inferieur gauche. im_ech = float(fft(shift(tf_im , -nx/2, -nx/2) * filtre, -1)) l_th_ech = float(fft(shift(tf_l_th, -nx/2, -nx/2) * filtre, -1)) im_ech = shift(im_ech , nx/2, nx/2) l_th_ech = shift(l_th_ech, nx/2, nx/2) ; Formulation d'avant le passage de tf_im et tf_l_th en entree de RH_CLEAN ; im_ech = float (fft(fft(im , 1) * filtre, -1)) ; ancien. ; l_th_ech = float (fft(fft(l_th, 1) * filtre, -1)) ; l_th_ech est deduit de l_th par filtrage, le filtre valant 1 ; dans sa bande passante. ; Rappels ; 1) l_th et im sont normalises dans MALC_im_2d de facon ; que leur totaux soient resp. 1 et le flux en sfu. ; 2) l_th_ech n'est plus normalise a 1 apres le 28 sept 99. ; 3) Les fonctions de soustraction et de remplacement dans ; CLEAN_ECHELLE sont respectivement : ; l_th_ech / max(l_th_ech) et i_fil / max(l_th_ech) ; depuis le 13 dec 00. ; 3) dans la saga delouisienne l_th_ech etait le 2eme ibeam. ; Controle du filtre, du lobe et de l'image filtr'es d'echelle. icon = 0 if (icon eq 1) then begin shade_surf, congrid(l_th , 512, 512, cubic=-0.5) window, 1 & wset, 1 shade_surf, congrid(l_th_ech, 512, 512, cubic=-0.5) window, 2 & wset, 2 shade_surf, congrid(im , 512, 512, cubic=-0.5) window, 3 & wset, 3 shade_surf, congrid(im_ech, 512, 512, cubic=-0.5) print, 't_th_ech : max =', max(l_th_ech), $ ' total = ', total(l_th_ech) stop endif ; ACCORD LE 19 AUG 99 AVEC LA VERSION JM. ; total(l_th_ech) = 1 si i_ech=3, 0 pour i_ech=0,1, 2. ; i_fil, l_th, l_th_ech sont centres en milieu de tableau. ; Fin de controle ; Calcul du max de l'image et d'une valeur rms representative du fond moyen ; (en vue de definir un seuil d'arret de l'iteration) im_ech_1 = im_ech if ((mode_clean eq 4) or (mode_clean eq 5)) then begin im_ech_1 = fltarr(np, np) im_ech_1(dom_c) = im_ech(dom_c) ; et le reste est nul. endif MAX_RMS, i_ech, im_ech_1, $ ; entree. max_im, rms_im ; sortie. ; max_im est le max de la valeur ; absolue de l'image filtree ; d'echelle. ; rms_im est calculee sur un domaine ; recouvrant un bon peu du soleil ; calme, a l'exclusion des sources ; compactes qu'il trouve tout seul ; sans utiliser dom_c (commentai- ; res dans MAX_RMS). max_im_ech (i_ech) = max_im rms_im_ech (i_ech) = rms_im ; Controle de max_im et rms_im dans MAX_RMS ; OBSOLESCENCE DEBUT ; Choix du mode de calcul du seuil d'arret (14 mai 02) ; Si l'image filtree d'echelle presente un maximum net on prend ; i_seuil(i_ech)=1, sinon on prend i_seuil(i_ech)=2. ; j_seuil = i_seuil(i_ech) ; D'avant le 14 mai 02. Sera ecras'e ; car il est prevu le 14 mai 02 de ; choisir i_seuil d'apres l'allure ; de l'image filtree d'echelle. ; OBSOLESCENCE FIN ; Calcul de "seuil_amp", quantite positive, valeur du max de la valeur ; absolue ou valeur rms de l'image residuelle sur la region netto- ; yee (champ total pour mode_clean = 1 2 ou 3, dom_c pour ; mode_clean = 4 ou 5) pour laquelle on arretera l'iteration ; dans CLEAN_ECHELLE. ; Ce critere d'arret est en competition avec d'autres criteres ; d'arret dans CLEAN_ECHELLE. ; Rappel : pour mode_clean=4 on soustrait les composantes clean ; interieure a dom_c, et l'image finale est l'image residuelle ; sale. Pour mode_clean=5 on nettoie une source intense de ses ; propres secondaires. Dans les deux cas il y a une source inten- ; se dans dom_c. Le max dans dom_c est en general le max absolu ; et il n'est pas besoin de specifier max(image(dom_c)). ; Rem : prendre le maximum sur dom_c equivaut a prendre le maximum ; absolu, puisque en pratique dom_c englobe une region particu- ; lieremt intense. On a mis dom_c pour etre coherent avec les ; definitions des maxima d'images residuelles manipulees dans ; CLEAN_ECHELLE. ; Amelioration (14 mai au 14 juin 02) : ; La dynamique sur une image nettoyee depend certainement de ; l'aspect de l'image. Il est intuitif que cette dynamique est ; plus importante quand l'image est dominee par une source com- ; pacte que quand c'est un magma complique. ; - dans le 1er cas on admet une fraction du max, mettons 1%, ; - dans le 2eme cas --------------------- de la valeur rms sur ; un domaine ou la source est importante (domaine plus petit ; que le soleil calme) mettons 5%. ; Dans le cas general on calcule la valeur rms de la source en ; excluant les sources compactes puissantes (ce que fair MAX_RMS) ; et on adopte comme seuil d'arret la plus grande de ces deux ; valeurs. ; Ca parait plus logique, pour les sources compliquees, que de ; comparer la moyenne rms du residu a celle de l'image initiale ; puisque la soustraction des composantes clean est essentielle- ; ment locale. ; ; En fait ce raisonnement intuitif ne vaut que pour des images ; non filtrees d'echelle. Une image non filtree peut consister en ; un fond de soleil calme bossele complexe plus ou moins domin'e ; par une ou plusieurs sources compactes localisees (beaucoup ; moins etendues que le SC). ; Dans les images filtrees d'echelle la base et les sources ; compactes ont des echelles moins separees (c'est d'ailleurs le ; but du filtrage d'echelle) : ; . pour i_ech=0, qui ne comporte que des tres bas harmoniques) ; l'image est tres large et on ne distingue pas entre base ; et source compacte. ; . pour i_ech=1 ou >1 il y a de forts secondaires et une source ; compacte est moins localisee que sur l'image non filtree. ; La procedure MAX_RMS (dont le but est de calculer les valeurs ; maximum et rms d'une image interferometrique filtree d'echelle) ; excise bien les sources compactes pour des images non filtrees ; et on obtient alors le max et une valeur rms typique du fond ; (source compacte exclue). ; Pour des images filtrees d'echelle MAX_RMS perd son efficacite ; Les essais du 14 juin 02 montrent que les valeurs max et rms ; sont essentiellement proportionnelles pour des sources consti- ; tuees d'un SC complexe (toujours la meme) et d'une source com- ; pacte dont on fait varier la brillance de 1 a 100 fois celle ; du SC. ; => pour des images filtrees d'echelle on ne peut prendre comme ; seuil d'arret qu'une fraction du max et non de la valeur ; moyenne rms (conclusion du 14 juin 02). ; Toutefois on prendra (automatiquement) un seuil 3 fois plus ; faible s'il existe une source compacte tres dominante. ; ; On prend toujours j_seuil=1 (depuis le 14 mai 02), cad on teste ; le max de abs(residu). ; ; Les valeurs de coeff_seuil choisies dans rh_dpnew servent a cal- ; culer les fractions du max et de la moyemme rms dont il est ; question plus haut. seuil_amp_max = coeff_seuil(1, i_ech) * max_im * freq/164 seuil_amp_rms = coeff_seuil(2, i_ech) * rms_im * freq/164 if (i_ech_max eq 0) then begin ; Il n'y a pas de filtrage d'echelle. Le seuil d'arret d'tera- ; tion peut etre defini par le maximum d'image ou la moyenne ; rms de l'image (excluant la contribution des sources com- ; pactes intenses), qui sont des quantites independantes. seuil_amp = seuil_amp_max > seuil_amp_rms if (seuil_amp_max gt seuil_amp_rms) then i_seuil_max_rms = 1 if (seuil_amp_max le seuil_amp_rms) then i_seuil_max_rms = 2 ; i_seuil_max_rms est transmis a CLEAN_ECHELLE pour impression ; de controle seulement. endif if (i_ech_max ge 1) then begin ; Les images sont filtrees d'echelle. Le maximum d'image et la ; moyenne rms fournies par MAX_RMS ne sont pas assez indepen- ; dantes et le seuil est defini par le maximum d'image seul. ; Le coeff_seuil defini dans RH_DPNEW est valable pour une ; source compacte tres dominante. Pour une source compacte ; plus faible on multiplie ce seuil par fac_mul compris entre ; 1 et 3. i_seuil_max_rms = 1 seuil_amp = fac_mul * seuil_amp_max ; fac_mul a ete calcul'e plus haut. endif ; Augmentation du seuil pour les echelles qui contiennent peu d'energie ; Modif du 25 sept 01 : quand l'image sale est lisse : ; . la valeur rms des images filtrees pour les petites echelles ; est faible, ; . le nbre des iterations est inutilement eleve. ; On tente de limiter ce nombre d'iterations. Le facteur ( )^0.3 ; reduit efficacement (mais pas trop...) quand la source est ; lisse. ; Augmenter l'exposant dans fac_seuil augmente le seuil et reduit ; le nombre d'iterations. Exposant usuel : 0.3 fac_seuil = (rms_im_ech(0) / rms_im_ech(i_ech))^0.3 seuil_amp = seuil_amp * fac_seuil j_seuil = 1 ; pour toutes les echelles apres le 15 mai 02. ; Controle de seuil_amp : voir apres l'obsolescence. ; OBSOLESCENCE DEBUT i_cadavre = 0 ; 1 pour ranimer le cadavre d'avant le 14 mai 02. if (i_cadavre eq 1) then begin max_im_ech = max(im_ech(dom_c), min=min_im_ech) ; cadavre. max_abs_im_ech = abs(max_im_ech) > abs(min_im_ech) ; ------- if (j_seuil eq 1) then begin ; seuil pr max image initiale. seuil_amp = coeff_seuil(1, i_ech) * max_abs_im_ech fac_seuil = 1 ; bidon endif if (j_seuil eq 2) then begin ; seuil pr rms image initiale mom_im_ech = moment (im_ech) ; La fonction "moment" calcule les 4 premiers moments de ; l'image filtree (non masquee depuis l'abandon des deloui- ; series). seuil_amp = coeff_seuil(2, i_ech) * $ sqrt(mom_im_ech(1) + mom_im_ech(0)^2) ; seuil_amp est la valeur quadratique moyenne, y compris ; dans le cas i_ech=0, quand le flux de l'image n'est pas ; nul. fac_seuil = (rms_im_ech(0) / rms_im_ech(i_ech))^0.3 ; augmenter l'exposant dans fac_seuil augmente le seuil et ; reduit le nombre d'iterations. Exposant usuel : 0.3 seuil_amp = seuil_amp * fac_seuil ; Modif du 25 sept 01 : ; Quand l'image sale est lisse : ; . la valeur rms des images filtrees pour les petites ; echelles est faible, ; . le nbre des iterations pour les petites echelles aug- ; mente beaucoup, de facon inutile. On tente ici de ; limiter ce nombre d'iterations. ; Remarque sur les param pour calculer seuil_amp : ; Les valeurs de coeff_seuil usuelles dans dpnew : ; coeff_seuil(1, *) = [ 5, 20, 20, 20, 20] * 0.001 ; coeff_seuil(2, *) = [10, 10, 20, 20, 20] * 0.010 ; donnent des nbres d'iterations voisins pour les images ; des differentes echelles avec : ; . une source complexe (ensemble de sources compactes) ; pour laquelle les rms_im_ech sont du meme ordre, ; . un exposant 1.0 (dans ce cas pratiquement inoperant ; puisqu'il s'applique a un nbre de l'ordre de 1). ; On admets donc qu'ils sont bien adaptes au cas de sour- ; ces compactes. Le facteur ( )^0.3 reduit efficace- ; ment (mais pas trop...) quand la source est lisse. endif ; fin du cas i_seuil = 2 endif ; fin du calcul anterieur au 14 mai 02. ; OBSOLESCENCE FIN ; Controle de "seuil_amp", avec visualisation de l'image filtree d'echelle icon = 0 if (icon eq 1) then begin ch_format ="(i2)" ch_ech = string(i_ech, format=ch_format) window, 0 shade_surf, objet_simule xyouts, 0.5, 0.97, /normal, alignment=0.5, 'Objet simule' window, 1 shade_surf, im xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Image initiale' window, 2 shade_surf, im_ech xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Image filtree d''echelle, filtrage ' + ch_ech print, ' ' print, 'RH_CLEAN calcul du seuil d''arret i_ech =', i_ech if (seuil_amp_rms ne 0) then begin ch_format ="(' seuils : max =', g8.2, ' rms =', g8.2, " +$ "' rapport =', f5.2, ' fac_seuil =', f5.2)" print, format=ch_format, seuil_amp_max, seuil_amp_rms, $ seuil_amp_max / seuil_amp_rms, fac_seuil endif if (seuil_amp_rms eq 0) then begin ch_format ="(' seuils : max =', g8.2, " +$ "' fac_seuil =', f5.2)" print, format=ch_format, seuil_amp_max, fac_seuil endif if (j_seuil eq 10) then begin ; anciennement 1 (rustine) print, 'seuil (defini pr max image) =' ch_format = "(' seuil =', f8.4, ' min =', f8.4," + $ "' max =', f8.4)" print, format=ch_format, seuil_amp, min_im_ech, max_im_ech endif if (j_seuil eq 20) then begin ; anciennement 2 (rustine) ch_format = "(' valeur rms =', f8.4, ' min =', " + $ "f8.4, ' max =', f8.4)" print, 'seuil (defini pr rms image) =', seuil_amp if ((i_ch_min eq 0) and (i_ch_max eq 0)) then begin print, 'Image sans filtrage d echelle :' print, format=ch_format, mom_im_ech(1)^0.5, $ min_im_ech, max_im_ech endif else begin ; filtrage d'echelle print, 'Image filtree d echelle :' print, format=ch_format, mom_im_ech(1)^0.5, $ min_im_ech, max_im_ech endelse endif stop endif ; fin de mise au point. ACCORD LE 19 AUG 99 AVEC LA VERSION JM. ; Calcul de "seuil_flux", quantite positive, valeur de la valeur absolue ; du flux de l'image residuelle pour laquelle on arretera l'ite- ; ration CLEAN dans CLEAN_ECHELLE. ; Ce critere d'arret est en competition avec d'autres criteres ; d'arret dans CLEAN_ECHELLE. seuil_flux = frac_flux * abs(total(im_ech)) ; Reduction du nbre de points sur im_ech, l_th_ech et i_fil (image filtre) ; et modification associee de dom_c (on tente ainsi de reduire ; le nombre des iterations dans CLEAN_ECHELLE) ; Rappel : dans FILTR_ECH les limites des domaines de filtrage sont ; des ellipses dont les axes (EW et NS) sont, pour i_ech de 0 a 4 ; (64, 50) (30, 30) (11, 16) (5, 8) (3, 4) if (i_ech eq 0) then reduc = 1 if (i_ech eq 1) then reduc = 1 if (i_ech eq 2) then reduc = 1 ;2 if (i_ech eq 3) then reduc = 1 ;4 if (i_ech eq 4) then reduc = 1 ;8 np_r = nx / reduc ; Rappel : quand congrid reduit le nbre de pts par un facteur en- ; tier n il prend un pt sur n, au contraire de rebin qui fait ; une moyenne avec les pts suivants. Il n'y a alors ni pb d'in- ; terpolation ni biais au bord droit du champ qui se trouve seu- ; lement un peu reduit. ; Par contre il y a interpolation quand la reduction du nbre de ; pts n'est pas entiere. im_ech_0 = im_ech & l_th_ech_0 = l_th_ech ; sauvegardes. i_fil_0 = i_fil ; ----------- if (reduc eq 1) then begin im_ech = im_ech_0 l_th_ech = l_th_ech_0 i_fil = i_fil_0 dom_c_ech = dom_c endif if (reduc gt 1) then begin ; ----------- im_ech = congrid(im_ech_0 , np_r, np_r) * reduc^2 l_th_ech = congrid(l_th_ech_0, np_r, np_r) * reduc^2 i_fil = congrid(i_fil_0 , np_r, np_r) vis_dom_c = intarr(nx, nx) ; pour visualiser dom_c. vis_dom_c(dom_c) = 1 vis_dom_c_ech = congrid(vis_dom_c, np_r, np_r) dom_c_ech = where(vis_dom_c_ech eq 1, nbre) ; Note du 4 jam 01 : ca fait sombrement suer de modifier dom_a. ; Donc on le passe sans reduction a CLEAN_ECHELLE. endif ; Appel a CLEAN_ECHELLE. ; Reduction du nombre d'iterations d'affilee pour affiner leur nbre total if (i_ech le 1) then jmax_ech = jmax if (i_ech eq 2) then jmax_ech = jmax / 2 if (i_ech eq 3) then jmax_ech = jmax / 4 if (i_ech eq 4) then jmax_ech = jmax / 8 if (jmax_ech lt 2) then jmax_ech = 2 ; precaution pour "moment". t0 = systime(1) ; temps en sec depuis le 1er jan 70. if (i_ech eq i_ech_min) then begin if (i_simul_soleil eq 0) then begin if (ch_pol eq 'non polar') then ch_format = $ "('RH_CLEAN : donnees reelles', f6.1, ' MHz non polar')" if (ch_pol eq 'polar') then ch_format = $ "('RH_CLEAN : donnees reelles', f6.1, ' MHz polar')" ; print, format=ch_format, freq endif if (i_simul_soleil gt 0) then begin if (ch_pol eq 'non polar') then ch_format = $ "('RH_CLEAN : donnees simulees', f6.1, ' MHz non polar')" if (ch_pol eq 'polar') then ch_format = $ "('RH_CLEAN : donnees simulees', f6.1, ' MHz polar')" ; print, format=ch_format, freq endif endif CLEAN_ECHELLE, $ ; entree: mode_clean, i_ech_min, i_ech_max, i_ech, $ im_ech, l_th_ech, i_fil, dom_c_ech, $ dom_a, $ itermax, jmax_ech, crit_arret, gain_iter, $ j_seuil, fac_crois_res, seuil_amp, $ i_seuil_max_rms, fac_seuil, seuil_flux, $ i_test_arret, $ ; i_seuil_max_rms et fac_seuil pour affichage. ; sortie: im_ech_s, im_ech_r ; Rappel 1 : im_ech, l_th_ech et i_fil sont ici definis sur le ; nbre de pts reduit. ; Rappel 2 : im_ech_s est "l'image synthetique", cad la somme ; des composantes clean habillees avec i_fil, sans addi- ; tion d'image residuelle. t1 = systime(1) ; temps en sec depuis le 1er jan 70. ch_1 = "(' RH_CLEAN, passage de CLEAN_ECHELLE : ', " ch_2 = "' sec')" ch_format = ch_1 + " f10.3," + ch_2 ; print, format=ch_format, t1 - t0 ; Controle : comparaison image sale et de l'image synthetique (somme des ; composantes clean habillees avec i_fil) filtr'ees d'echelle. icon = 0 if (icon eq 1) then begin print, 'RH_CLEAN : i_ech =', i_ech shade_surf, congrid(im_ech , 512, 512, cubic=-0.5) window, 1 & wset, 1 shade_surf, congrid(im_ech_r, 512, 512, cubic=-0.5) window, 2 & wset, 2 shade_surf, congrid(im_ech_s, 512, 512, cubic=-0.5) window, 3 & wset, 3 shade_surf, congrid(im_ech_s + im_ech_r, 512, 512, cubic=-0.5) print, 'image sale : min, max=', min(im_ech), max(im_ech) print, 'image synthetique: min, max=', min(im_ech_s), max(im_ech_s) stop endif ; Fin de controle ; Restauration du nbre de points initial (pour sommer les images syntheti- ; ques filtr'ees d'echelle) if (reduc eq 1) then begin ; rien a faire ! endif if (reduc gt 1) then begin im_ech_s = RH_INTER_FFT(im_ech_s, nx) im_ech_r = RH_INTER_FFT(im_ech_r, nx) i_fil = RH_INTER_FFT(i_fil , nx) endif ; Calcul de l'image propre finale pour un domaine d'echelle ; Suivant les valeurs de mode_clean : ; 1 (clean usuel), recherche des composantes clean sur toute l'ima- ; ge interferometrique. L'image finale est la somme de : ; . image synthetique filtree d'echelle, ; . image sale residuelle -------------. ; 2 idem mais on n'ajoute pas le residu a l'image synthetique. ; 3 recherche des composantes clean sur toute l'image. On ne resti- ; tue dans l'image synthetique que les composantes clean dont la ; position est exterieure a dom_c. L'image finale est la somme de ; . image synthetique filtree d'echelle (incomplete), ; . image sale residuelle -------------. ; On nettoie ainsi toute l'image et on supprine l'image synthe- ; tique d'un centre d'orage intense. ; 4 recherche des composantes clean sur le domaine dom_c pour les ; soustraire. L'image finale est l'image residuelle. ; Cela revient a debarasser l'image d'un centre d'orage localise ; et de ses secondaires, et le restant n'etant pas nettoy'e. ; 5 pour obtenir rapidement une image d'un orage intense interieur ; a dom_c, debarasse seulement de ses propres secondaires : on ; restreint la recherche des composantes a dom_c. Le reste n'est ; pas nettoy'e et ajout'e a l'image synthetique partielle. if (mode_clean eq 1) then begin im_p = im_p + im_ech_s + im_ech_r ; d'apres Thomson p 345 endif if (mode_clean eq 2) then begin ; image synthetique im_p = im_p + im_ech_s ; sans le residu. endif if (mode_clean eq 3) then begin im_p = im_p + im_ech_s + im_ech_r ; d'apres Thomson p 345 endif if (mode_clean eq 4) then begin ; on ne garde que le residu. im_p = im_p + im_ech_r endif if (mode_clean eq 5) then begin ; comme pour 1 ou 3. im_p = im_p + im_ech_s + im_ech_r endif im_r = im_r + im_ech_r endfor ; fin de boucle sur echelles (indice i_ech). ; sigma_filtre ; Impression des valeurs quadratiques moyennes de images filtrees d'chelle i_impression = 0 if (i_impression eq 1) then begin print, 'RH_CLEAN : valeurs quadratiques moyennes d''mages ' + $ 'filtrees d''echelle :' print, ' (determinant les facteurs mult ci-dessus)' ch_format = "(' image non filtree : ', f6.3)" print, format=ch_format, rms_im_nf for i_ech=i_ech_min, i_ech_max do begin ch_format = "(' image filtree ', i1, ' : ', f6.3)" print, format=ch_format, i_ech, rms_im_ech(i_ech) endfor print, ' ' endif ; Lissage de l'image finale (moyenne avec les 6 voisins immediats ou filtrage ; gaussien). if (i_lissage eq 1) then begin im_1 = shift(im_p, 1, 0) ; im a droite jusqu'au 29 nov 00. im_2 = shift(im_p, 1, 1) im_3 = shift(im_p, 0, 1) im_4 = shift(im_p, -1, 1) im_5 = shift(im_p, -1, 0) im_6 = shift(im_p, -1, -1) moy_v = (im_1 + im_2 + im_2 + im_3 + im_4 + im_5 + im_6) / 6 im_p = poids_c * im_p + (1 - poids_c) * moy_v ; rappel : usuel poids_c = 0.7 endif ; Calcul du filtre dans le plan (u, v) if (i_lissage eq 0) then begin ; jamais le cas fil_lis = replicate(1, nx, nx) endif if (i_lissage gt 1) then begin ; toujours le cas. Filtrage gaussien. ; Rappel de RH_MALC_IM_2D : usage de FFT directe au sens d'IDL. ; (drapeau -1) pour calculer l'image. vu = indgen(nx) - nx/2 v1 = replicate(1, nx) mu = vu # v1 ; matrice des abscisses dans le plan uv. mv = transpose(mu) if (i_lissage eq 2) then begin ; filtrage gaussien de revolution. dist_c = shift(dist(nx), nx/2, nx/2) ; distance au centre sig = larg_2(0) / sqrt(2 * alog(2)) arg = dist_c^2 / (2 * sig^2) fil_lis = exp(-arg) endif if (i_lissage eq 3) then begin ; filtrage gaussien elliptique. ; larg_2 est la 1/2 largeur a mi-puissance du filtre (en nbre ; d'harmoniques. sigu = larg_2(0) / sqrt(2 * alog(2)) sigv = larg_2(1) / sqrt(2 * alog(2)) arg = mu^2 / (2 * sigu^2) + mv^2 / (2 * sigv^2) fil_lis = exp(-arg) endif fil_lis_0 = fil_lis ; sauvegarde. ; Raccord progressif a zero du filtre en bord du domaine uv. ; raccord lineaire sur un bordure de largeur l_b i_saut = 0 ; 0 sauf pour essais. if (i_saut eq 1) then goto, a1 ; print, 'RH_CLEAN : adoucissement des bords du filtre de lissage final' i_raccord = 2 ; 1 pour raccord lineaire de l_b a np ; 2 pour multiplication par une cloche qui ; s'annule au bord de la couronne externe ; de filtrage. if (i_raccord eq 1) then begin l_b = 32 fac_bord = findgen(l_b) / (l_b - 1) ; de 0 a 1 . fac_bord_rev = reverse(fac_bord) ; de 1 a 0 . fac_x = replicate(1., np) fac_x( 0 : l_b - 1) = fac_bord fac_x(np - l_b : np- 1 ) = fac_bord_rev fac_y = fac_x fac_tour = fac_x # fac_y fil_lis = fil_lis * fac_tour endif if (i_raccord eq 2) then begin ; la couronne externe est elliptique. Par simplicite on defi- ; nit un rayon moyen de sa limite externe (ou elle s'annule) ; pour definir une cloche de revolution. if (i_ech_max eq 0) then begin ; C'est qu'on a pris i_filtr_ech = 0 dans rh_dpnew. ; On definit les axes de la limite du raccord de la meme ; facon que pour la couromnne la plus externe dans le cas ; du filtrage d'echelle. ax_e0 = [61.8, 52.0] endif ax_e0_m = (ax_e0(0) + ax_e0(1)) / 2 dist_c = shift(dist(nx), nx/2, nx/2) ; distance au centre ar = dist_c / ax_e0_m ar(where(ar ge 1)) = 1 cloche = (1 - ar^4)^2 ; shade_surf, cloche ; stop fil_lis = fil_lis * cloche endif a1: J_M = 0 ; Controle du filtre icon = 0 if (icon eq 1) then begin im_fil_lis = float(fft(shift(fil_lis, -nx/2, -nx/2), -1)) im_fil_lis = shift(im_fil_lis, nx/2, nx/2) window, 0 n1 = nx/2-64 & n2 = nx/2+63 ; reduction a 64 harm shade_surf, fil_lis(n1:n2, n1:n2) xyouts, 0.5, 0.87, /normal, alignment=0.5, $ 'filtre de lissage sur -+ 64 harm (raccorde a zero)' window, 1 n_z = 3 ; 1/2 largeur du zoom. ch_n_z = string(format="(i2)", n_z) n1 = nx/2-n_z & n2 = nx/2+n_z absc = - n_z + indgen(2*n_z + 1) plot, absc, im_fil_lis(n1:n2, nx/2) xyouts, 0.5, 0.87, /normal, alignment=0.5, $ 'coupe image du filtre de lissage sur -+' + ch_n_z + $ ' pixels' window, 2 shade_surf, fil_lis_0 ; - fil_lis xyouts, 0.5, 0.87, /normal, alignment=0.5, $ 'Filtre initial (non raccode a zero)' ; 'Filtre initial - filtre raccorde a zero' stop endif ; Filtrage (pour la clarte on manipule des TF quasi_reelles) tf_im_p = shift(fft(shift(im_p, -np/2, -np/2), 1), np/2, np/2) im_p = shift(fft(shift(tf_im_p * fil_lis, -np/2, -np/2), -1), $ np/2,np/2) im_p = float(im_p) if (i_simul_soleil gt 0) then begin tf_o = tf_objet_simule ; notation abregee. objet_simule_lisse = $ shift(fft(shift(tf_o * fil_lis, -np/2, -np/2), -1), np/2,np/2) ; pour controle dans RH_CALC_IMAGE_HELIO. objet_simule_lisse = float(objet_simule_lisse) endif else objet_simule_lisse = 0 ; bidon ; Comparaison objet simul'e et image finale liss'es par le meme filtre ; (doit se faire dans le if du filtrage gaussien) icon = 0 if (icon eq 1) then begin nz1 = 46 ; zoom de -nz1 a +nz1 autour de l'origine. i_s_t = 1 ; 1 pour shade_surf, 2 pour tcvsl. n1 = np/2 - nz1 & n2 = np/2 + nz1 ch_format = "('zoom de ', i3, ' a ', i3)" ch_zoom = string(format=ch_format, -nz1, nz1) window, 0 osl = congrid(objet_simule_lisse(n1:n2, n1:n2), 512,512,cubic=-0.5) if (i_s_t eq 1) then shade_surf, osl if (i_s_t eq 2) then tvscl , osl xyouts, 0.50, 0.97, /normal, alignment=0.5, $ 'Objet simule lisse par le filtre final, ' + ch_zoom window, 1 isl = congrid(im_p(n1:n2, n1:n2), 512, 512, cubic=-0.5) if (i_s_t eq 1) then shade_surf, isl if (i_s_t eq 2) then tvscl , isl xyouts, 0.50, 0.97, /normal, alignment=0.5, $ 'Image propre lissee par le filtre final, ' + ch_zoom stop endif endif ; fin du filtrage gaussien. ; Calcul de la TF du lobe synthetique ; C'est le produit de "sigma_filtre" (somme des filtres de filtrage ; d'echelle) par le filtre de lissage final. ; Il s'agit ici du lobe dans l'espace intyerferometrique, et il faudra ; le transformer en lobe heliographique dans RH_CALC_IMAGE_HELIO. tf_lobe_synt = sigma_filtre * fil_lis lobe_synt = float (shift (fft (shift ($ tf_lobe_synt, -nx/2, -nx/2), -1), nx/2, nx/2) ) ; Controle du lobe synthetique et de sa TF icon = 0 if (icon eq 1) then begin help, tf_lobe_synt, lobe_synt j_sh_tv = 2 ; trace tf_lobe_synt if (j_sh_tv eq 1) then begin window, 0 loadct, 0 shade_surf, abs(tf_lobe_synt) endif if (j_sh_tv eq 2) then begin window, 0, xsize=512, ysize=512 loadct, 13 tvscl, congrid(abs(tf_lobe_synt), 512, 512, cubic=-0.5) endif loadct, 0 xyouts, 0.5, 0.87, /normal, alignment=0.5, charsize=1.3, $ 'abs ( TF (lobe synthetique interferometrique) )' ; trace lobe_synt if (j_sh_tv eq 1) then begin window, 1 loadct, 0 shade_surf, lobe_synt endif if (j_sh_tv eq 2) then begin window, 1, xsize=512, ysize=512 loadct, 13 tvscl, congrid(lobe_synt, 512, 512, cubic=-0.5) endif loadct, 0 xyouts, 0.5, 0.87, /normal, alignment=0.5, charsize=1.3, $ 'lobe synthetique interferometrique' stop endif ; Controle du remplissage du plan uv et de l'effet du filtrage final ; (simulation seulement) icon = 0 if (icon eq 1) then begin a_ecrete = 0.95 ; niveau d'ecretage TF objet et image propre. nh = 25 ; zooming de -nh a + nh harmoniques ; TF objet simul'e i_s_f = 1 ; 1 pour multiplication par sigma_filtre if (i_s_f eq 0) then tf_o_f = tf_o if (i_s_f eq 1) then tf_o_f = tf_o * sigma_filtre dom_e = where(abs(tf_o_f) gt a_ecrete*abs(tf_o_f(np/2, np/2)), ipe) window, 0 tf_o_f_e = tf_o_f if (ipe gt 0) then tf_o_f_e(dom_e) = a_ecrete * abs(tf_o_f(np/2, np/2)) shade_surf, abs(tf_o_f_e(np/2-nh:np/2+nh, np/2-nh:np/2+nh)) xyouts, 0.5, 0.87, /normal, alignment=0.5, charsize=1.3, $ 'abs(TF(objet simule))' ; TF image sale tf_i_s = shift(fft(im, 1), np/2, np/2) window, 1 shade_surf, abs(tf_i_s(np/2-nh:np/2+nh, np/2-nh:np/2+nh)) xyouts, 0.5, 0.87, /normal, alignment=0.5, $ 'abs(TF(image sale))' window, 2 tvscl, congrid(abs(tf_i_s(np/2-nh:np/2+nh, np/2-nh:np/2+nh)), $ 512, 512, cubic=-0.5) xyouts, 0.5, 0.87, /normal, alignment=0.5, $ 'abs(TF(image sale))' stop ; TF image propre non lissee dom_e = where(abs(tf_im_p) gt a_ecrete*abs(tf_im_p(np/2, np/2)), ipe) tf_im_e = tf_im_p if (ipe gt 0) then tf_im_e(dom_e) = a_ecrete * abs(tf_im_p(np/2, np/2)) window, 0 shade_surf, abs(tf_im_e(np/2-nh:np/2+nh, np/2-nh:np/2+nh)) xyouts, 0.5, 0.87, /normal, alignment=0.5, $ 'abs(TF(image propre non lissee))' window, 1 tvscl, congrid(abs(tf_im_p(np/2-nh:np/2+nh, np/2-nh:np/2+nh)), $ 512, 512, cubic=-0.5) xyouts, 0.5, 0.87, /normal, alignment=0.5, $ 'abs(TF(image propre non lissee))' window, 2 n1 = np/2 - nh & n2 = np/2 + nh ymax = 1.05 * ( max(abs(tf_o (n1:n2, np/2 ))) > $ max(abs(tf_im_p(n1:n2, np/2 ))) > $ max(abs(tf_im_p(n1:n2, np/2 ))) ) absc = indgen(2*nh + 1) - nh plot, absc, abs(tf_o (n1:n2, np/2 )), xstyle=1, yrange=[0, ymax], $ ystyle=1 oplot, absc, abs(tf_im_p(n1:n2, np/2 )), linestyle=2 oplot, absc, abs(tf_im_p(n1:n2, np/2+1)), linestyle=1 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'TF coupes //u. objet ___, im propre v=0 - - et v=1...' stop ; TF (image propre - objet simul'e) uv_0 = shift(dist(nx), nx/2, nx/2) ; distance au centre du plan uv dom_ext = where(uv_0 gt 40, ip) diff = tf_im_p - tf_o_f ; difference avec la tf de l'objet tronquee if (ip gt 0) then diff(dom_ext) = 0 window, 0 shade_surf, abs(diff(np/2-nh:np/2+nh, np/2-nh:np/2+nh)) xyouts, 0.5, 0.87, /normal, alignment=0.5, $ 'abs(TF(image propre) - TF(objet simule)) pour (tf_im_p ne 0)' window, 1 tvscl, congrid(abs(diff(np/2-nh:np/2+nh, np/2-nh:np/2+nh)), $ 512, 512, cubic=-0.5) xyouts, 0.5, 0.87, /normal, alignment=0.5, $ 'abs(TF(image propre) - TF(objet simule)) pour (tf_im_p ne 0)' stop ; TF objet simule liss'e pour comparaison avec image lissee. tf_o_lis = tf_o * fil_lis if (i_s_f eq 0) then tf_o_lis_f = tf_o_lis if (i_s_f eq 1) then tf_o_lis_f = tf_o_lis * sigma_filtre dom_e = where( $ abs(tf_o_lis_f) gt a_ecrete*abs(tf_o_lis_f(np/2, np/2)), ipe) tf_o_lis_f_e = tf_o_lis_f if (ipe gt 0) then tf_o_lis_f_e(dom_e) = $ a_ecrete * abs(tf_o_lis_f(np/2, np/2)) window, 0 shade_surf, abs(tf_o_lis_f_e(np/2-nh:np/2+nh, np/2-nh:np/2+nh)) xyouts, 0.5, 0.87, /normal, alignment=0.5, $ 'abs(TF(objet simule lisse))' ; TF image propre lissee tf_im_lis = tf_im_p * fil_lis dom_e = where(abs(tf_im_lis) gt a_ecrete*abs(tf_im_lis(np/2, np/2)), $ ipe) window, 1 tf_im_lis_e = tf_im_lis if (ipe gt 0) then tf_im_lis_e(dom_e) = $ a_ecrete * abs(tf_im_lis(np/2, np/2)) shade_surf, abs(tf_im_lis_e(np/2-nh:np/2+nh, np/2-nh:np/2+nh)) xyouts, 0.5, 0.87, /normal, alignment=0.5, $ 'abs(TF(image propre lissee))' ; TF (image propre lissee) - TF (objet simul'e liss'e) diff = (tf_im_p - tf_o_f) * fil_lis ; difference avec la tf de ; l'objet tronquee. window, 2 shade_surf, abs(diff(np/2-nh:np/2+nh, np/2-nh:np/2+nh)) xyouts, 0.5, 0.87, /normal, alignment=0.5, $ 'abs(TF(image propre lissee) - TF(objet simule lisse))' ; window, 5 ; shade_surf, fil_lis ; xyouts, 0.5, 0.87, /normal, alignment=0.5, $ ; 'filtre de lissage' stop ; Arret le 18 avr 02 : A faire a la rentree: ; . comparer TF(image propre) avec TF(objet simul'e) ; . reflections a creuser : ; 1) Si on suppose qu'on retrouve les memes composantes clean ; a toutes les echelles, ma fonction de remplacement n'est pas ; vraiment une gaussienne puisque c'est la TF d'une gaussienne ; (le filtrage final) par la somme des filtres, laquelle ; est bornee a la bordure de la couronne externe, et en plus ; n'est pas constante sur le glacis exterieur de la couronne ; externe. A la difference de JM, dont la somme des filtres ; etait directement une gaussienne. Il n'est donc pas etonnant ; qu'il y ait des secondaires negatifs. ; L'idee au depart etait (?) de se passer de filtrage final et ; on cherchait a eviter les transitions brusques. ; Au 17 avril 02 la TF de l'image propre s'annule (aux arrondis ; pres) au bord de la couronne externe et est nulle au dela. ; Il est donc inutile de faire un raccordement a zero progres- ; sif du filtre final au bord de la couronne ou au bord du do- ; maine de Fourier (usage de i_raccord). Et d'un. ; En plus la TF de (fonction de remplacement * filtre final) ne ; peut etre une gaussienne a cause du glacis exterieur de ; la couronne externe. ; Il vaudrait peut etre mieux avoir des couronnes sans glacis, ; couvrant tout le domaine de Fourier (donc une zone externe ; vaste si on prends a_red_fil=0.5 par exemple) et avoir un ; filtrage final ne depassant pas trop le pave central. ; 2) des couronnes annulaires exterieures au pave et ne contenant ; pas les barres de E2 et NS45 seraient peu couvertes et donc ; peu contraintes. Il semble donc qu'il vaille mieux que la ; couronne qui ne contient plus du tout le pave DOIVE contenir ; les barres pour etre un peu couverte. La consequence est que ; l'avant derniere couronne doit contenir au moins le bord ; du pave pour etre un peu contrainte. Avec i_ech_max=3 ; le pave devra contenir le disque, le 1er anneau et la moitie ; du 2eme. Le domaine exterieur devra contenir la moitie du ; 2eme anneau et tout le 3eme, qui englobe les barres. ; Vu ainsi les couronnes devront etre plates, on ne pourra plus ; beaucoup jouer sur a_red_fil, mais seulement sur larg_2. ; 3) sur le filtre final : reflechir s'il doit etre quasi nul a ; +- 64. Sans doute puisqu'il ne doit pas trop depasser le ; pave. ; 4) question : faut-il que la couronned la plus externe contienne ; un peu du pave pour etre mieux contrainte ? ; 5) refaire tous les essais en jouant sur les couronnes et le ; filtrage final. endif ; Controle: image sale initiale, image propre finale lissee, image residuelle icon = 0 if (icon eq 1) then begin ; Image originale max_o = max(im, min=min_o) ch_format = "('RH_CLEAN image originale : max =', f5.3, " + $ "' min =', f6.3)" print, format=ch_format, max_o, min_o window, 0 tvscl, congrid(im , 512, 512, cubic=-0.5) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Image originale' ; Image propre lissee max_p = max(im_p, min=min_p) frac = 1 ; utilise pour visualiser dom_c. i_ecrete = 0 if (i_ecrete eq 0) then im_p_e = im_p if (i_ecrete eq 1) then begin frac = 0.02 dom_ecr = where(im_p gt frac*max_p, ipe) im_p_e = im_p if (ipe gt 0) then im_p_e(dom_ecr) = frac * max_p endif ch_format = "('RH_CLEAN image propre : max =', f5.3, " + $ "' min =', f6.3)" print, format=ch_format, max_p, min_p window, 1 ; Visualisation du contour du domaine dom_c (pour tvscl seulement) im_dom_c = fltarr(nx, nx) im_dom_c(dom_a) = 0.0 * frac * max_p ; usuel 0.5 . 0 pour shade_surf tvscl , congrid(im_p_e + im_dom_c, 512, 512, cubic=-0.5) ; shade_surf, congrid(im_p_e + im_dom_c, 512, 512, cubic=-0.5) if (i_ecrete eq 0) then xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Image propre' if (i_ecrete eq 1) then xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Image propre ecretee' ; Image residuelle max_r = max(im_r, min=min_r) ch_format = "('RH_CLEAN image residuelle : max =', f5.3, " + $ "' min =', f6.3)" print, format=ch_format, max_r, min_r window, 2 ; Visualisation du contour du domaine dom_c im_dom_c = fltarr(nx, nx) im_dom_c(dom_a) = 0.7 * max_r ; usuel 0.5 tvscl, congrid(im_r + im_dom_c, 512, 512, cubic=-0.5) xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Image residuelle' stop endif ; Controle: coupes 1D (initialement a la suite du controle precedent) icon = 0 if (icon eq 1) then begin ; Coupes 1D. Transient : max en (78, 58) iy = 64 ; 58 pour transient. izoom = 63 ; trace de -izoom a +izoom. absc = -izoom + indgen(2*izoom + 1) window, 0 frac = 0.99 ; fraction d'ecretage image initiale. amax = max(im, min=amin) ymax = frac * amax ymin = frac * amin plot, absc, im(np/2-izoom:np/2+izoom, iy), xrange=[-izoom, izoom], $ xstyle=1, yrange=[-ymax, ymax], ystyle=1 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Image originale : coupe EW' window, 1 frac = 0.99 ; fraction d'ecretage image propre. amax = max(im_p, min=amin) ; rappel :im_p est lissee. ymax = frac * amax ymin = frac * amin plot, absc, im_p(np/2-izoom:np/2+izoom, iy), xrange=[-izoom, izoom], $ xstyle=1, yrange=[-ymax, ymax], ystyle=1 oplot, absc, im_r(np/2-izoom:np/2+izoom, iy), linestyle=1 xyouts, 0.5, 0.97, /normal, alignment=0.5, $ 'Images propre ___ et residuelle ... : coupes EW' stop ; Definitions des TF sale et propre centrees et zoomees nhz = 12 ; nbre d'harm zoomes. n1 = nx/2 - nhz & n2 = nx/2 + nhz tf_s_c = shift(fft(im, 1), nx/2, ny/2) ; TF sale centree. tf_s_c = tf_s_c(nx/2 - nhz : nx/2 + nhz, nx/2 - nhz : nx/2 + nhz) a_s_c = abs(tf_s_c) tf_s_c = tf_s_c + 0.00001 p_s_c = 180/!pi * atan(imaginary(tf_s_c), float(tf_s_c)) tf_p_c = shift(fft(im_p, 1), nx/2, ny/2) ; TF propre centree. tf_p_c = tf_p_c(nx/2 - nhz : nx/2 + nhz, nx/2 - nhz : nx/2 + nhz) a_p_c = abs(tf_p_c) tf_p_c = tf_p_c + 0.00001 p_p_c = 180/!pi * atan(imaginary(tf_p_c), float(tf_p_c)) stop ; Traces TF sale et propre centrees et zoomees window, 0 & wset, 0 shade_surf, congrid(a_s_c, 512, 512, cubic=-0.5) xyouts, 0.5, 0.97, /normal, alignment=0.5, 'abs(TF sale)' window, 1 & wset, 1 shade_surf, congrid(p_s_c, 512, 512, cubic=-0.5) xyouts, 0.5, 0.97, /normal, alignment=0.5, 'phase(TF sale), deg' window, 2 & wset, 2 shade_surf, congrid(a_p_c, 512, 512, cubic=-0.5) xyouts, 0.5, 0.97, /normal, alignment=0.5, 'abs(TF propre)' window, 3 & wset, 3 shade_surf, congrid(p_p_c, 512, 512, cubic=-0.5) xyouts, 0.5, 0.97, /normal, alignment=0.5, 'phase(TF propre), deg' stop endif ; Fin de controle. ; Test final du clean dans le cas soleil simul'e. Comparaison entre : ; . objet reel soleil (simule par RH_SIMUL_HARM_N), ; . image interferometrique sale, directement produite par le RH, ; . ----------------------- propre lissee avec dom_c (en option). icon = 1 if ((i_simul_soleil gt 0) and (icon eq 1)) then begin soleil = objet_simule tf_p = tf_im_p ; pour ne pas alterer tf_im_p plus bas. tf_s = tf_im ; ------------------- tf_im -------- ; tf_p = shift(fft(shift(im_p, -nx/2, -nx/2), 1), nx/2, nx/2) ; tf_s = shift(fft(shift(im , -nx/2, -nx/2), 1), nx/2, nx/2) ; Zooming des TFs icici izoom = 1 nz = 8 ; domaine de zoom : -nz a +nz. NZ DOIT ETRE PAIR. if (izoom eq 1) then begin tf_p = tf_p (nx/2-nz : nx/2+nz, ny/2-nz : ny/2+nz) tf_s = tf_s (nx/2-nz : nx/2+nz, ny/2-nz : ny/2+nz) endif ; Visualisation de dom_c i_dom_c = 0 ; 1 pour visualiser dom_c sur image propre lissee. ; Incompatible avec izoom=1 ci-dessous. ; Zooming de l'objet et des images (utile seulement pour une seule source ; compacte simulee au centre du disque) izoom = 0 nz = 8 ; domaine de zoom : -nz a +nz. NZ DOIT ETRE PAIR. if (izoom eq 0) then begin sol_z = soleil im_s_z = im ; sale. im_p_z = im_p ; propre lissee. endif if (izoom eq 1) then begin sol_z = soleil (nx/2-nz : nx/2+nz, ny/2-nz : ny/2+nz) im_s_z = im (nx/2-nz : nx/2+nz, ny/2-nz : ny/2+nz) ; sale im_p_z = im_p (nx/2-nz : nx/2+nz, ny/2-nz : ny/2+nz) ; propre endif ; Definition des modules et passage a 512 pts a_tf_p = abs(tf_p) a_tf_p_512 = congrid(a_tf_p, 512, 512, cubic=-0.5) a_tf_s = abs(tf_s) a_tf_s_512 = congrid(a_tf_s , 512, 512, cubic=-0.5) soleil_512 = congrid(sol_z , 512, 512, cubic=-0.5) im_ini_512 = congrid(im_s_z , 512, 512, cubic=-0.5) im_p_512 = congrid(im_p_z , 512, 512, cubic=-0.5) ; Calcul et impression des flux flux_o = total(soleil) flux_i = total(im) flux_p = total(im_p) print, 'RH_CLEAN : comparaison totaux des images initiales et propres' ch_format = "(' total (soleil) =', f6.2)" print, format=ch_format, flux_o ch_format = "(' total (image initiale) =', f6.2)" print, format=ch_format, flux_i ch_format = "(' total (----- finale ) =', f6.2)" print, format=ch_format, flux_p ; Ecretage soleil initial i_ecret= 0 ; 0 ou 1 if (i_ecret eq 1) then begin frac = 0.50 max_im = max(abs(soleil_512), ip) dom = where(soleil_512 ge frac * max_im) soleil_512(dom) = frac * soleil_512(ip) endif ; Ecretage image initiale i_ecret= 0 if (i_ecret eq 1) then begin frac = 0.50 max_im = max(abs(im_ini_512), ip) dom = where(im_ini_512 ge frac * max_im) im_ini_512(dom) = frac * im_ini_512(ip) endif ; Ecretage image propre i_ecret= 0 if (i_ecret eq 1) then begin frac = 0.50 max_im = max(abs(im_p_512), ip) dom = where(im_p_512 ge frac * max_im) im_p_512(dom) = frac * im_p_512(ip) endif ; Ecriture soleil simule n_fich = 'sol_sim.gif' write_gif, n_fich, bytscl(soleil_512) ; Affichage soleil reel, image sale et image propre window, 0 ; shade_surf, soleil_512 loadct, 13 tvscl , soleil_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, charsize=1.3, $ 'Soleil reel simule' window, 1 loadct, 0 shade_surf, im_ini_512 ; tvscl , im_ini_512 if (i_pave_plein eq 0) then ch = ' (pave lacunaire)' if (i_pave_plein gt 0) then ch = ' (pave plein)' xyouts, 0.5, 0.97, /normal, alignment=0.5, charsize=1.3, $ 'Image RH sale initiale' + ch window, 2 ; visualisation du domaine dom_c sur l'image interferometrique propre toto = max(abs(im_p), ip) i_col = (ip mod npew) ; pour origine en bas a gauche. i_lig = (ip / npew) ; ---------------------------- ; I MOD J is equal to the remainder when I is divided by J. max_im = im_p(i_col, i_lig) im_dom_c = fltarr(npew, npew) im_dom_c(dom_a) = 0.3 * max_im ; usuel 0.8 1.5 ? im_dom_c_512 = congrid(im_dom_c, 512, 512, cubic = 0.5) loadct, 0 if (i_dom_c eq 0) then begin shade_surf, im_p_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, charsize=1.3, $ 'Image propre finale apres lissage' + ch endif if (i_dom_c eq 1) then begin shade_surf, im_p_512 + im_dom_c_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, charsize=1.3, $ 'Image propre finale apres lissage' + ch + ' avec dom_c' endif window, 3 loadct, 13 ; Incrustation de l'echelle de couleurs en bas a gauche npc = 201 ; longueur de l'echelle de couleur. hpc = 22 ; hauteur ----------------------- p_max = max(im_p_512, min=p_min) & da_p = p_max - p_min ramp_n = findgen(npc) / (npc - 1) ramp_p = p_min + 1.05 * da_p * ramp_n ; il manque le haut de l'echelle si on ne prend pas 1.05 band_p = ramp_p # replicate(1, hpc) band_p (*, hpc-5 : hpc-1) = 0 ; Incrustation d'une echelle de pourcentage dans la bande vmp = 1.05 * p_max a_m = 0.43 ; bleu tres clair avec rainbow. for km = 0, 10 do begin ; boucle sur les ticks. if (km eq 0) then begin ; elargiss a 2 pixels k1 = 20 * km k2 = 20 * km + 1 endif if ( (km gt 0) and (km lt 10) ) then begin ; larg 3 pixels k1 = 20 * km - 1 k2 = 20 * km + 1 endif if (km eq 10) then begin ; elargissement a 2 pixels. k1 = 20 * km - 1 k2 = 20 * km endif lm = 2 ; longueur des ticks ("lt" interdit). if ((km eq 0) or (km eq 5) or (km eq 10) ) then lm = 4 band_p (k1 : k2, hpc - 1 - 4 : hpc - 1 - 4 + lm) = a_m * vmp endfor im_p_512(0, 0) = band_p if (i_dom_c eq 0) then begin tvscl , im_p_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, charsize=1.3, $ 'Image propre finale apres lissage' + ch endif if (i_dom_c eq 1) then begin tvscl , im_p_512 + im_dom_c_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, charsize=1.3, $ 'Image propre finale apres lissage' + ch + ' avec dom_c' endif ; Ecriture image propre (existe aussi dans CALC_IMAGE_HELIO) ch_code = string(format='(i2)', fix(100*a_red_fil)) ; a_red_fil est le facteur d'homothetie pour reduire l'extension ; des zones de filtrage d'echelle, defini dans dpnew. ; n_fich = 'im_' + ch_code + '.gif' n_fich = 'im_tata.gif' write_gif, n_fich, bytscl(im_p_512) ; sur /home/mercier/cleanself/. ; ; ex : read_gif, 'im_68.gif', image, r, g, b pour a_red_fil=0.68. stop ; Trace des TF sale et propre (peu utile avec le remplissage du pave) window, 0 shade_surf, a_tf_s_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, charsize=1.3, $ 'abs(TF image sale)' window, 1 shade_surf, a_tf_p_512 xyouts, 0.5, 0.97, /normal, alignment=0.5, charsize=1.3, $ 'abs(TF image propre)' ; Comparaison des TF sale et propre (si la TF sale est nulle, on l'inter- ; pole la ou elle est nulle parallelement a l'axe x) ; L'interpolation peut ne pas etre tres bonne pres d l'origine si ; la TF observee y a un gros pic (soleil calme complexe). On peut ; alors avoir un ecart en % important, alors que la comparaison ; de l'image propre lissee et du soleil simul'e initial est ; acceptable. dom_nul = where(a_tf_s lt 0.0001*max(a_tf_s), ic) if (izoom eq 0) then nt = nx if (izoom eq 1) then nt = 2*nz + 1 diff_tf = tf_p - tf_s if (ic gt 0) then begin ; interpolation de diff_tf. iy = dom_nul / nt ix = dom_nul - iy * nt diff_tf(ix, iy) = 0.5 * (diff_tf(ix-1, iy) + diff_tf(ix+1, iy)) ; pas de debordement d'indice si nz est pair. endif a_diff = abs(diff_tf) window, 2 shade_surf, 100 * a_diff / max(abs(tf_p)) xyouts, 0.5, 0.97, /normal, alignment=0.5, charsize=1.3, $ 'abs(TF propre - TF sale) / max(abs(TF propre)), (%)' xyouts, 0.5, 0.94, /normal, alignment=0.5, charsize=1.3, $ '(avec interpolation la ou TF sale=0)' window, 3 shade_surf, 100 * a_diff / abs(tf_p) xyouts, 0.5, 0.97, /normal, alignment=0.5, charsize=1.3, $ 'abs(TF propre - TF sale) / abs(TF propre), (%)' xyouts, 0.5, 0.94, /normal, alignment=0.5, charsize=1.3, $ '(avec interpolation la ou TF sale=0)' print, 'RH_CLEAN : fin de procedure' stop endif ; Controle icon = 0 if (icon eq 1) then begin print, 'max(im) =', max(im), ' max(im_p) =', max(im_p) tvscl, congrid (im, 512, 512, cubic=-0.5) window, 1 wset, 1 tvscl, congrid (im_p, 512, 512, cubic=-0.5) stop endif ; Controle: comparaison im_p en sortie de RH_CLEAN icon = 0 if (icon eq 1) then begin window, 0 & wset, 0 shade_surf, congrid(im_p, 512, 512, cubic=-0.5) ; im_p_c = shift(im_p, n/2, n/2) window, 1 & wset, 1 ; shade_surf, congrid(im_p_c, 512, 512, cubic=-0.5) shade_surf, congrid(im_p , 512, 512, cubic=-0.5) stop window, 2 & wset, 2 tvscl, congrid(im_p, 512, 512, cubic=-0.5) im_p_c = shift(im_p, n/2, n/2) window, 3 & wset, 3 tvscl, congrid(im_p , 512, 512, cubic=-0.5) stop ; Conclusion: il semble que im_p soit "super-nettoyee", avec seulement ; les composantes CLEAN, sans ajout de l'image residuelle, et centree ; dans les coins du tableau !!! ; Donc on prend reg=0 dans RH_MALC_IM_2D et on saute de force la regu ; larisation du CLEAN dans REGCLEAN. Ainsi le CLEAN se termine avec ; RH_CLEAN. endif end ; fin de RH_CLEAN