;+ ;Project : SOHO-LASCO/EIT ; ;Name : cr_rem_array ; ;Purpose : Remove cosmic rays from LASCO/EIT images ; ;Explanation : This procedure removes cosmic rays from an array ; previously defined in the idl session by first ; using statistical techniques to identify regions ; containing cosmic rays (or stars) and then replaces ; the pixel values in that region with new values ; determined by taking the average of column and row ; fits across the region based on neighboring points. ; Statistical noise based upon the quality of the ; fit to the neighboring data is added to the fit ; to give realistic local variation for aesthetic ; purposes. ; ; Syntax : result=cr_rem_array(image,telescope, ; BASE_IMAGE=base_image, ; BLOCK=block, ; GROUP_SIZE=groupsize, ; N_SIGMA=n_sigma, ; CX=cx,CY=cy, ; RMIN=rmin,RMAX=rmax, ; RFILE=rfile, ; USE_ROI=use_roi, ; ORDER=order) ; ; Inputs : image : The array name of a previously read-in image ; telescope : An integer denoting the telescope used for the ; image. C1=1, C2=2, C3=3, EIT=4 ; ; Keywords : BASE_IMAGE : Image to subtract from the input image array if ; differencing is desired ; BLOCK : Size in pixels of the larger neighborhoods used ; to determine the statistical characteristics ; which determine which smaller neighborhoods have ; possible cosmic rays. Default is 32. ; GROUP_SIZE : Size in pixels of the smaller neighborhoods ; which are examined for cosmic rays. ; Default is 8. ; N_SIGMA : Statistical requirement used to determine ; candidate smaller neighborhoods withing a larger ; neighborhood. The maximum value of a given ; smaller neighborhood must deviate ; from the median value of the means of the ; smaller neighborhoods comprising a larger ; neighborhood in order to be considered a ; candidate. ; ORDER : The order of the polynomial fit to use to model ; the data in the smaller neighborhoods being ; corrected. Default is 2. Errors may occur for ; values higher than 3. For smaller neighborhoods ; abutting the occulter or abutting telemetry ; dropout regions, a linear fit is used. ; ; *** The program is designed to handle the occulter and vignetted outer edge ; of the circular FOV of the instruments. The high intensity gradients ; at these edges leads to poor cosmic ray correction when candidate ; pixel groups are identified there. The following keywords should be ; included when using this feature. ; CX,CY : Coordinates of the center of the occulter in ; the image. Only required if the occulter region ; is to be ignored in the removal process. This ; is done in processing LASCO images to avoid ; problems arising from having large regions ; with values drastically different ; from neighboring valid data. In EIT images it is ; done to avoid processing disk data because some ; bright features on the disk are ; indistinguishable from cosmic rays using these ; statistical methods. ; RMIN,RMAX : The minimum and maximum radial values to ; consider for the processing. The minimum value ; should be slightly larger than the radius of ; the occulter in the image. The maximum should ; be slightly smaller than the radius of the ; unvignetted field of view. IF RMIN is defined ; but RMAX is not, the entire field of view ; outside RMIN will be corrected. ; RFILE : The name of a file containing the radial ; distance of each pixel in the image from ; the center coordinates CX, CY. If the file ; does not exist, it will be generated. If the ; file exists, but the size of the image or the ; center coordinates have changed, the file is ; regenerated. ; USE_ROI : Sometimes one might want to exclude a region from ; cosmic ray removal processing, such as when a bright ; star or a comet appears in the field of view. To ; do this, set this keyword and then use the mouse ; buttons to set the ROI (=region of interest). The left ; mouse button draws line segments, the middle button ; removes them if necessary, and the right button ; closes the region. See the description of DEFROI ; in the IDL manual. ; ; Calls : generate_r_matrix ; ; Restrictions : ; ; Side effects : ; ; Category : Image processing ; ; Prev. Hist. : None ; ; Written : Norm Moulton, NRL, Sept. 1996 ; ; Modified : ; ; Version : ;- ;______________________________________________________________________________ ; FUNCTION cr_rem_array,on_band,telescope,$ BASE_IMAGE=base_image,$ BLOCK=block,$ GROUP_SIZE=group_size,$ N_SIGMA=n_sigma,$ CX=cx,CY=cy,$ RMIN=rmin,RMAX=rmax,$ RFILE=rfile,$ USE_ROI=use_roi,$ ORDER=order common processed_image,new_image common detection_array,good_bad_invalid_array common fixing_images,row_image,col_image ;###HANDLE KEYWORDS### if not(keyword_set(ORDER)) then fit_order=2 $ else fit_order=order if not(keyword_set(RFILE)) then r_filename='default_rfile.dat' $ else r_filename=rfile if not(keyword_set(cx)) then ctrx=0 else ctrx=cx if not(keyword_set(cy)) then ctry=0 else ctry=cy if not(keyword_set(rmin)) then minr=0 else minr=rmin if not(keyword_set(rmax)) then maxr=0 else maxr=rmax if not(keyword_set(n_sigma)) then stat_factor=6 $ else stat_factor=n_sigma if not(keyword_set(group_size)) then nbrhd_size=8 $ else nbrhd_size=group_size if not(keyword_set(base_image)) then off_band=[0] $ else off_band=base_image if not(keyword_set(block)) then blocksize=32 $ else blocksize=block ;###Create Image Arrays for Program### if (n_elements(off_band) gt 1) then image=on_band-off_band $ else image=on_band new_image=image ;###Get or Generate Appropriate R Array### s=size(image) sx=s(1) sy=s(2) if ((minr gt 0) OR $ (maxr gt 0)) then r=generate_r_matrix(r_filename,sx,sy,ctrx,ctry) if ((minr gt 0) AND (maxr eq 0)) then $ MAX R UNDEFINED maxr=max(r) ;###Replace Bad Columns with Neighboring columns### if not(keyword_set(telescope)) then begin print,'Not fixing bad rows or columns for this image...' endif $ else $ if (telescope eq 1) then begin ;C1, cols 0, 1023 are bad print,'C1 image' image(0,*)=image(1,*) image(sx-1,*)=image(sx-2,*) new_image(0,*)=image(1,*) new_image(sx-1,*)=image(sx-2,*) endif $ else $ if (telescope eq 4) then begin ; C2 print,'C2 image' endif $ else $ if (telescope eq 3) then begin ; C3, rows 0, 1022,1023, col 0 print,'C3 image' image(0,*)=image(1,*) new_image(0,*)=image(1,*) image(*,0)=image(*,1) new_image(*,0)=image(*,1) image(*,sy-1)=image(*,sy-3) new_image(*,sy-1)=image(*,sy-3) image(*,sy-2)=image(*,sy-3) new_image(*,sy-2)=image(*,sy-3) endif $ else $ if (telescope eq 4) then begin ; EIT, rows 0,1 are bad print,'EIT image' image(*,0)=image(*,2) image(*,1)=image(*,2) new_image(*,0)=image(*,2) new_image(*,1)=image(*,2) endif $ else begin print,'INVALID TELESCOPE DESIGNATION' print,'Not fixing bad rows or columns for this image...' endelse tvscl,new_image ;###HANDLE USE_ROI### ;Define the patch and make a copy, then paste it back over after fixing is done if keyword_set(USE_ROI) then begin roi=defroi(sx,sy) patch=new_image(roi) on_band(roi)=0 ; This forces these pixels to be labelled invalid for fitting endif ;###Set Up Arrays for Processing Image### n_subsets_x=blocksize/nbrhd_size n_subsets_y=n_subsets_x sub_mean=fltarr(n_subsets_x,n_subsets_y) sub_sdev=fltarr(n_subsets_x,n_subsets_y) sub_max=intarr(n_subsets_x,n_subsets_y) sub_min=intarr(n_subsets_x,n_subsets_y) img_min_x_pix=intarr(n_subsets_x,n_subsets_y) img_min_y_pix=intarr(n_subsets_x,n_subsets_y) img_max_x_pix=intarr(n_subsets_x,n_subsets_y) img_max_y_pix=intarr(n_subsets_x,n_subsets_y) n_subsets_total_x=sx/nbrhd_size n_subsets_total_y=sy/nbrhd_size subset_good=intarr(n_subsets_total_x,n_subsets_total_y) subset_valid=intarr(n_subsets_total_x,n_subsets_total_y) subset_good(*,*)=1 subset_valid(*,*)=1 good_bad_invalid_array=intarr(sx,sy) good_bad_invalid_array(*,*)=0 nxblocks=s(1)/blocksize nyblocks=s(2)/blocksize ;++++++++++++++++++++++++++++++++ ; BEGIN MAIN PART OF PROGRAM ;++++++++++++++++++++++++++++++++ print,'Doing neighborhood-sampling...' print,'Finding bad and invalid subneighborhoods...' for in=0,nxblocks-1 do begin ;NEIGHBORHOOD X COUNTER for jn=0,nyblocks-1 do begin ;NEIGHBORHOOD Y COUNTER ;###SET THE NEIGHBORHOOD### neighborhood=image(in*blocksize:(in+1)*blocksize-1,$ jn*blocksize:(jn+1)*blocksize-1) ;###GET SUBSET STATISTICS### for i=0,n_subsets_x-1 do begin for j=0,n_subsets_y-1 do begin subset=neighborhood(i*nbrhd_size:(i+1)*$ nbrhd_size-1,$ j*nbrhd_size:(j+1)*nbrhd_size-1) ;###Get Statistics if not a Flat region (e.g. a DROPOUT)### if ((max(subset)-min(subset)) ne 0) then begin substats=moment(subset,sdev=subsdev) sub_mean(i,j)=substats(0) sub_sdev(i,j)=subsdev sub_max(i,j)=max(subset) sub_min(i,j)=min(subset) endif endfor;j endfor;i ;###ANALYZE SUBSET STATISTICS### subsdev_moment=0 mean_subset_sdev=0 if ((max(subset)-min(subset)) ne 0) then begin subsdev_moment=moment(sub_sdev) mean_subset_sdev=median(sub_sdev) endif ;###IDENTIFY INVALID SUBNEIGHBORHOODS (OCCULTER AND DROPOUT)### for i=0,n_subsets_x-1 do begin for j=0,n_subsets_y-1 do begin ;###Identify Pixels in Subneighborhood### img_min_x_pix(i,j)=in*blocksize+$ i*nbrhd_size img_min_y_pix(i,j)=jn*blocksize+$ j*nbrhd_size img_max_x_pix(i,j)=img_min_x_pix(i,j)+$ nbrhd_size-1 img_max_y_pix(i,j)=img_min_y_pix(i,j)+$ nbrhd_size-1 ;###Handle Occulter and Max FOV### if ((minr gt 0) OR $ (maxr gt 0)) then begin r_values=r(img_min_x_pix(i,j):img_max_x_pix(i,j),$ img_min_y_pix(i,j):img_max_y_pix(i,j)) if ((min(r_values) le minr) OR $ (max(r_values) ge maxr)) then begin subset_valid(n_subsets_x*in+i,n_subsets_y*jn+j)=0 endif $ else begin subset_valid(n_subsets_x*in+i,n_subsets_y*jn+j)=1 endelse endif; minr ;###Handle Dropout Regions### zeros_count_on=0 zeros_count_off=0 junk=where((on_band(img_min_x_pix(i,j):$ img_max_x_pix(i,j),$ img_min_y_pix(i,j):$ img_max_y_pix(i,j)) eq 0),zeros_count_on) if (n_elements(off_band) gt 1) then begin junk=where((off_band(img_min_x_pix(i,j):$ img_max_x_pix(i,j),$ img_min_y_pix(i,j):$ img_max_y_pix(i,j)) eq 0),zeros_count_off) endif if ((zeros_count_off ne 0) OR $ (zeros_count_on ne 0)) then begin subset_valid(n_subsets_x*in+i,n_subsets_y*jn+j)=0 endif endfor;j endfor;i ;###IDENTIFY BAD SUBNEIGHBORHOODS### for i=0,n_subsets_x-1 do begin for j=0,n_subsets_y-1 do begin if ((sub_max(i,j) gt sub_mean(i,j)+$ stat_factor*mean_subset_sdev) $ ;Max too high OR (sub_min(i,j) lt sub_mean(i,j)-$ stat_factor*mean_subset_sdev)) $ ;Min too low then begin subset_good(n_subsets_x*in+i,n_subsets_y*jn+j)=0 endif $ else begin subset_good(n_subsets_x*in+i,n_subsets_y*jn+j)=1 endelse endfor;j endfor;i endfor;jn endfor;in print,'Setting group parameters for fixing...' print,'Counted ',n_elements(where(subset_valid eq 0)),' invalid groups (droupout/occulter).' print,'Counted ',n_elements(where(subset_good eq 0)),' bad groups (probable cosmic rays)' ;###SET GOOD/BAD/VALID FLAGS FOR EACH PIXEL### for i=0,n_subsets_total_x-1 do begin for j=0,n_subsets_total_y-1 do begin if ((subset_good(i,j) eq 0) OR $ (subset_valid(i,j) eq 0)) then begin subset_min_x_pix=i*nbrhd_size subset_min_y_pix=j*nbrhd_size subset_max_x_pix=subset_min_x_pix+(nbrhd_size-1) subset_max_y_pix=subset_min_y_pix+(nbrhd_size-1) if (subset_good(i,j) eq 0) then $ good_bad_invalid_array(subset_min_x_pix:$ subset_max_x_pix,$ subset_min_y_pix:$ subset_max_y_pix)=1 if (subset_valid(i,j) eq 0) then $ good_bad_invalid_array(subset_min_x_pix:$ subset_max_x_pix,$ subset_min_y_pix:$ subset_max_y_pix)=2 endif;subset endfor;j endfor;i tv,good_bad_invalid_array*127 ;###SET NEW IMAGE VALUES = 0 FOR BAD PIXELS new_image(where(good_bad_invalid_array eq 1))=0 ;###DO ROW CURVEFITTING### print,'doing row curvefitting' row_indices=indgen(sx) ;Indices into the row row_image=lonarr(sx,sy) row_image(*,*)=0 for i=0,sy-1 do begin row=image(*,i) row_flags=good_bad_invalid_array(*,i) start_group_flag=0 for k=0,sx-1 do begin if (start_group_flag eq 0) then begin ;look for ;new bad group if (row_flags(k) eq 1) then begin ;Start of group start_group_flag = 1 min_group_pixel=k endif;row_flags endif $ ;start_group_flag else $ if (start_group_flag eq 1) then begin ;look for ;end of group if ((row_flags(k) eq 0) OR $ (k eq (sx -1))) then begin ;End of group+1 start_group_flag = 0 if (row_flags(k) eq 0) then max_group_pixel=k-1 $ else max_group_pixel=k ;**SORT PIXELS IN TERMS OF ** ;**DISTANCE FROM THE GROUP CTR** group_ctr_pt=(max_group_pixel-$ min_group_pixel)/2+min_group_pixel sorted_x=sort(abs(row_indices-group_ctr_pt)) ;**FIND THE INDICES TO FIT** indices_to_fit=row_indices(min_group_pixel:$ max_group_pixel) indices_to_fit=indices_to_fit(where(row_flags($ indices_to_fit) eq 1)) ;Set desired range for fit on either side of group group_invalid_count=n_elements(where(row_flags($ min_group_pixel:$ max_group_pixel) eq 2,count)) desired_fitting_range=((max_group_pixel-$ min_group_pixel+1)-count) < 16 ;Find the nearest bad/invalid on low side if (min_group_pixel gt 0) then begin min_indices=where(row_flags(0:$ (min_group_pixel-1)) ne 0,count) if (count gt 0) then begin min_side_delta=min((min_group_pixel-1)-$ row_indices(min_indices)) near_inval_pixel_min=min_group_pixel-1-$ min_side_delta endif $ else near_inval_pixel_min=0 endif $ else near_inval_pixel_min=-1 ;Find the nearest bad/invalid on high side if (max_group_pixel lt (sx-1)) then begin temp_indices=where($ (row_indices gt max_group_pixel) AND $ (row_flags ne 0),count) if (count gt 0) then begin max_indices=row_indices(temp_indices) max_side_delta=min(row_indices(max_indices)-$ (max_group_pixel+1)) near_inval_pixel_max=max_group_pixel+1+$ max_side_delta endif $ else near_inval_pixel_max=(sx-1) endif $ else near_inval_pixel_max=sx ;Determine the possible ranges for the fit if (near_inval_pixel_min eq -1) then $ possible_min_range=0 $ else $ possible_min_range=(min_group_pixel-1)-(near_inval_pixel_min+1)+1 if (near_inval_pixel_max eq sx) then $ possible_max_range=0 $ else $ possible_max_range=(near_inval_pixel_max-1)-(max_group_pixel+1)+1 ;print,'possible_min_range:',possible_min_range,$ ;Determine the pixel values to use for fitting CASE 1 of (possible_min_range ge desired_fitting_range) AND $ (possible_max_range ge desired_fitting_range) : begin ;print,'Case 0' usable_min=row_indices(((min_group_pixel-1)-$ (desired_fitting_range)+1):(min_group_pixel-1)) usable_max=row_indices((max_group_pixel+1):$ ((max_group_pixel+1)+(desired_fitting_range)-1)) usable_x=[usable_min,usable_max] fit=poly_fit(usable_x,row(usable_x),$ fit_order,$ yfit,yband,fit_sigma,correl_matrix) factor=double(1.0) for p=0,fit_order do begin row_image(indices_to_fit,i)=$ long(double(row_image(indices_to_fit,$ i))+$ fit(p)*factor) factor=factor*double(indices_to_fit) endfor ;###ADD NOISE TO FIT### result=randomn(seed,n_elements(indices_to_fit)) noise=result*fit_sigma row_image(indices_to_fit,i)=$ row_image(indices_to_fit,i)+noise end (possible_min_range lt desired_fitting_range) AND $ (possible_max_range ge desired_fitting_range) : begin ;print, 'Case 1' if ((min_group_pixel gt 0) AND $ (possible_min_range gt 0)) then begin usable_min=row_indices(((min_group_pixel-1)-$ (possible_min_range)+1):(min_group_pixel-1)) lower_max_limit=max_group_pixel+1 upper_max_limit=((max_group_pixel+1)+$ (desired_fitting_range)-1)+$ (desired_fitting_range-possible_min_range) < $ ((max_group_pixel+1)+(possible_max_range)) usable_max=row_indices(lower_max_limit:upper_max_limit) usable_x=[usable_min,usable_max] endif $ else begin lower_max_limit=max_group_pixel+1 upper_max_limit=((max_group_pixel+1)+$ (2*desired_fitting_range)-1)< $ ((max_group_pixel+1)+(possible_max_range)) usable_max=row_indices(lower_max_limit:upper_max_limit) usable_x=usable_max endelse if (possible_min_range gt 0) then begin fit=poly_fit(usable_x,row(usable_x),$ fit_order,$ yfit,yband,fit_sigma,correl_matrix) factor=double(1.0) for p=0,fit_order do begin row_image(indices_to_fit,i)=$ long(double(row_image(indices_to_fit,$ i))+$ fit(p)*factor) factor=factor*double(indices_to_fit) endfor ;###ADD NOISE TO FIT### result=randomn(seed,n_elements(indices_to_fit)) noise=result*fit_sigma row_image(indices_to_fit,i)=$ row_image(indices_to_fit,i)+noise endif $ else begin fit=poly_fit(usable_x,row(usable_x),$ 1,$ yfit,yband,fit_sigma,correl_matrix) row_image(indices_to_fit,i)=$ fix(fit(0)+fit(1)*$ float(indices_to_fit)) ;###ADD NOISE TO FIT### result=randomn(seed,n_elements(indices_to_fit)) noise=result*fit_sigma row_image(indices_to_fit,i)=$ row_image(indices_to_fit,i)+noise endelse end (possible_min_range ge desired_fitting_range) AND $ (possible_max_range lt desired_fitting_range) : begin ;print,'Case 2' if ((max_group_pixel lt (sx-1)) AND $ (possible_max_range gt 0)) then begin usable_max=row_indices((max_group_pixel+1):$ ((max_group_pixel+1)+$ (possible_max_range)-1)) upper_min_limit=min_group_pixel-1 lower_min_limit=((min_group_pixel-1)-$ (desired_fitting_range)+1)-$ (desired_fitting_range-possible_max_range) > $ ((min_group_pixel-1)-(possible_min_range)) usable_min=row_indices(lower_min_limit:upper_min_limit) usable_x=[usable_min,usable_max] endif $ else begin upper_min_limit=min_group_pixel-1 lower_min_limit=((min_group_pixel-1)-$ (2*desired_fitting_range)+1) > $ ((min_group_pixel-1)-(possible_min_range)) usable_min=row_indices(lower_min_limit:upper_min_limit) usable_x=usable_min endelse if (possible_max_range gt 0) then begin fit=poly_fit(usable_x,row(usable_x),$ fit_order,$ yfit,yband,fit_sigma,correl_matrix) factor=double(1.0) for p=0,fit_order do begin row_image(indices_to_fit,i)=$ long(double(row_image(indices_to_fit,$ i))+$ fit(p)*factor) factor=factor*double(indices_to_fit) endfor ;###ADD NOISE TO FIT### result=randomn(seed,n_elements(indices_to_fit)) noise=result*fit_sigma row_image(indices_to_fit,i)=$ row_image(indices_to_fit,i)+noise endif $ else begin fit=poly_fit(usable_x,row(usable_x),$ 1,$ yfit,yband,fit_sigma,correl_matrix) row_image(indices_to_fit,i)=$ fix(fit(0)+fit(1)*$ float(indices_to_fit)) ;###ADD NOISE TO FIT### result=randomn(seed,n_elements(indices_to_fit)) noise=result*fit_sigma row_image(indices_to_fit,i)=$ row_image(indices_to_fit,i)+noise endelse end (possible_min_range lt desired_fitting_range) AND $ (possible_max_range lt desired_fitting_range) : begin ;print,'Case 3' if ((min_group_pixel gt 0) AND $ (possible_min_range gt 0)) then begin usable_min=row_indices(((min_group_pixel-1)-$ (possible_min_range)+1):$ (min_group_pixel-1)) endif $ else usable_min=[-1] if ((max_group_pixel lt (sy-1)) AND $ (possible_max_range gt 0)) then begin usable_max=row_indices((max_group_pixel+1):$ (max_group_pixel+1)+(possible_max_range)) endif $ else usable_max=[-1] CASE 1 of (usable_min(0) ge 0) AND $ (usable_max(0) ge 0) : usable_x=[usable_min,usable_max] (usable_min(0) lt 0) AND $ (usable_max(0) ge 0) : usable_x=usable_max (usable_min(0) ge 0) AND $ (usable_max(0) lt 0) : usable_x=usable_min (usable_min(0) lt 0) AND $ (usable_max(0) lt 0) : print,'no usable points for fit' else: print,'invalid case setting usable_x in case 3' endcase if (n_elements(usable_x) lt desired_fitting_range) then begin fit=poly_fit(usable_x,row(usable_x),$ 1,$ yfit,yband,fit_sigma,correl_matrix) row_image(indices_to_fit,i)=fit(0)+fit(1)*indices_to_fit ;###ADD NOISE TO FIT### result=randomn(seed,n_elements(indices_to_fit)) noise=result*fit_sigma row_image(indices_to_fit,i)=$ row_image(indices_to_fit,i)+noise endif $ else begin fit=poly_fit(usable_x,row(usable_x),$ fit_order,$ yfit,yband,fit_sigma,correl_matrix) factor=double(1.0) for p=0,fit_order do begin row_image(indices_to_fit,i)=$ long(double(row_image(indices_to_fit,i))+$ fit(p)*factor) factor=factor*double(indices_to_fit) endfor ;###ADD NOISE TO FIT### result=randomn(seed,n_elements(indices_to_fit)) noise=result*fit_sigma row_image(indices_to_fit,i)=$ row_image(indices_to_fit,i)+noise endelse end else: print,'Undefined Case' endcase endif; row_flags endif; start_group_flag endfor ;k endfor;i ;###DO COLUMN CURVEFITTING### print,'doing col curvefitting' col_indices=indgen(sy) col_image=lonarr(sx,sy) col_image(*,*)=0 for i=0,sx-1 do begin col=image(i,*) col_flags=good_bad_invalid_array(i,*) start_group_flag=0 for k=0,sy-1 do begin if (start_group_flag eq 0) then begin ;look for ;new bad group if (col_flags(k) eq 1) then begin ;Start of group start_group_flag = 1 min_group_pixel=k endif;col_flags endif $ ;start_group_flag else $ if (start_group_flag eq 1) then begin ;look for ;end of group if ((col_flags(k) eq 0) OR $ (k eq (sy - 1))) then begin ;End of group+1 or end of column start_group_flag = 0 if (col_flags(k) eq 0) then max_group_pixel=k-1 $ else max_group_pixel=k ;**SORT PIXELS IN TERMS OF ** ;**DISTANCE FROM THE GROUP CTR** group_ctr_pt=(max_group_pixel-$ min_group_pixel)/2+min_group_pixel sorted_x=sort(abs(col_indices-group_ctr_pt)) ;**FIND THE INDICES TO FIT** indices_to_fit=col_indices(min_group_pixel:$ max_group_pixel) indices_to_fit=indices_to_fit(where(col_flags($ indices_to_fit) eq 1)) ;Set desired range for fit on either side of group group_invalid_count=n_elements(where(col_flags($ min_group_pixel:$ max_group_pixel) eq 2,count)) desired_fitting_range=(max_group_pixel-$ min_group_pixel+1)-count < 16 ;Find the nearest bad/invalid on low side if (min_group_pixel gt 0) then begin min_indices=where(col_flags(0:$ (min_group_pixel-1)) ne 0,count) if (count gt 0) then begin min_side_delta=min((min_group_pixel-1)-$ col_indices(min_indices)) near_inval_pixel_min=min_group_pixel-1-$ min_side_delta endif $ else near_inval_pixel_min=0 endif $ else near_inval_pixel_min=-1 ;Find the nearest bad/invalid on high side if (max_group_pixel lt (sy-1)) then begin temp_indices=where($ (col_indices gt max_group_pixel) AND $ (col_flags ne 0),count) if (count gt 0) then begin max_indices=col_indices(temp_indices) max_side_delta=min(col_indices(max_indices)-$ (max_group_pixel+1)) near_inval_pixel_max=max_group_pixel+1+$ max_side_delta endif $ else near_inval_pixel_max=(sy-1) endif $ else near_inval_pixel_max=sy ;Determine the possible ranges for the fit if (near_inval_pixel_min eq -1) then $ possible_min_range=0 $ else $ possible_min_range=(min_group_pixel-1)-(near_inval_pixel_min+1)+1 if (near_inval_pixel_max eq sy) then $ possible_max_range=0 $ else $ possible_max_range=(near_inval_pixel_max-1)-(max_group_pixel+1)+1 ;Determine the pixel values to use for fitting CASE 1 of (possible_min_range ge desired_fitting_range) AND $ (possible_max_range ge desired_fitting_range) : begin ;print,'Col ',strcompress(string(i)),' Case 0' usable_min=col_indices(((min_group_pixel-1)-$ (desired_fitting_range)+1):(min_group_pixel-1)) usable_max=col_indices((max_group_pixel+1):$ ((max_group_pixel+1)+(desired_fitting_range)-1)) usable_x=[usable_min,usable_max] fit=poly_fit(usable_x,col(usable_x),$ fit_order,$ yfit,yband,fit_sigma,correl_matrix) factor=double(1.0) for p=0,fit_order do begin col_image(i,indices_to_fit)=$ long(double(col_image(i,indices_to_fit))+$ fit(p)*factor) factor=factor*double(indices_to_fit) endfor ;###ADD NOISE TO FIT### result=randomn(seed,n_elements(indices_to_fit)) noise=result*fit_sigma col_image(i,indices_to_fit)=$ col_image(i,indices_to_fit)+noise end (possible_min_range lt desired_fitting_range) AND $ (possible_max_range ge desired_fitting_range) : begin ;print, 'Column ',strcompress(string(i)),' Case 1' if ((min_group_pixel gt 0) AND $ (possible_min_range gt 0)) then begin usable_min=col_indices(((min_group_pixel-1)-$ (possible_min_range)+1):(min_group_pixel-1)) lower_max_limit=max_group_pixel+1 upper_max_limit=((max_group_pixel+1)+$ (desired_fitting_range)-1)+$ (desired_fitting_range-possible_min_range) < $ ((max_group_pixel+1)+(possible_max_range)) usable_max=col_indices(lower_max_limit:upper_max_limit) usable_x=[usable_min,usable_max] endif $ else begin lower_max_limit=max_group_pixel+1 upper_max_limit=((max_group_pixel+1)+$ (2*desired_fitting_range)-1)< $ ((max_group_pixel+1)+(possible_max_range)) usable_max=col_indices(lower_max_limit:upper_max_limit) usable_x=usable_max endelse if (possible_min_range gt 0) then begin fit=poly_fit(usable_x,col(usable_x),$ fit_order,$ yfit,yband,fit_sigma,correl_matrix) factor=double(1.0) for p=0,fit_order do begin col_image(i,indices_to_fit)=$ long(double(col_image(i,indices_to_fit))+$ fit(p)*factor) factor=factor*double(indices_to_fit) endfor ;###ADD NOISE TO FIT### result=randomn(seed,n_elements(indices_to_fit)) noise=result*fit_sigma col_image(i,indices_to_fit)=$ col_image(i,indices_to_fit)+noise endif $ else begin fit=poly_fit(usable_x,col(usable_x),$ 1,$ yfit,yband,fit_sigma,correl_matrix) col_image(i,indices_to_fit)=$ fix(fit(0)+fit(1)*$ float(indices_to_fit)) ;###ADD NOISE TO FIT### result=randomn(seed,n_elements(indices_to_fit)) noise=result*fit_sigma col_image(i,indices_to_fit)=$ col_image(i,indices_to_fit)+noise endelse end (possible_min_range ge desired_fitting_range) AND $ (possible_max_range lt desired_fitting_range) : begin if ((max_group_pixel lt (sy-1)) AND $ (possible_max_range gt 0)) then begin usable_max=col_indices((max_group_pixel+1):$ ((max_group_pixel+1)+$ (possible_max_range)-1)) upper_min_limit=min_group_pixel-1 lower_min_limit=((min_group_pixel-1)-$ (desired_fitting_range)+1)-$ (desired_fitting_range-possible_max_range) > $ ((min_group_pixel-1)-(possible_min_range)) usable_min=col_indices(lower_min_limit:upper_min_limit) usable_x=[usable_min,usable_max] endif $ else begin upper_min_limit=min_group_pixel-1 lower_min_limit=((min_group_pixel-1)-$ (2*desired_fitting_range)+1) > $ ((min_group_pixel-1)-(possible_min_range)) usable_min=col_indices(lower_min_limit:upper_min_limit) usable_x=usable_min endelse if (possible_max_range gt 0) then begin fit=poly_fit(usable_x,col(usable_x),$ fit_order,$ yfit,yband,fit_sigma,correl_matrix) factor=double(1.0) for p=0,fit_order do begin col_image(i,indices_to_fit)=$ long(double(col_image(i,indices_to_fit))+$ fit(p)*factor) factor=factor*double(indices_to_fit) endfor ;###ADD NOISE TO FIT### result=randomn(seed,n_elements(indices_to_fit)) noise=result*fit_sigma col_image(i,indices_to_fit)=$ col_image(i,indices_to_fit)+noise endif $ else begin fit=poly_fit(usable_x,col(usable_x),$ 1,$ yfit,yband,fit_sigma,correl_matrix) col_image(i,indices_to_fit)=$ fix(fit(0)+fit(1)*$ float(indices_to_fit)) ;###ADD NOISE TO FIT### result=randomn(seed,n_elements(indices_to_fit)) noise=result*fit_sigma col_image(i,indices_to_fit)=$ col_image(i,indices_to_fit)+noise endelse end (possible_min_range lt desired_fitting_range) AND $ (possible_max_range lt desired_fitting_range) : begin ;print,'Column '+strcompress(string(i),/remove_all)+' Case 3' if ((min_group_pixel gt 0) AND $ (possible_min_range gt 0)) then begin usable_min=col_indices(((min_group_pixel-1)-$ (possible_min_range)+1):$ (min_group_pixel-1)) endif $ else usable_min=[-1] if ((max_group_pixel lt (sy-1)) AND $ (possible_max_range gt 0)) then begin usable_max=col_indices((max_group_pixel+1):$ (max_group_pixel+1)+(possible_max_range)) endif $ else usable_max=[-1] CASE 1 of (usable_min(0) ge 0) AND $ (usable_max(0) ge 0) : usable_x=[usable_min,usable_max] (usable_min(0) lt 0) AND $ (usable_max(0) ge 0) : usable_x=usable_max (usable_min(0) ge 0) AND $ (usable_max(0) lt 0) : usable_x=usable_min (usable_min(0) lt 0) AND $ (usable_max(0) lt 0) : print,'no usable points for fit' else: print,'invalid case setting usable_x in case 3' endcase if (n_elements(usable_x) lt desired_fitting_range) then begin fit=poly_fit(usable_x,col(usable_x),$ 1,$ yfit,yband,fit_sigma,correl_matrix) col_image(i,indices_to_fit)=fit(0)+fit(1)*indices_to_fit ;###ADD NOISE TO FIT### result=randomn(seed,n_elements(indices_to_fit)) noise=result*fit_sigma col_image(i,indices_to_fit)=$ col_image(i,indices_to_fit)+noise endif $ else begin fit=poly_fit(usable_x,col(usable_x),$ fit_order,$ yfit,yband,fit_sigma,correl_matrix) factor=double(1.0) for p=0,fit_order do begin col_image(i,indices_to_fit)=$ long(double(col_image(i,indices_to_fit))+$ fit(p)*factor) factor=factor*double(indices_to_fit) endfor ;###ADD NOISE TO FIT### result=randomn(seed,n_elements(indices_to_fit)) noise=result*fit_sigma col_image(i,indices_to_fit)=$ col_image(i,indices_to_fit)+noise endelse end else: print,'Undefined Case' endcase endif; col_flags endif; start_group_flag endfor ;k endfor;i new_image=new_image+(row_image+col_image)/2 ;###HANDLE THE ROI### if keyword_set(USE_ROI) then begin new_image(roi)=patch endif ; keyword_set(use_roi) return,new_image END; FUNCTION