import sys,os
# sys.path.append("/Users/yangkai/Software/rawpy")
from PIL import Image
import numpy as np
import rawpy
import rawpy.enhance
import matplotlib.pyplot as plt
from scipy.fft import fft2, ifft2, fftshift
import glob
import cv2 as cv
from skimage import data, filters
from tqdm import tqdm
from scipy import interpolate
from scipy.interpolate import splrep, BSpline
from scipy import signal,stats
from scipy.signal import convolve
from skimage import restoration
from skimage import color
from skimage.transform import hough_circle, hough_circle_peaks
from skimage.feature import canny
from skimage.draw import circle_perimeter
from skimage.filters import unsharp_mask
from scipy import ndimage
def get_filtered(image, cutoffs, squared_butterworth=True, order=3.0, npad=0):
"""Lowpass and highpass butterworth filtering at all specified cutoffs.
Parameters
----------
image : ndarray
The image to be filtered.
cutoffs : sequence of int
Both lowpass and highpass filtering will be performed for each cutoff
frequency in `cutoffs`.
squared_butterworth : bool, optional
Whether the traditional Butterworth filter or its square is used.
order : float, optional
The order of the Butterworth filter
Returns
-------
lowpass_filtered : list of ndarray
List of images lowpass filtered at the frequencies in `cutoffs`.
highpass_filtered : list of ndarray
List of images highpass filtered at the frequencies in `cutoffs`.
"""
lowpass_filtered = []
highpass_filtered = []
for cutoff in cutoffs:
lowpass_filtered.append(
filters.butterworth(
image,
cutoff_frequency_ratio=cutoff,
order=order,
high_pass=False,
squared_butterworth=squared_butterworth,
npad=npad,
)
)
highpass_filtered.append(
filters.butterworth(
image,
cutoff_frequency_ratio=cutoff,
order=order,
high_pass=True,
squared_butterworth=squared_butterworth,
npad=npad,
)
)
return lowpass_filtered, highpass_filtered
def cross_correlation_using_fft(x, y):
f1 = fft2(x)
f2 = fft2(np.flipud(np.fliplr(y)))
cc = np.real(ifft2(f1 * f2))
return fftshift(cc) # Shift zero frequency component to the center of the spectrum
def compute_shift(x, y):
# Assuming x and y are the same size
correlation = cross_correlation_using_fft(x, y)
maxima = np.unravel_index(np.argmax(correlation), correlation.shape)
mid_points = np.array([np.floor(dim/2) for dim in x.shape])
shifts = maxima - mid_points
return shifts
# Example usage:
# x and y are two images (2D arrays) that you want to compare
# They need to be the same size and preferably normalized or preprocessed similarly
x = np.random.rand(256, 256) # Simulated image
y = np.roll(x, shift=(5, 10), axis=(0, 1)) # Simulated shifted image of x
shifts = compute_shift(x, y)
print("Estimated shifts:", shifts)
Estimated shifts: [ -6. -11.]
exptime=[1/6.0]
exptime_all=[1/20.,1/13.0,1/20.0,1/13.0,1/10.0,1/10.,0.1,1/8.,1/13.,1/13.,1/13.,1/13.,\
1/13.,1/13.,1/13.,1/50.0,1/6.,1/6.,1/6.,1/6.,1/6.,1/6.,0.25,0.25,0.25,0.2,0.1,1/6.,0.2,\
0.2,0.2,0.2,0.2,0.2,0.2,1/15.,0.2,1/3.,1/3.,1/3.,1/3.,1/3.,1/3.,1/3.,1/6.,0.2,0.2,0.2,\
0.2,0.2,0.2,0.25,0.2,0.2,1./3.,1/3.,0.1,1/13.,1/15.,1/15.,1/15.,1/15.,1/15.,1/15.,1/15.,1/15.,1/15.,1/15.,1/15.]
reffile=glob.glob('./reference/2A5A9259.CR3')
files=glob.glob("./selected/*CR3")
files.sort()
for path in reffile:
with rawpy.imread(path) as raw:
refimage= raw.postprocess(highlight_mode=rawpy.HighlightMode.Reconstruct(3),output_bps=16,use_camera_wb=True,four_color_rgb=True)
# refimage= raw.postprocess(gamma=(1,1), highlight_mode=rawpy.HighlightMode.Reconstruct(3),output_bps=16,use_camera_wb=True)
refimage.min(),refimage.max(),refimage.shape
(0, 65535, (4498, 6742, 3))
plt.imshow(refimage/refimage.max())
<matplotlib.image.AxesImage at 0x7fb35964cd60>
for path in reffile:
with rawpy.imread(path) as raw:
refimage= raw.postprocess(gamma=(1,1), highlight_mode=rawpy.HighlightMode.Reconstruct(3),output_bps=16,use_camera_wb=True)
hdr=refimage*exptime[0]
for ix in tqdm(range(len(files))):
path=files[ix]
with rawpy.imread(path) as raw:
tmp = raw.postprocess(gamma=(1,1), highlight_mode=rawpy.HighlightMode.Reconstruct(3),output_bps=16,use_camera_wb=True)
shift=compute_shift(refimage[...,0],tmp[...,0])
tmp1=np.copy(tmp)
for ix in range(3):
yy=np.roll(tmp[...,ix], shift=(int(shift[0]),int(shift[1])), axis=(0, 1))
tmp1[...,ix]=yy
hdr+=tmp1*exptime_all[ix]
100%|██████████| 69/69 [08:04<00:00, 7.03s/it]
def cal_center(imgin):
sunr=355
ggg=color.rgb2gray(imgin)
edges = canny(ggg, sigma=5, low_threshold=0.02,high_threshold=0.1)
hough_res = hough_circle(edges, sunr)
accums, cx, cy, radii = hough_circle_peaks(hough_res, [sunr-10,sunr+10], total_num_peaks=2)
x=np.arange(imgin.shape[0])
y=np.arange(imgin.shape[1])
yy,xx=np.meshgrid(y,x)
xc = np.nanmean(xx[edges])
yc = np.nanmean(yy[edges])
rad=np.sqrt((xx-xc)**2+(yy-yc)**2)
return xc,yc,rad,xx,yy
def remove_bg(imgin,radIn):
radbin=np.linspace(0,4000,4000)
bg_r=np.zeros(len(radbin)-1)
bg_g=np.zeros(len(radbin)-1)
bg_b=np.zeros(len(radbin)-1)
for ix in tqdm(range(1,4000)):
mask= (radIn>radbin[ix-1]) & (radIn<radbin[ix])
bg_r[ix-1]=np.median(1.0*imgin[...,0][mask])
bg_g[ix-1]=np.median(1.0*imgin[...,1][mask])
bg_b[ix-1]=np.median(1.0*imgin[...,2][mask])
rr=0.5*(radbin[1:]+radbin[:-1])
fr = interpolate.interp1d(rr, np.log10(bg_r), kind='cubic',fill_value=np.nan,bounds_error=False)
fg = interpolate.interp1d(rr, np.log10(bg_g), kind='cubic',fill_value=np.nan,bounds_error=False)
fb = interpolate.interp1d(rr, np.log10(bg_b), kind='cubic',fill_value=np.nan,bounds_error=False)
bgrgb=np.zeros_like(imgin)
bgrgb[:,:,0]=10**fr(radIn.flatten()).reshape(bgrgb.shape[0],bgrgb.shape[1])
bgrgb[:,:,1]=10**fg(radIn.flatten()).reshape(bgrgb.shape[0],bgrgb.shape[1])
bgrgb[:,:,2]=10**fb(radIn.flatten()).reshape(bgrgb.shape[0],bgrgb.shape[1])
# bgrgb=10.0**bgrgb
tmp=imgin/bgrgb
tmp/=np.nanmax(tmp)
tmp[~np.isfinite(tmp)]=0.0
return tmp,bgrgb,fr,fg,fb
rad.min(),rad.max()
(nan, nan)
refimage.min(),refimage.max()
(0, 65535)
ggg.min(),ggg.max()
(0.0, 0.9999694819562066)
ggg=color.rgb2gray(refimage)
edges = canny(ggg, sigma=5, low_threshold=0.02,high_threshold=0.1)
refimage.min(),refimage.max(),ggg.min(),ggg.max()
(0, 65535, 0.0, 0.9999694819562066)
kernal=np.array([[1,1,1],[1,-8,1],[1,1,1]])
edges=ndimage.convolve(ggg,kernal)
plt.imshow(edges,vmax=0.001,vmin=-0.001,cmap='bwr')
plt.xlim([3500,5000])
plt.ylim([1500,3000])
(1500.0, 3000.0)
mask=edges>0.005
plt.imshow(mask,vmax=1,vmin=0,cmap='bwr')
plt.xlim([3500,5000])
plt.ylim([1500,3000])
(1500.0, 3000.0)
xc,yc,rad,xx,yy=cal_center(refimage)
clearimg,bgrgb,fr,fg,fb=remove_bg(refimage,rad)
100%|██████████| 3999/3999 [01:38<00:00, 40.69it/s] /var/folders/rr/2n3c73qx6r141zscx2v__r9w0000gn/T/ipykernel_56177/656547796.py:30: RuntimeWarning: invalid value encountered in cast bgrgb[:,:,0]=10**fr(radIn.flatten()).reshape(bgrgb.shape[0],bgrgb.shape[1]) /var/folders/rr/2n3c73qx6r141zscx2v__r9w0000gn/T/ipykernel_56177/656547796.py:31: RuntimeWarning: invalid value encountered in cast bgrgb[:,:,1]=10**fg(radIn.flatten()).reshape(bgrgb.shape[0],bgrgb.shape[1]) /var/folders/rr/2n3c73qx6r141zscx2v__r9w0000gn/T/ipykernel_56177/656547796.py:32: RuntimeWarning: invalid value encountered in cast bgrgb[:,:,2]=10**fb(radIn.flatten()).reshape(bgrgb.shape[0],bgrgb.shape[1]) /var/folders/rr/2n3c73qx6r141zscx2v__r9w0000gn/T/ipykernel_56177/656547796.py:34: RuntimeWarning: divide by zero encountered in divide tmp=imgin/bgrgb /var/folders/rr/2n3c73qx6r141zscx2v__r9w0000gn/T/ipykernel_56177/656547796.py:34: RuntimeWarning: invalid value encountered in divide tmp=imgin/bgrgb /var/folders/rr/2n3c73qx6r141zscx2v__r9w0000gn/T/ipykernel_56177/656547796.py:35: RuntimeWarning: invalid value encountered in divide tmp/=np.nanmax(tmp)
# xarr=np.linspace(0,1000,1000)
plt.plot(10**fr(rad[2100,:]))
# plt.plot(rad[2100,:])
[<matplotlib.lines.Line2D at 0x7fb36ae0ab20>]
sharp_img = unsharp_mask(clearimg, radius=2, amount=10)
/opt/anaconda3/envs/works/lib/python3.8/site-packages/skimage/_shared/utils.py:348: RuntimeWarning: Images with dimensions (M, N, 3) are interpreted as 2D+RGB by default. Use `multichannel=False` to interpret as 3D image with last dimension of length 3. return func(*args, **kwargs)
clearimg.max(),clearimg.min(),sharp_img.max(),sharp_img.min()
(1.0, 0.0, 1.0, 0.0)
# plt.subplot(121)
plt.plot(bgrgb[2100,:])
plt.yscale('log')
# plt.subplot(122)
# plt.imshow(rad)
clearimg.min(),clearimg.max(),sharp_img.min(),sharp_img.max(),bgrgb.min(),bgrgb.max()
(0.0, 0.0, 0.0, 0.0, 0, 13800)
plt.figure(figsize=(20,10))
plt.subplot(131)
plt.imshow(refimage/np.nanmax(refimage)*1e1)
plt.subplot(132)
plt.imshow(clearimg*1e3)
plt.subplot(133)
plt.imshow(refimage/bgrgb*1e1)
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). /var/folders/rr/2n3c73qx6r141zscx2v__r9w0000gn/T/ipykernel_56177/1272407391.py:9: RuntimeWarning: divide by zero encountered in divide plt.imshow(refimage/bgrgb*1e1) /var/folders/rr/2n3c73qx6r141zscx2v__r9w0000gn/T/ipykernel_56177/1272407391.py:9: RuntimeWarning: invalid value encountered in divide plt.imshow(refimage/bgrgb*1e1) Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
<matplotlib.image.AxesImage at 0x7fb319f2a490>
np.nanmin(aa),np.nanmax(aa)
(0.0, 1872.4285714285713)
aa=np.divide(refimage,bgrgb)
aa[~np.isfinite(aa)]=0.0
plt.imshow(aa/np.nanmax(aa)*1e3)
plt.xlim([3500,5000])
plt.ylim([1500,3000])
/var/folders/rr/2n3c73qx6r141zscx2v__r9w0000gn/T/ipykernel_56177/3957215381.py:1: RuntimeWarning: divide by zero encountered in divide aa=np.divide(refimage,bgrgb) /var/folders/rr/2n3c73qx6r141zscx2v__r9w0000gn/T/ipykernel_56177/3957215381.py:1: RuntimeWarning: invalid value encountered in divide aa=np.divide(refimage,bgrgb) Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
(1500.0, 3000.0)
sunr=355
ggg=color.rgb2gray(hdr)
edges = canny(ggg, sigma=8, low_threshold=10000,high_threshold=50000)
hough_res = hough_circle(edges, sunr)
accums, cx, cy, radii = hough_circle_peaks(hough_res, [sunr-10,sunr+10], total_num_peaks=2)
x=np.arange(hdr.shape[0])
y=np.arange(hdr.shape[1])
yy,xx=np.meshgrid(y,x)
xc = np.nanmean(xx[edges])-5
yc = np.nanmean(yy[edges])
# x=np.arange(hdr.shape[0])
# y=np.arange(hdr.shape[1])
# yy,xx=np.meshgrid(y,x)
# xc,yc=2220,4210
# xc,yc=cy[0],cx[0]
# xc,yc=2223,4520
rad=np.sqrt((xx-xc)**2+(yy-yc)**2)
# mask0=np.exp(-(rad-6*sunr)**2/(0.5*(10*sunr)**2)**2)
# mask0=1/(rad/sunr+1)**4
# mask1=((np.tanh((rad/sunr - 3))+1)*0.5)**2
# mask=mask0*mask1
plt.figure(figsize=(10,4))
plt.imshow(5*hdr/hdr.max())
plt.plot([yc],[xc],'-x')
plt.contour(rad,levels=[20,100,300,355,400])
plt.xlim([3500,4800])
plt.ylim([1700,2700])
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
(1700.0, 2700.0)
radbin=np.linspace(0,4000,4000)
bg_r=np.zeros(len(radbin)-1)
bg_g=np.zeros(len(radbin)-1)
bg_b=np.zeros(len(radbin)-1)
for ix in tqdm(range(1,4000)):
mask= (rad>radbin[ix-1]) & (rad<radbin[ix])
bg_r[ix-1]=np.median(hdr[...,0][mask])
bg_g[ix-1]=np.median(hdr[...,1][mask])
bg_b[ix-1]=np.median(hdr[...,2][mask])
100%|██████████| 3999/3999 [01:44<00:00, 38.11it/s]
rr=0.5*(radbin[1:]+radbin[:-1])
plt.plot(rr/sunr,bg_r,'--r',label='R')
plt.plot(rr/sunr,bg_g,'--g',label='G')
plt.plot(rr/sunr,bg_b,'--b',label='B')
plt.legend()
plt.yscale('log')
rr=0.5*(radbin[1:]+radbin[:-1])
fr = interpolate.interp1d(rr, np.log10(bg_r), kind='cubic',fill_value=np.nan,bounds_error=False)
fg = interpolate.interp1d(rr, np.log10(bg_g), kind='cubic',fill_value=np.nan,bounds_error=False)
fb = interpolate.interp1d(rr, np.log10(bg_b), kind='cubic',fill_value=np.nan,bounds_error=False)
bgrgb=np.zeros_like(hdr)
for ix in range(bgrgb.shape[0]):
bgrgb[ix,:,0]=fr(rad[ix,:])
bgrgb[ix,:,1]=fg(rad[ix,:])
bgrgb[ix,:,2]=fb(rad[ix,:])
bgrgb=10**bgrgb
plt.plot(bgrgb[2000,10:])
plt.yscale('log')
sx, sy = 21,21
X, Y = np.ogrid[0:sx, 0:sy]
psf = stats.norm.pdf(np.sqrt((X - sx/2)**2 + (Y - sy/2)**2), 0, 10)
tmp=hdr/bgrgb
tmp/=np.nanmax(tmp)
tmp[~np.isfinite(tmp)]=0.0
# deconvolved_sun=np.copy(tmp)
# deconvolved_sun[...,0]= restoration.richardson_lucy(deconvolved_sun[...,0], psf, num_iter=30)
# deconvolved_sun[...,1]= restoration.richardson_lucy(deconvolved_sun[...,1], psf, num_iter=30)
# deconvolved_sun[...,2]= restoration.richardson_lucy(deconvolved_sun[...,2], psf, num_iter=30)
# plt.figure(figsize=(20,10))
# plt.subplot(121)
# plt.imshow(tmp*4e1)
# plt.xlim([3000,5500])
# plt.ylim([1000,3500])
# # plt.axis('off')
# # plt.savefig("eclipse_hdr.png", bbox_inches='tight',dpi=300)
# plt.subplot(122)
# plt.imshow(deconvolved_sun*4e1)
# plt.xlim([3000,5500])
# plt.ylim([1000,3500])
# # plt.axis('off')
# # plt.savefig("eclipse_hdr.png", bbox_inches='tight',dpi=300)
sharp_img = unsharp_mask(tmp, radius=8, amount=10)
plt.figure(figsize=(20,10))
plt.subplot(121)
plt.imshow(tmp*4e1)
plt.xlim([3000,5500])
plt.ylim([1000,3500])
# plt.axis('off')
# plt.savefig("eclipse_hdr.png", bbox_inches='tight',dpi=300)
plt.subplot(122)
plt.imshow(sharp_img*2e1)
plt.xlim([3000,5500])
plt.ylim([1000,3500])
# plt.axis('off')
# plt.savefig("eclipse_hdr.png", bbox_inches='tight',dpi=300)
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
(1000.0, 3500.0)
from skimage.restoration import (
denoise_tv_chambolle,
denoise_bilateral,
denoise_wavelet,
estimate_sigma,
)
aa=denoise_tv_chambolle(tmp, weight=0.1, channel_axis=-1)
sharp_img = unsharp_mask(tmp, radius=2, amount=10)
plt.figure(figsize=(10,10))
plt.imshow(sharp_img*2e1)
plt.xlim([2500,5500])
plt.ylim([500,4000])
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
(500.0, 4000.0)
# aa=denoise_tv_chambolle(result_3, weight=0.1, channel_axis=-1)
# bb=denoise_bilateral(result_3, sigma_color=0.05, sigma_spatial=15, channel_axis=-1)
# cc=denoise_wavelet(result_3, channel_axis=-1, rescale_sigma=True)
# dd=denoise_tv_chambolle(result_3, weight=0.2, channel_axis=-1)
# ee=denoise_bilateral(result_3, sigma_color=0.1, sigma_spatial=15, channel_axis=-1)
# ff=denoise_wavelet(result_3, channel_axis=-1, convert2ycbcr=True, rescale_sigma=True)
color.rgb2gray(tmp).min(),color.rgb2gray(tmp).max()
(0.0, 0.48040574223597254)
sharp_img = unsharp_mask(color.rgb2gray(tmp), radius=20, amount=10)
plt.figure(figsize=(10,10))
plt.imshow((sharp_img),vmax=5e-2,vmin=0)
plt.xlim([2500,5500])
plt.ylim([500,4000])
(500.0, 4000.0)