PRO bovelet,image,seuil,pas,n_seuil,fenetre,image_lab ;,choix ; ; le 29 Avril 2003 ; programme fait par Thierry Roudier ; ; OBJET: reproduire la segmentation du papier de ; Bovelet et Wiehr Solar Physics 2001,201,13 ; ; Seuil = premier seuil a donner (First threshold to give) ; pas = pas de variation du seuil (step for between thresholds) ; n_seuil= nombre de seuils desires ( number of threshold) ; fenetre=fenetre de lissage =resolution souhaitee (smoothing window) ; choix = intensite ou vitesse (choice between intensity or velocity) ;pour notre exemple faire : ( example) ;bovelet,0.135,0.005,3,2,1 ; pour l'intensite ; ;bovelet,0.07,0.01,3,2,2 ; pour la vitesse ; ;OUTPUT: labelled image IMAGE_LAB (in a save file 'image_lab') ;------------------------------------------------------------- ;if (choix eq 1) then begin ; image=readfits('int.fits') ;endif else begin ; image=readfits('vit_z.fits') ;endelse dim=size(image) tvscl,image,dim(1)+1,0 ;premier seuil ;seuil=0.135 ; pour l'image en intensite image_bin=image*0 onrecommence: if (fenetre gt 1) then image=smoothe(image,fenetre) w1=where(image gt seuil,countseuil) if (countseuil le 0) then begin print,' descendez le seuil, seuil=' read,seuil goto,onrecommence endif print,' ======================> Premier seuil =',seuil image_bin(w1)=1.0 image_lab = LABEL_REGION(image_bin) image_ref=image_lab for nn=0,n_seuil-1 do begin seuil2=seuil-pas*(nn+1) print,' ======================> Nouveau seuil =',seuil2 print,' numero maximun=', max(image_lab) image_bin=image*0. w2=where(image gt seuil2) image_bin(w2)=10000. image_bin(w1)=image_lab(w1) image_bin_ref=image_bin nombre=0 While (nombre lt 100) do begin ; Donne la population et les membres de chaque region: h = HISTOGRAM(image_lab, REVERSE_INDICES=r) ; Pour chaque region FOR i=1, N_ELEMENTS(h)-1 DO BEGIN ; PRINT, 'Region ', i, ' Population = ', h[i] ; on recherche les coordonnees pour savoir la valeur wi=where(image_bin eq i,counti) x = wi mod dim[1] y = wi / dim[1] for j=0,counti-1 do begin x_ori= x(j) & y_ori= y(j) ;pour les x positifs if(x_ori+3 le dim(1)-1 and y_ori-3 ge 0 and y_ori+3 le dim(2)-1 ) then begin if( image_bin(x_ori+1,y_ori) eq i) then goto,suivant1 if (image_bin(x_ori+1,y_ori) eq 10000.) then begin if (image_bin(x_ori+3,y_ori) lt 10000. and image_bin(x_ori+3,y_ori) gt 0 and image_bin(x_ori+3,y_ori) ne i) then begin goto,suivant1 endif else begin if (image_bin(x_ori+1,y_ori+3) lt 10000. and image_bin(x_ori+1,y_ori+3) gt 0 and image_bin(x_ori+1,y_ori+3) ne i) then begin goto,suivant1 endif if (image_bin(x_ori+1,y_ori-3) lt 10000. and image_bin(x_ori+1,y_ori-3) gt 0 and image_bin(x_ori+1,y_ori-3) ne i) then begin goto,suivant1 endif image_lab(x_ori+1,y_ori)=i endelse endif endif suivant1: ;pour les x negatifs if(x_ori-3 ge 0 and y_ori-3 ge 0 and y_ori+3 le dim(2)-1 ) then begin if( image_bin(x_ori-1,y_ori) eq i) then goto,suivant2 if (image_bin(x_ori-1,y_ori) eq 10000.)then begin if (image_bin(x_ori-3,y_ori) lt 10000. and image_bin(x_ori-3,y_ori) gt 0 and image_bin(x_ori-3,y_ori) ne i) then begin goto,suivant2 endif else begin if (image_bin(x_ori-1,y_ori+3) lt 10000. and image_bin(x_ori-1,y_ori+3) gt 0 and image_bin(x_ori-1,y_ori+3) ne i) then begin goto,suivant2 endif if (image_bin(x_ori-1,y_ori-3) lt 10000. and image_bin(x_ori-1,y_ori-3) gt 0 and image_bin(x_ori-1,y_ori-3) ne i) then begin goto,suivant2 endif image_lab(x_ori-1,y_ori)=i endelse endif endif suivant2: ;pour les y positifs if(y_ori+3 le dim(2)-1 and x_ori-3 ge 0 and x_ori+3 le dim(1)-1 ) then begin if( image_bin(x_ori,y_ori+1) eq i) then goto,suivant3 if (image_bin(x_ori,y_ori+1) eq 10000.) then begin if (image_bin(x_ori,y_ori+3) lt 10000. and image_bin(x_ori,y_ori+3) gt 0 and image_bin(x_ori,y_ori+3) ne i) then begin goto,suivant3 endif else begin if (image_bin(x_ori+3,y_ori+1) lt 10000. and image_bin(x_ori+3,y_ori+1) gt 0 and image_bin(x_ori+3,y_ori+1) ne i) then begin goto,suivant3 endif if (image_bin(x_ori-3,y_ori+1) lt 10000. and image_bin(x_ori-3,y_ori+1) gt 0 and image_bin(x_ori-3,y_ori+1) ne i) then begin goto,suivant3 endif image_lab(x_ori,y_ori+1)=i endelse endif endif suivant3: ;pour les y negatifs if(y_ori-3 ge 0 and x_ori-3 ge 0 and x_ori+3 le dim(1)-1) then begin if( image_bin(x_ori,y_ori-1) eq i) then goto,suivant4 if (image_bin(x_ori,y_ori-1) eq 10000.) then begin if (image_bin(x_ori,y_ori-3) lt 10000. and image_bin(x_ori,y_ori-3) gt 0 and image_bin(x_ori,y_ori-3) ne i) then begin goto,suivant4 endif else begin if (image_bin(x_ori+3,y_ori-1) lt 10000. and image_bin(x_ori+3,y_ori-1) gt 0 and image_bin(x_ori+3,y_ori-1) ne i) then begin goto,suivant4 endif if (image_bin(x_ori-3,y_ori-1) lt 10000. and image_bin(x_ori-3,y_ori-1) gt 0 and image_bin(x_ori-3,y_ori-1) ne i) then begin goto,suivant4 endif image_lab(x_ori,y_ori-1)=i endelse endif endif suivant4: endfor ENDFOR win=where(image_lab gt 0) image_bin =image_bin_ref image_bin(win)= image_lab(win) tvscl, image_lab ; pour l'arret des iterations dans le while if (nombre gt 0.) then begin sum=0. FOR i=1, N_ELEMENTS(h)-1 DO sum= h[i]- h_ref[i] +sum if (sum eq 0.) then begin print, 'nombre d''iterations necessaires =',nombre goto,finito endif endif h_ref=h nombre=nombre+1 endwhile ; finito: ; on detecte les nouveaux granules a ce seuil. wfin=where(image_bin_ref gt 0) image_dessus=image_bin_ref*0. image_dessus(wfin)=1. image_dessus_lab = LABEL_REGION(image_dessus) h_des = HISTOGRAM(image_dessus_lab, REVERSE_INDICES=r) FOR i=1, N_ELEMENTS(h_des)-1 DO BEGIN ; PRINT, 'Region ', i, ' Population = ', h[i] ; on recherche les coordonnees pour savoir la valeur w_des=where(image_dessus_lab eq i,count_des) if ( total(image_lab(w_des)) eq 0.) then image_lab(w_des)=5000.+ i ENDFOR ; on relabellise pour avoir des numeros continus ; possibilite d'ameliorer les segmentions ; on effectue une erosion pour bien separer les granules. ; utile pour les images du Pic de la sequence de 3h. ;fac=1 ; taille de l'element structurant pour l'ouverture ; ;S = REPLICATE(1, fac, fac, 1) ;creation de l'element structurant ;image_lab= ERODE(image_lab, S) ; operateur d'ouverture. hsortie = HISTOGRAM(image_lab, REVERSE_INDICES=r) ; Pour chaque region numer =1. & image_lab_temp=image_lab*0. FOR mm=1, N_ELEMENTS(hsortie)-1 DO BEGIN if ( hsortie(mm) gt 0.) then begin ; print,' numero detecte =',mm wsor=where( image_lab eq mm,countmm) if (countmm le 0) then goto,suivant image_lab_temp(wsor)=numer numer=numer+1 endif suivant: ENDFOR image_lab=image_lab_temp endfor print,' numero maximun final =', max(image_lab) ;save,image_lab,file='image_lab' end