gaus_func.py

import os
from scipy import signal
import numpy as np
import matplotlib.pyplot as plt
import math
from PIL import Image
from scipy.special import comb

def gaussian(x, sigma):
    """ Return the normalized Gaussian with standard deviation sigma. """
    c = np.sqrt(2 * np.pi)
    return np.exp(-0.5 * (x / sigma)**2) / sigma / c

def gaussian2(x, sigma1, mu):
  #sigma1=1
  #mu in [2,3,5,10,20]
  #value = (mu / (2 * sigma1 * (1 / mu))) * np.exp(-np.abs(x) ** mu)
  value = (mu / (2 * sigma1 * (1 / mu))) * np.exp(-np.abs(x) ** mu)

  return value

def gaussian3(x, sigma, mu):
    """ Return the normalized Gaussian with standard deviation sigma. """
    c = np.sqrt(2 * np.pi)
    return np.exp(-0.5 * ((x - mu)/ sigma)**2) / sigma / c

def super_gaussian(x, amplitude=1.0, center=0.0, sigma=1.0, expon=2.0):
    """super-Gaussian distribution
    super_gaussian(x, amplitude, center, sigma, expon) =
        (amplitude/(sqrt(2*pi)*sigma)) * exp(-abs(x-center)**expon / (2*sigma**expon))
    """
    sigma = max(1.e-15, sigma)
    return ((amplitude/(np.sqrt(2*np.pi)*sigma))
            * np.exp(-abs(x-center)**expon / 2*sigma**expon))

#def turkey_window(x,alpha):
#    window = signal.windows.tukey(x, alpha)
#    return window

def draw_one_curve(y,im):
    print('Making image with curve')
######
    hist = im.convert('L').histogram()
    hist = (np.array(hist) / np.max(np.array(hist)))
    #hist=np.array(hist)
    hist = hist * 255
########################

    x = np.linspace(0, 255, 256)
    x = x.astype(int)

    fig = plt.figure()
    fig.patch.set_facecolor((.5, .5, .5))
    ax = fig.add_subplot()
#    f, ax = plt.subplots(1)
    ax.set_facecolor((.5, .5, .5))

    #plot histogram
    plt.plot(hist, color=(.2, .2, .2), alpha=0.3)
    plt.fill(hist, color=(.2, .2, .2), alpha=0.3)


    plt.plot(x, y, marker='o', color=(.9,.9,.9), markerfacecolor=(.1,.1,.1), markersize=3, linewidth=1, markeredgewidth=.5)
    plt.grid()
    ax.set_aspect('equal', adjustable='box')

    yt = ax.get_yticks()
    yt = []
    yt = np.append(yt, [0, 63, 127, 195, 255])
    ax.set_yticks(yt)
    # ax.set_yticklabels(np.round(yt,1))
    ax.set_yticklabels(yt.astype(int))

    ax.set_xticks(yt)
    ax.set_xticklabels(yt.astype(int))

    plt.xlim(0, 255)
    plt.ylim(0, 255)
    #plt.show()
    #plt.savefig('curve.png', bbox_inches='tight', dpi=(250))
    img = fig2img(fig)
    plt.close()
    return img

def fig2img(fig):
    """Convert a Matplotlib figure to a PIL Image and return it"""
    import io
    buf = io.BytesIO()
    fig.savefig(buf, bbox_inches='tight', dpi=(250))
    buf.seek(0)
    img = Image.open(buf)
    return img

def smoothclamp(x, mi, mx): return mi + (mx-mi)*(lambda t: np.where(t < 0 , 0, np.where( t <= 1 , 3*t**2-2*t**3, 1 ) ) )( (x-mi)/(mx-mi) )

def smoothstep(x, x_min=0, x_max=1, N=1):
    #x = np.clip((x - x_min) / (x_max - x_min), 0, 1)

    result = 0
    for n in range(0, N + 1):
         result += comb(N + n, n) * comb(2 * N + 1, N - n) * (-x) ** n

    result *= x ** (N + 1)

    return result

def make_step(width):
    x2 = np.linspace(0, 255, 256)
    x2 = x2.astype(int)

    if width>len(x2):
        wigth=len(x2)
    width=int(width)
    xstart=len(x2)/2
    xstart=int(xstart)
    #y2=x2
    y2=[1]*len(x2)
    d1=int(xstart-width/2)
    d2=int(xstart+width/2)
    #y2[:]=1
    #tw=turkey_window(d1*2,1)
    #tw=signal.windows.gaussian(d1*2+1, 30)
    #tw=signal.windows.tukey(d1*2+1, 1)
    tw=signal.windows.blackmanharris(d1*2+1)
    for step in range(0,d1,1):
       y2[step]=tw[step]

    tw_step = d1
    for step in range(d2,len(x2),1):
        y2[step]=tw[tw_step]
        tw_step=tw_step+1

    #more smoothing for curve
    y2 = signal.savgol_filter(y2, 20, 1),  # order of fitted polynomial
    y2 = np.squeeze(y2)

    y2 = y2 - np.min(y2)
    y2 = y2 / np.max(y2)
    y2 = y2 * 255
    y2 = y2.astype(int)

    print(f'xstart: {xstart}')
    print(f'y2: {y2}')
    #smoothed_step=smoothstep(y2, 0, 1, 1)
    #smoothed_step = smoothclamp(y2, 50, 128)
    #smoothed_step =math.erf(y2)
    #smoothclamp
    print(f'smoothed_step: {y2}')
    return y2

print('Working dir:')
print(os.getcwd())


im=Image.open("test.jpeg") # PIL image



#for beta in [2,3,10,20]:
#for window_size in [10,50,100,150,200]:
#      g1=make_step(window_size)
#      img=draw_one_curve(g1, im)

lut_x = [7, 57, 77, 97, 147]
lut_y = [0, 190, 255, 190, 0]

g1=make_step(100)
#img=draw_one_curve(np.array(lut_x), im)
#img.show()
#img = Image.fromarray(img,'L')

#a = np.array(img.convert('L'))
#r, g, b = img.split()
#
import splines
import cv2

#from helper import plot_spline_2d, plot_tangent_2d
points1 = [
    (-1, -0.5),
    (0, 2.3),
    (1, 1),
    (4, 1.3),
    (3.8, -0.2),
    (2.5, 0.1),
]
#s1 = splines.CatmullRom(points1, endconditions='closed')
#fig, ax = plt.subplots()
#a=s1.grid()
#img.save("curve_single.png")
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
from scipy.signal import savgol_filter
import scipy.interpolate as spi

# Example input data (replace with your actual data)
#x_data = np.array([0, 50, 100, 150, 200])
#y_data = np.array([10, 25, 40, 35, 20])
#"[3] STRICT":
#x_data = [7, 57, 77, 97, 147]
#y_data = [0, 190, 255, 190, 0]


#x_data = [0, 31, 51, 71, 121]
#y_data = [72, 190, 255, 190, 0]
x_data = [0, 10, 23, 38, 56, 75, 95, 120,255]
y_data = [255, 198, 141, 77, 35, 12, 2, 0,0]

# Interpolate to 0-255 values of x
x_interp = np.linspace(0, 255, 256)
#‘linear’, ‘nearest’, ‘nearest - up’, ‘zero’, ‘slinear’, ‘quadratic’, ‘cubic’, ‘previous’, or ‘next’.‘zero’, ‘slinear’, ‘quadratic’ and ‘cubic’
f = interp1d(x_data, y_data, kind='quadratic', fill_value="extrapolate")  # Use cubic interpolation

y_interp = f(x_interp)
#intfunc = spi.interp1d(x_interp,y_data,fill_value="extrapolate")
#y_interp = intfunc(x_new)
y_interp[y_interp < 0] = 0
y_interp[y_interp > 255] = 255
# Smooth the interpolated data
y_smooth = savgol_filter(y_interp, window_length=50, polyorder=2)
y_smooth[y_smooth<0] =0
y_smooth[y_smooth>255] = 255

# Plot the results
plt.figure(figsize=(8, 6))

plt.plot(x_interp, y_interp, '-', label='Interpolated Data', marker='o')
print(f'max_y: {max(y_interp)}')
plt.plot(x_interp, y_smooth, '-', label='Smoothed Data', marker='x')
plt.plot(x_data, y_data, '-', label='Original Data', marker='^', markersize=10)
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Interpolation and Smoothing')
plt.legend()
plt.grid(True)
#plt.show()

def catmull_rom_one_point(x, v0, v1, v2, v3):
    """Computes interpolated y-coord for given x-coord using Catmull-Rom.

    Computes an interpolated y-coordinate for the given x-coordinate between
    the support points v1 and v2. The neighboring support points v0 and v3 are
    used by Catmull-Rom to ensure a smooth transition between the spline
    segments.
    Args:
        x: the x-coord, for which the y-coord is needed
        v0: 1st support point
        v1: 2nd support point
        v2: 3rd support point
        v3: 4th support point
    """
    c1 = 1. * v1
    c2 = -.5 * v0 + .5 * v2
    c3 = 1. * v0 + -2.5 * v1 + 2. * v2 -.5 * v3
    c4 = -.5 * v0 + 1.5 * v1 + -1.5 * v2 + .5 * v3
    return (((c4 * x + c3) * x + c2) * x + c1)

def catmull_rom(p_x, p_y, res):
    """Computes Catmull-Rom Spline for given support points and resolution.

    Args:
        p_x: array of x-coords
        p_y: array of y-coords
        res: resolution of a segment (including the start point, but not the
            endpoint of the segment)
    """
    # create arrays for spline points
    x_intpol = np.empty(res*(len(p_x)-1) + 1)
    y_intpol = np.empty(res*(len(p_x)-1) + 1)

    # set the last x- and y-coord, the others will be set in the loop
    x_intpol[-1] = p_x[-1]
    y_intpol[-1] = p_y[-1]

    # loop over segments (we have n-1 segments for n points)
    for i in range(len(p_x)-1):
        # set x-coords
        x_intpol[i*res:(i+1)*res] = np.linspace(
            p_x[i], p_x[i+1], res, endpoint=False)
        if i == 0:
            # need to estimate an additional support point before the first
            y_intpol[:res] = np.array([
                catmull_rom_one_point(
                    x,
                    p_y[0] - (p_y[1] - p_y[0]), # estimated start point,
                    p_y[0],
                    p_y[1],
                    p_y[2])
                for x in np.linspace(0.,1.,res, endpoint=False)])
        elif i == len(p_x) - 2:
            # need to estimate an additional support point after the last
            y_intpol[i*res:-1] = np.array([
                catmull_rom_one_point(
                    x,
                    p_y[i-1],
                    p_y[i],
                    p_y[i+1],
                    p_y[i+1] + (p_y[i+1] - p_y[i]) # estimated end point
                ) for x in np.linspace(0.,1.,res, endpoint=False)])
        else:
            y_intpol[i*res:(i+1)*res] = np.array([
                catmull_rom_one_point(
                    x,
                    p_y[i-1],
                    p_y[i],
                    p_y[i+1],
                    p_y[i+2]) for x in np.linspace(0.,1.,res, endpoint=False)])


    return (x_intpol, y_intpol)

res=50
lut_x = [0,58, 108, 128, 148, 198,255]
lut_y = [0,0, 190, 255, 190, 0,0]

x_intpol, y_intpol = catmull_rom(lut_x, lut_y, res)
# fancy plotting
#plt.figure()
#plt.scatter(lut_x, lut_y)
#plt.plot(x_intpol, y_intpol)
#plt.show()
#==============================================================================

#from PIL import Image, ImageCms

#brightness_lab=1.5
#image_lab = im.convert('RGB', 'Lab')
#image_lab = im.convert('RGB').convert('Lab')

# Convert to Lab colourspace
# srgb_p = ImageCms.createProfile("sRGB")
# lab_p  = ImageCms.createProfile("LAB")
#
# rgb2lab = ImageCms.buildTransformFromOpenProfiles(srgb_p, lab_p, "RGB", "LAB")
# lab2rgb = ImageCms.buildTransformFromOpenProfiles(lab_p, srgb_p,  "LAB","RGB")
# Lab = ImageCms.applyTransform(im, rgb2lab)
#
# L, a, b = Lab.split()

# L.save('L.png')
# a.save('a.png')
# b.save('b.png')
#L.show()

#image_lab=Lab
# Adjust the L (Lightness) channel
#image_lab[:, :, 0] = image_lab[:, :, 0] * brightness_lab

# Clip values to valid range (0-100 for L)
#image_lab[:, :, 0] = np.clip(image_lab[:, :, 0], 0, 100)

# Convert back to RGB
#image_rgb = image_lab.convert(image_lab, 'Lab', 'RGB')
#modified_image = Image.fromarray(np.uint8(image_rgb * 255))  # back to Pillow
#modified_image.show()

#https://pyimagesearch.com/2015/10/05/opencv-gamma-correction/

img = cv2.imread('test.jpeg')
value=100

gamma=2.2
invGamma = 1.0 / gamma

table = np.array([((i / 255.0) ** invGamma) * 255
                  for i in np.arange(0, 256)]).astype("uint8")
# apply gamma correction using the lookup table



RGB = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
#RGB=cv2.add(RGB, value)
#RGB=cv2.pow(RGB, gamma)
RGB=cv2.LUT(RGB, table)
#cv2.pow
cv2.imwrite('image_rgb_br+100.jpg',cv2.cvtColor(RGB, cv2.COLOR_RGB2BGR))

LAB = cv2.cvtColor(img, cv2.COLOR_BGR2LAB)
# Split into L, a, b channels
l, a, b = cv2.split(LAB)

# Adjust L channel (lightness)
#l = cv2.add(l, value)
#l = cv2.pow(l, gamma)
l=cv2.LUT(l, table)
#return cv2.LUT(image, table)

# Clip L channel values to valid range (0-255)
l = np.clip(l, 0, 255)

# Merge channels back into LAB image
lab_image = cv2.merge((l, a, b))
LAB=lab_image
print(np.min(LAB))
print(np.max(LAB))



BGR = cv2.cvtColor(LAB, cv2.COLOR_LAB2BGR)
# img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
cv2.imwrite('image_lab_br+100.jpg',BGR)
cv2.imshow('preview', BGR)
cv2.waitKey(0)