;+ ;NAME: ; CONTIG_PIXELS ;PROJECT: ; Image manipulation, NASA/GSFC, June 2004 ; James McAteer ;PURPOSE: ; To provide a skeletal image of an AR boundary. ; Utilises combination of threshold, contour and image size. ; IMPORTANT! Need to be careful to account for projection effects ; and image calibration. ; ;USAGE: ; result=contig_pixels(image) ; result=contig_pixels(image,threshold=threshold,lim_wr=lim_wr,contig=contig ; ,/quiet,/nodisp) ; result=contig_pixels(image,threshold=threshold,lim_wr=lim_wr,contig=contig,/neg) ; ;INPUTS: ; IMAGE - the input image of an active region (2D array) ; ;KEYWORD INPUTS: ; THRESH - the level of the contours (FLOAT or INTEGER),in Gauss. Default: 50G ; LIM_WR - the maximum linear distance from the image maximum. CAREFUL: If this is small, the output skeletal ; image will have straight edges. Default: 1/2 smallest image dimension ; CONTIG - the minimum number of contiguous pixels to define a contour. Default: 20 pixels ; CENTXY - define center of region. Default: Position of image maximum. ; QUIET - no messages displayed. Default: Messages on ; NODISP - no images displayed. Deffult: Display on ; POS - only return positive polarity of active region. ; NEG - only return negative polarity of active region. ; ;OUTPUTS: ; RESULT - a skeletal image of the AR boundary. ; ;KNOWN BUGS: ; This is really more to do with CONTOUR, which we I can't hack directly. ; There's a problem with the /CLOSE keyword to CONTOUR, which means ; that it doesn't always close the contour. ; Also have to be careful for 'doughnut' shapes ; as the contig requirement will mean it will probably only catch the outside contour. ; ;HISTORY: ; v1.0 - Written as getstuff.pro by David Williams & James McAteer, QUB (30/10/01) after ; much blood, sweat & tears. Originally a quiet Sun program for solar network ; v1.1 - JMA added IMAGEMAX keyword for RATIOLC usability (31/10/01) ; v1.2 - DRW added QUIET keyword to shut the thing up during looped ; use (2/11/01). ; v1.3 - JMA added ability to return value of -1 if the contours can't be ; drawn, OR if the threshold is too large. Stops the programme from ; bombing out (5/11/01). ; v1.4 - JMA added in loop to return valur of -1 and explain why: when the threshold ; is too large the previous addition (v1.3) was confusing. ; Also added added in a loop which corrects ; for the fact that every value of y between mny and mxy may not be taken ; v2.0 - JMA (NASA/GSFC) adapted program into new format, contig_pixels.pro ; v2.1 - JMA added in ability to include all of one contour, as long as one value of ; it is inside the lim_wr. Also added in a default return of a geometric ; shape (rectangle) if the threshold is set too high. ; v2.2 - JMA added in choose_params. Allow for program to overide contig & lim_wr ; keywords as fraction of largest contour. ; v2.3 - JMA added to choose_params. the smallest contour is defined as wher 3 contours ; of the same size occur. ; v2.4 - JMA hardcoded in a maximum lim_wr of 50 pixels. ;- FUNCTION contig_pixels,im,thresh=thresh,lim_wr=lim_wr,contig=contig,$ centxy=centxy,quiet=quiet,nodisp=nodisp,pos=pos,neg=neg,$ choose_params=choose_params, choose_params2=choose_params2 ini_im=im IF KEYWORD_SET(neg) THEN im(where(im GE 0)) = 0 IF KEYWORD_SET(pos) THEN im(where(im LE 0)) = 0 im=abs(im) simage=SIZE(im) IF KEYWORD_SET(choose_params) THEN BEGIN CONTOUR,im,levels=thresh,path_xy=xy,path_info=info,$ /path_data_coords,/closed IF N_ELEMENTS(info) EQ 0 THEN GOTO, no_contours ;account for no contours at all largest_contour=max(info.n) IF NOT KEYWORD_SET(quiet) THEN print,'Largest contours:',arr2str(largest_contour,/trim),' pixels' IF largest_contour LT 50 THEN GOTO, no_contours ;keep a 20 pixel threshold (50/e=20) sorted=reverse(sort(info.n)) sorted2=sorted for i=0,N_ELEMENTS(sorted)-1 do sorted2(i)=(info.n)(sorted(i)) for i=0,N_ELEMENTS(sorted2)-2 do sorted(i)=sorted2(i)-sorted2(i+1) for i=0,N_ELEMENTS(sorted)-2 do $ if sorted(i) eq 0 AND sorted(i+1) eq 0 then begin sz=sorted2(i) goto, found_size endif found_size: lim_wr = sz contig = sz ;lim_wr=50 < FIX((0.2)*largest_contour) ;1/5 fall off from largest contour size ;contig=FIX((0.2)*largest_contour) ;1/5 fall off from largest contour size GOTO, ch_params ENDIF IF KEYWORD_SET(choose_params2) THEN BEGIN CONTOUR,im,levels=thresh,path_xy=xy,path_info=info,$ /path_data_coords,/closed ;IF N_ELEMENTS(info) EQ 0 THEN GOTO, no_contours ;account for no contours at all IF N_ELEMENTS(info) LT 10 THEN GOTO, no_contours ;if there only a few contours, data is probably screwed lim_wr = FIX( (total(info(0:9).n))/50.) IF lim_wr LT 20 THEN lim_wr=20;miniumum of 20 pixels IF lim_wr GT 50 THEN lim_wr=50;maximum of 50 pixels ;contig = 2*lim_wr contig = lim_wr ;print,lim_wr,contig largest_contour=max(info.n) IF largest_contour LT 20 THEN GOTO, no_contours ;keep a 20 pixel threshold (50/e=20) GOTO, ch_params ENDIF ;Have a default wander limit of 1/2 smallest dimension IF KEYWORD_SET(lim_wr) THEN lim_wr=lim_wr ELSE lim_wr=((simage)(1) < (simage)(2))/2. -1 ;Have a default contiguous pixel number of 20 pixels IF KEYWORD_SET(contig) THEN contig=contig ELSE contig=20. ch_params: ;Have a default threshold of 50G IF KEYWORD_SET(thresh) THEN thresh=abs(thresh) ELSE thresh=50. IF NOT KEYWORD_SET(quiet) THEN BEGIN print,'Wander limit: ',arr2str(lim_wr,/trim),' pixels' print,'Minimum contiguous: ',arr2str(contig,/trim),' pixels' print,'Threshold Field: ',arr2str(thresh,/trim),' G' ENDIF ;Have a default image center of image maximum IF KEYWORD_SET(centxy) THEN BEGIN x=centxy(0) y=centxy(1) IF NOT KEYWORD_SET(quiet) THEN print,'Setting AR center: ',arr2str(FIX(x),/trim),', ',arr2str(FIX(y),/trim) ENDIF ELSE BEGIN maxim = WHERE(im eq MAX(im)) printxy,im,maxim,x,y IF NOT KEYWORD_SET(quiet) THEN print,'Image max: ',arr2str(FIX(MAX(im)),/trim),$ ' G occurs at ',arr2str(FIX(x),/trim),', ',arr2str(FIX(y),/trim) ENDELSE ;define region of interest as position of image max +/- wander limit, or edges of image x0=0 > FIX(x-(2*lim_wr)) y0=0 > FIX(y-(lim_wr)) x1=simage(1)-1 < FIX(x+(2*lim_wr)) y1=simage(2)-1 < FIX(y+lim_wr) ; If /NODISP is set, then don't display IF NOT KEYWORD_SET(nodisp) THEN BEGIN ;wset,2 loadct,0,/silent mult,2,1 plot_image,ini_im,tit='Original Image',min=-500,max=500,/nosquare loadct,3,/silent contour,im,levels=thresh,/over,/noerase,c_color=128 draw_boxcorn,x0,y0,x1,y1,/data verline,x horline,y END ;check to see if the threshold is set too high largest_value=( max( im( X0:X1,y0:y1)))(0) IF (thresh GT largest_value) or (n_elements(x) ne 1) THEN BEGIN ;print,'WARNING! Highest pixel value of ',arr2str(FIX(largest_value),/trim)$ ; ,'G, and threshold of ',arr2str(FIX(thresh),/trim), 'G .NO PIXELS CONTOURED' no_contours: ;print,'WARNING! NO PIXELS CONTOURED' null_im=im null_im(*,*)=0 ;sz=size(null_im,/dim) ;null_im( fix(sz(0)/4.):3*fix(sz(0)/4.), fix(sz(1)/4.):3*fix(sz(1)/4.))=1000 RETURN,null_im ENDIF ; Work out the path of the contours around pixels with values greater than THRESH ; Then store the details of these contours in a structure called INFO CONTOUR,im,levels=thresh,path_xy=xy,path_info=info,$ /path_data_coords,/closed perim=fltarr(simage(1),simage(2)) ;work through each contour, saving it to an perimeter image if it consists of enough pixels ;AND has at least one pixel inside the lim_wr region FOR i=0,(size(info))(1)-1 DO BEGIN ; for each contour.... sz_contour=info(i).n ;...calculate size of contour... start=info(i).offset ;...calculate starting position... FOR pixel_index=start,start+sz_contour-1 DO BEGIN ;...work through each pixel... ;...if any pixel appears inside the lim_wr region... IF xy(0,pixel_index) GT X0 AND xy(0,pixel_index) LT X1 AND xy(1,pixel_index) GT Y0 $ AND xy(1,pixel_index) LT Y1 THEN BEGIN ;...then save that contour to a binary image... IF sz_contour GT contig THEN perim(xy(0,start:start+sz_contour-1),xy(1,start:start+sz_contour-1))=1 GOTO,next_contour ;...go and find the next contour ENDIF next_contour: ENDFOR ENDFOR IF NOT KEYWORD_SET(nodisp) THEN BEGIN loadct,0,/silent ;wset,2 plot_image,perim,tit='Skeletal Image',/nosquare draw_boxcorn,x0,y0,x1,y1,/data verline,x horline,y mult,1,1 ENDIF ;restore original data and window preferences. im=ini_im mult,1,1 RETURN,perim END