;+ ; NAME: ; P2SW_IMAGE_FILTER ; ; PURPOSE: ; This routine generates radial filters with slow azimuthal variation to account for large-scale ; changes in EUV irradiance from different parts of the Sun. These filters can be used to reduce ; the strong brightness gradient in large field-of-view EUV solar images to enhance visibility ; of structure and improve overall display quality. It was developed for SWAP images, but has been used ; successfully on GOES/SUVI (including large-FOV mosaics; see https://doi.org/10.1038/s41550-021-01427-8), ; STEREO/EUVI, and SolO/EUI/FSI, though these latter instruments may require some degree of manual adjustment. ; This filter will work fine on a single image, but sometimes if the filter is intended for a long ; movie, it works better to generate a median stack of all the images to be used and then use that ; to make the filter. It has been tested pretty extensively on SWAP and SUVI data, and to a limited ; extent on other datasets, but it might not work well on FITS files with nonstandard metadata. ; ; CATEGORY: ; PROBA2/SWAP Data Analysis Tools. ; ; CALLING SEQUENCE: ; filter = P2SW_IMAGE_FILTER(index, image, xcen = xcen, ycen = ycen, solrad = solrad, roll = roll, $ ; blur_factor = blur_factor, range = range, prefilter_blur = prefilter_blur, $ ; disk_smoothing = disk_smoothing, radial_smoothing = radial_smoothing, $ ; tapering_constant = tapering_constant, crush_factor = crush_factor, $ ; /fill, /w304, /gauss_blur, /regularize_edge, /taper_filter, /observer_earth) ; ; INPUTS: ; index: FITS header index structure from mreadfits (or similar) ; ; image: the image to use to generate the filter ; ; OPTIONAL INPUTS: ; xcen/ycen: center of the solar disk in pixels in x and y respectively (default is CRPIX1 and CRPIX2) ; ; solrad: solar radius in pixels (default is to calculate Rsun/CDELT1) ; ; roll: orientation of the solar north pole in degrees (default is CROTA or CROTA1/2 if available) ; (A future version of this code should also be able to obtain this from CDi_j or PCi_j matrices) ; ; blur_factor: an extra amount of blurring to reduce the local contrast-enhancing effects of the filter if ; desired. Unsmoothed filters generated from single images can also include weird artifacts ; so it is worth setting a value here (default is a 16-pixel boxcar smoothing.) ; If gauss_blur is set, this value sets the sigma of the gaussing smoothing kernel. ; ; range: the angle in degrees over which the filter is smoothed azimuthally (default is 25) ; ; prefilter_blur: width of a boxcar smoothing applied to the image before generating the filter. This is ; often unnecessary, but can be useful for some images (default is 0; no smoothing) ; ; disk_smoothing: amount to smooth the disk filter image after computing the median of features on the disk. ; This really helps reduce the contrast of features in the filtered image and prevent unwanted, ; artifacts. It is often useful EVEN if the filter itself is blurred. (default is 5-pixel boxcar) ; ; radial_smoothing: amount to smooth the radial direction when generating the filter. This suppresses radial banding, ; which often occurs when strong, localized features like AR loops appear above the limb. ; (default value is 8-pixel boxcar) ; ; tapering_parameter: the value to which to taper the edges of the filter to prevent over-enchancement of noise ; in the faint outer parts of the image. This only does anything if /taper_filter is set, which ; is strongly recommended. Default behavior is to try to guess at the right value. So far this ; seems to work ok, but this tapering requires some finesse and intuition, so the defaults ; might not work very well. If all else fails, just use trial and error and pick values near the ; mean of the darkest parts of the image. ; ; crush_factor: amount to reduce the intensity of the filter by raising it to a power. Doing this helps reduce the amount ; of filtering of the inner corona, which makes things look more "natural" by preserving some degree of ; brightness falloff in the images. You can use any number you want, but values between 0.5-1.0 work best. ; bigger than 1.0 and you will get crazy results. (default is 0.75) ; ; KEYWORD IMPUTS: ; fill: patch areas with missing or bad data with the values of their neighbors when generating the filter. ; This is not needed for every image, but sometimes really makes a difference if there are artifacts ; near the corners. If your filter looks bad in the outer reaches of the image, try setting this keyword. ; ; w304: there is a little bit of logic here that finds the limb location, which is needed to merge the disk part of ; the filter and the off-limb part. This logic doesn't work that well in the case of 304 images, where the limb ; brightens much less than in other (optically thin) EUV channels. This uses a different approach that seems ; to work better for 304. It's still experimental. If you set this keyword on a non-304 image, you will probably ; not enjoy the results very much. ; ; gauss_blur: use a gaussian smoothing instead of boxcar when smoothing the final filter. This often yields better looking ; images, but it's super slow, especially if you want a lot of blurring. Good to leave off if you are just ; trying to dial in your filter settings. ; ; regularize_edge: sometimes the interpolation used to fill image corners yields some ugly artifacts are caused by ; a sudden jump in filter strength. This is some logic that tries to fix that some. Don't set until ; you are pretty sure you need this. ; ; taper_filter: add an offset to the filter to reduce the amount of brightening in faint features. Setting this is a ; good idea (unless you plan to monkey with the filter after you generate it or something). ; ; observer_earth: if set, assume the observatory is near earth and compute solar parameters using this info instead ; of header information. ; ; OUTPUTS: ; filter: array with the filter you are going to use to normalize your image. Has the same dimensions as the input image. ; ; HOW TO USE YOUR FILTER: ; There's no hard and fast rule, but the filters are designed to reduce dynamic range by dividing the image(s) ; of interest for display. An approach that woks pretty well is to make a filtered image as in line 1 below and then ; divide the image and rescale it using the code in line 2. ; ; IDL> filter = p2sw_image_filter(index, image, /gauss, blur = 4, /taper_filter) ; IDL> filtered_image = bytscl( (image/filter)^0.25, 0.5, 2.0) ; ; You will probably have to experiment to get all the settings right, but raising the filtered image to 1/4 power ; really helps. The exact range to scale over to display the image can vary widely, but the values above have worked ; reasonably well on several different images from SWAP and SUVI. ; ; LIMITATIONS: ; Lots of things in here are still somehwat experimental. Some limited testing indicates this is not a great method ; for visible light images (both coronagraph and eclipses), but there's nothing that says you can't try it. ; We make no promises you will love the results of any particular filter, but usually with some experimentation ; you can get a decent result for all but the noisiest images. ; ; MODIFICATON HISTORY: ; Written by: Daniel B. Seaton, Royal Observatory of Belgium (Original implementation 2013) ; September 2022: First full version with several performance improvements ; Developed some documentation finally too! (dbs) ; January 2023: First formal release (dbs) ;- function p2sw_image_filter, index, image, xcen = xcen, ycen = ycen, solrad = solrad, blur_factor = blur_factor, roll = roll, $ fill = fill, range = range, w304 = w304, gauss_blur = gauss_blur, $ regularize_edge = regularize_edge, prefilter_blur = prefilter_blur, $ disk_smoothing = disk_smoothing, radial_smoothing = radial_smoothing, $ taper_filter = taper_filter, tapering_constant = tapering_constant, $ crush_factor = crush_factor, observer_earth = observer_earth ;; get metadata tag names tags = tag_names(index) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; configure keywords & defaults from optional input ;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; image properties if ~n_elements(xcen) then xcen = index.crpix1 - 1 if ~n_elements(ycen) then ycen = index.crpix2 - 1 ;; need a date to do subsequent calculations if (where(tags.matches('date_d$obs', /fold_case)))[0] ne -1 then $ obs_date = index.date_d$obs else obs_date = index.date_obs ;; images handle solar radii in too many ways, so we try to estimate it from first principles ;; if the observatory distance is available. This is crude but good enough for most observatories. if ~n_elements(solrad) then begin $ if (where(tags.matches('dsun_obs', /fold)))[0] ne -1 then begin rsun_1_au = 959.63D au_m = 149597870691D dsun_obs = index.dsun_obs solrad = (rsun_1_au / (dsun_obs/au_m))/index.cdelt1 endif else $ ;; if not, we just guess that we are at Earth, a strategy that is very ;; likely to fair poorly for stereo/solo data, but dsun_obs is usually present solrad = get_rb0p(obs_date, /radius, /quiet)/index.cdelt1 endif if ~n_elements(roll) then begin crota_position = (where(tags.matches('CROTA', /fold_case)))[0] if crota_position ne -1 then begin crota_name = tags[crota_position] status = execute('roll = -index.' + crota_name) endif else roll = 0. endif im_size = size(image) ;; processing configurables if ~n_elements(range) then range = 25. if ~n_elements(prefilter_blur) then prefilter_blur = 0 if ~n_elements(disk_smoothing) then disk_smoothing = 5 if ~n_elements(radial_smoothing) then radial_smoothing = 8 if ~n_elements(blur_factor) then blur_factor = 16 if ~n_elements(crush_factor) then crush_factor = 0.75 gauss_blur = keyword_set(gauss_blur) fill = keyword_set(fill) w304 = keyword_set(w304) regularize_edge = keyword_set(regularize_edge) taper_filter = keyword_set(taper_filter) observer_earth = keyword_set(observer_earth) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; configure required variables ;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; coordinate and distance arrays mask = fltarr(im_size[1], im_size[2]) + 1 xarr = findgen(im_size[1], im_size[2]) mod im_size[1] yarr = rotate(findgen(im_size[2], im_size[1]) mod im_size[2], 1) dist_arr = sqrt((xarr - xcen[0])^2. + (yarr - ycen[0])^2) ;; solar radius variables rsun_km = 6.955e5 rsun_px = solrad rsun = rsun_px * index.cdelt1 ;; observer location variables if (where(tags.matches('hglt_obs', /fold)))[0] ne -1 then B0 = index.hglt_obs $ else if observer_earth then result_sun = get_sun(obs_date, he_lat = B0) else B0 = 0. B0_radians = B0 * !pi/180. if (where(tags.matches('hgln_obs', /fold)))[0] ne -1 then obs_lon = index.hgln_obs $ else if observer_earth then result_sun = get_sun(obs_date, he_lon = obs_lon) else obs_lon = 0. obs_lon_radians = obs_lon * !pi/180. ;; make mask ;; This is used to separate out the disk image from the off-disk image because ;; they are processed differently. Must be sensitive to roll angle here lest ;; we end up with an image with the solar rotation misalligned to the image orientation mask[where(dist_arr gt solrad[0])] = !values.f_nan masked_image = image * mask if roll ne 0. then masked_image = rot(masked_image, roll, 1., xcen, ycen, cubic = -0.5) ;;;;;;;;;;;;;;;;;;;;;;;;; ;; Process Disk Filter ;; ;;;;;;;;;;;;;;;;;;;;;;;;; ;; here we are going to transform the disk image into a rectangular lat/lon map so we can ;; median across latitudes to make artificial coronal hole and AR tracks ;; define the map size map_size = [round(rsun_px), round(rsun_px)] ;; set up the image coordinate systems ;; methods largely from W. T. Thompson (2006), A&A 449: 791 ;; DOI: 10.1051/0004-6361:20054262 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Heliocentric Cartesian <-> Stonyhurst Heliographic ;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Heliocentric Cartesian image coords in km (x, y, z) x0 = (xarr - xcen) * index.cdelt1/rsun * rsun_km y0 = (yarr - ycen) * index.cdelt2/rsun * rsun_km ;; Thompson, p.799: if z is missing, assume r = r_sun z0 = sqrt(rsun_km^2 - x0^2 - y0^2) ;; make a map of longitudes over the disk ;; Thompson eq 12 lon_img = -(obs_lon_radians + atan(z0 * cos(B0_radians) - y0 * sin(B0_radians), x0) - !pi/2) lat_img = asin( (y0 * cos(B0_radians) + z0 * sin(B0_radians))/rsun_km ) ;; use the lat/lon array to figure out which pixels in the deprojected image we want to pull from lon_from_map = (lon_img + !pi/2) * map_size[0]/!pi lat_from_map = (lat_img + !pi/2) * map_size[1]/!pi ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Set up lat/lon rectangular solar map ;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; map coords run from center longitude +/- 90 degrees (pi/2 in radians) ;; and -90 to +90 (-pi/2 to +pi/2) in longitude lon_arr = (findgen(map_size[0], map_size[1]) mod map_size[0])/map_size[0] * !pi - !pi/2 + obs_lon_radians lat_arr = rotate((findgen(map_size[1], map_size[0]) mod map_size[1]), 1)/map_size[1] * !pi - !pi/2 ;; Thompson eq 11 x_map = rsun * cos(lat_arr) * sin(lon_arr - obs_lon_radians)/index.cdelt1 + xcen y_map = rsun * ( sin(lat_arr) * cos(B0_radians) - $ cos(lat_arr) * cos(lon_arr - obs_lon_radians) * sin(B0_radians))/index.cdelt2+ ycen ;; coordinate transforms are a pain if the image is not oriented with solar north up ;; so we take the easy way out and just rotate to north up rather than trying to adapt ;; coordinates if roll ne 0 then image_rolled = rot(image, roll, 1, xcen, ycen, cubic = -0.1, missing = 0, /pivot) $ else image_rolled = image ;; interpolate the solar image into the flat map img_proj = interpolate(image_rolled, x_map, y_map, cubic = -0.5) ;; make a median image across longitude ;; includes option to smooth the median slice a bit to suppress artifacts disk_med = median(img_proj[20:-21, *], dim = 1) smooth_disk_med = smooth(disk_med, disk_smoothing, /edge_truncate) med_proj = img_proj * 0 for col = 0, map_size[0] - 1 do med_proj[col, *] = smooth_disk_med ;; map the median smoothed image back to image, including re-rolling it med_deproj = interpolate(med_proj, lon_from_map, lat_from_map) med_deproj = rot(med_deproj, -roll, 1, xcen, ycen, /pivot, cubic = -0.5) ;; we begin setting up the filter here by replacing the actual disk with our stripes prefilter = image prefilter(where(dist_arr le solrad[0] - 0)) = med_deproj(where(dist_arr le solrad[0] - 0)) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Process Off-Disk Filter ;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; transform the image into theta/r space im_max = max( sqrt( (findgen(im_size[1]) - xcen)^2 + (findgen(im_size[2]) - ycen)^2)) unwrap = p2sw_polar_map(prefilter, 360, 300, xcen, ycen, 0, 360, 0, im_max, cubic = -0.5, missing = !values.f_nan) if regularize_edge then begin unwrap_nofill = p2sw_polar_map(prefilter, 360, 300, xcen, ycen, 0, 360, 0, im_max, cubic = -0.5, missing = !values.f_nan) reg_mask = ~finite(unwrap_nofill) mask_height = fltarr(360) for n = 0, 359 do mask_height[n] = min(where(reg_mask[n, *])) above_mask_arr = fltarr(360) for n = 0, 359 do if mask_height[n] ne -1 then above_mask_arr[n] = mean(unwrap[n, mask_height[n]: (mask_height[n] + 5)<299]) below_mask_arr = fltarr(360) for n = 0, 359 do if mask_height[n] ne -1 then below_mask_arr[n] = mean(unwrap[n, (mask_height[n]-6) > 0: (mask_height[n] -1)<299]) fix_factor = mean( (below_mask_arr/above_mask_arr)[where((mask_height ne -1) AND (above_mask_arr ge 0) AND (below_mask_arr ge 0))], /nan ) unwrap_fix = unwrap * float(1 - reg_mask) + unwrap * float(reg_mask) * fix_factor unwrap = gauss_smooth(unwrap_fix, 1, /edge_truncate, /nan) endif ;; before we make the filter, we can blur things a bit to help smooth out artifacts if prefilter_blur[0] ne 0 then unwrap = gauss_smooth(unwrap, prefilter_blur, /nan, /edge_truncate) ;; negative values in the filter are a real problem, so we set them to NANs to exclude them from the calculation unwrap(where(unwrap le 0)) = !values.f_nan ;; to avoid edge effects, we pad the transformed array on both sides so we can run the median filter (which is wide) ;; all the way across the image and keep going around at the edges where position angle reaches 0 again. we will toss ;; out the padding once we compute the radial filter, and join the edges when we map back to image coordinates. padded_unwrap = [unwrap, unwrap, unwrap] ;; if fill is set, we will set regions of the image where bad values occurred to the average of their ;; sane valued neighbors. i've had mixed results using this approach, but it sometimes help ;; fix issues with the filter due to excessive regions where a filter cannot be computed if fill then begin fill_unwrap = padded_unwrap ;; generate a list of values to be considered list = where(~finite(padded_unwrap), count) if count gt 0 then begin list_ind = array_indices(padded_unwrap, list) ;; we need only to apply our fill in the main array, and will repad with the filled version at the end sublist = list[where(list_ind[0, *] ge 360 and list_ind[0, *] le 719, subcount)] sub_ind = array_indices(padded_unwrap, sublist) ;; for every value we are going to patch.. for n = 0, subcount - 1 do begin ;; check for neighbors where there are ok values px = sub_ind(*, n) good_left = where(finite(padded_unwrap[0:px[0], px[1]]), count_left) good_right = where(finite(padded_unwrap[px[0]:*, px[1]]), count_right) + px[0] ;; if you find some pixels to work with, you can make a patch if count_left gt 0 AND count_right gt 0 then $ fill_unwrap[sublist[n]] = (padded_unwrap[max(good_left), px[1]] + padded_unwrap[min(good_right), px[1]])/2. endfor endif ;; repad with the filled values padded_unwrap = [fill_unwrap[360:719, *], fill_unwrap[360:719, *], fill_unwrap[360:719, *]] endif ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Generate the filter in radial space ;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; prepare the filter, to start, we copy the padded, unwrapped image filter_unwrap_full = padded_unwrap filter_unwrap = padded_unwrap ;; run across the unwrapped image with a median filter to create a filter that preserves radial falloff ;; but smooths the azimuthal information over a range of about 50 degrees (this could be tweaked or allowed to vary by keyword) for n = 0, 1079 do filter_unwrap_full(n, *) = median(padded_unwrap((n-range)>0:(n+range)<1079, *), dim = 1) ;; if the image we are working with is a 304 image, we have to find the limb using a different technique ;; than if it is a normal EUV image if w304 then begin profile = abs(deriv(median(filter_unwrap_full, dim = 1))) filt_max = where(max(profile) eq profile) endif else $ ;; find the limb by searching for the place where the filter (on average) is at a max ;; this is the limit of where we will fill in the inner filter image filt_max = where(max(median(filter_unwrap_full, dim = 1), /nan) eq median(filter_unwrap_full, dim = 1)) if n_elements(filt_max) ne 1 then filt_max = min(filt_max) ;; here's where we start to construct the filter for n = 0, 1079 do filter_unwrap(n, filt_max:*) = smooth(median(padded_unwrap((n-range)>0:(n+range)<1079, filt_max:*), dim = 1), $ radial_smoothing, /edge_truncate) filter_unwrap = filter_unwrap(360:719, *) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Finalize filter in normal coordinates ;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; filter = p2sw_polar_map(filter_unwrap, im_size[1], im_size[2], xcen, ycen, 0, 360, 0, im_max, cubic = 1, /inverse) ;; smoothing the filter can help remove artifacts or tune the scale of small features you enhance ;; but too much blurring leads to problems if blur_factor ne 0 then $ if gauss_blur then filter = gauss_smooth(filter, blur_factor, /edge_mirror, /nan) $ else filter = smooth(filter, blur_factor, /edge_mirror, /nan) ;;;;;;;;;;;;;;;;;;;;;;;; ;; Taper if requested ;; ;;;;;;;;;;;;;;;;;;;;;;;; if taper_filter then begin if ~n_elements(tapering_constant) then begin ;; this method seems to do ok, mostly tapering_constant = median(filter[where(filter gt 0)]) ;; experimented with histogram-based approaches, but median works ok ; hist = histogram(filter[where(filter gt 0)], binsize = 0.01, locations = xpos) ; cdf = total(hist, /cumulative)/n_elements(image) ; tapering_constant = xpos[max(where(cdf le 0.5))] endif filter = filter + tapering_constant endif ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Crush filter for more natural results ;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; filter = filter^crush_factor return, filter end