In [217]:
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
In [2]:
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
In [3]:
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.]
In [4]:
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.]
In [5]:
reffile=glob.glob('./reference/2A5A9259.CR3')
files=glob.glob("./selected/*CR3")
files.sort()
In [6]:
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)
In [7]:
refimage.min(),refimage.max(),refimage.shape
Out[7]:
(0, 65535, (4498, 6742, 3))
In [8]:
plt.imshow(refimage/refimage.max())
Out[8]:
<matplotlib.image.AxesImage at 0x7fb35964cd60>
In [9]:
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)
In [10]:
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]
In [296]:
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
In [187]:
rad.min(),rad.max()
Out[187]:
(nan, nan)
In [188]:
refimage.min(),refimage.max()
Out[188]:
(0, 65535)
In [195]:
ggg.min(),ggg.max()
Out[195]:
(0.0, 0.9999694819562066)
In [259]:
ggg=color.rgb2gray(refimage)
edges = canny(ggg, sigma=5, low_threshold=0.02,high_threshold=0.1)
In [202]:
refimage.min(),refimage.max(),ggg.min(),ggg.max()
Out[202]:
(0, 65535, 0.0, 0.9999694819562066)
In [225]:
kernal=np.array([[1,1,1],[1,-8,1],[1,1,1]])
edges=ndimage.convolve(ggg,kernal)
In [260]:
plt.imshow(edges,vmax=0.001,vmin=-0.001,cmap='bwr')
plt.xlim([3500,5000])
plt.ylim([1500,3000])
Out[260]:
(1500.0, 3000.0)
In [231]:
mask=edges>0.005
plt.imshow(mask,vmax=1,vmin=0,cmap='bwr')
plt.xlim([3500,5000])
plt.ylim([1500,3000])
Out[231]:
(1500.0, 3000.0)
In [291]:
xc,yc,rad,xx,yy=cal_center(refimage)
In [297]:
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)
In [298]:
# xarr=np.linspace(0,1000,1000)
plt.plot(10**fr(rad[2100,:]))
# plt.plot(rad[2100,:])
Out[298]:
[<matplotlib.lines.Line2D at 0x7fb36ae0ab20>]
In [300]:
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)
In [266]:
clearimg.max(),clearimg.min(),sharp_img.max(),sharp_img.min()
Out[266]:
(1.0, 0.0, 1.0, 0.0)
In [299]:
# plt.subplot(121)
plt.plot(bgrgb[2100,:])
plt.yscale('log')
# plt.subplot(122)
# plt.imshow(rad)
In [303]:
clearimg.min(),clearimg.max(),sharp_img.min(),sharp_img.max(),bgrgb.min(),bgrgb.max()
Out[303]:
(0.0, 0.0, 0.0, 0.0, 0, 13800)
In [305]:
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).
Out[305]:
<matplotlib.image.AxesImage at 0x7fb319f2a490>
In [311]:
np.nanmin(aa),np.nanmax(aa)
Out[311]:
(0.0, 1872.4285714285713)
In [315]:
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).
Out[315]:
(1500.0, 3000.0)
In [13]:
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)
In [ ]:
In [14]:
x=np.arange(hdr.shape[0])
y=np.arange(hdr.shape[1])
yy,xx=np.meshgrid(y,x)
In [124]:
xc = np.nanmean(xx[edges])-5
yc = np.nanmean(yy[edges])
In [125]:
# 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
In [126]:
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).
Out[126]:
(1700.0, 2700.0)
In [127]:
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)
In [128]:
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]
In [129]:
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')
In [130]:
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)
In [131]:
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
In [132]:
In [133]:
plt.plot(bgrgb[2000,10:])
plt.yscale('log')
In [134]:
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)
In [135]:
# 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)
In [140]:
sharp_img = unsharp_mask(tmp, radius=8, amount=10)
In [141]:
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).
Out[141]:
(1000.0, 3500.0)
In [109]:
from skimage.restoration import (
    denoise_tv_chambolle,
    denoise_bilateral,
    denoise_wavelet,
    estimate_sigma,
)
In [46]:
aa=denoise_tv_chambolle(tmp, weight=0.1, channel_axis=-1)
In [60]:
sharp_img = unsharp_mask(tmp, radius=2, amount=10)
In [61]:
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).
Out[61]:
(500.0, 4000.0)
In [44]:
# 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)
In [68]:
color.rgb2gray(tmp).min(),color.rgb2gray(tmp).max()
Out[68]:
(0.0, 0.48040574223597254)
In [79]:
sharp_img = unsharp_mask(color.rgb2gray(tmp), radius=20, amount=10)
In [80]:
plt.figure(figsize=(10,10))
plt.imshow((sharp_img),vmax=5e-2,vmin=0)
plt.xlim([2500,5500])
plt.ylim([500,4000])
Out[80]:
(500.0, 4000.0)
In [ ]: