SamWavelet.py

import pywt
import numpy as np
import cv2 as cv
from matplotlib import pyplot as plt
import math



def adjust_image_size(image):
    """
    Entree : une image sous forme de matrice
    Sortie : l'image avec des dimensions divisibles par 2
    """
    row, col = image.shape
    # S'assurer que les dimensions sont divisibles par 2 (pour éviter les erreurs lors de la transformation)
    new_row = row if row % 2 == 0 else row - 1
    new_col = col if col % 2 == 0 else col - 1
    return image[:new_row, :new_col]

def bruit(imag,var):
    """
    Entree : img : image sous forme de matrice
             var : variance du bruit gaussien
    Sortie : image bruitée par bruit gaussien de variance var
    """
    row, col = imag.shape  # récupération des dimensions de la matrice de l'image
    mean = 0
    sigma = var**0.5
    gauss = np.random.normal(mean, sigma, (row, col))  # on génère du bruit gaussien
    image_bruit = imag + gauss  # construction de l'image bruitée par addition des matrices
    image_bruit = np.clip(image_bruit, 0, 255).astype(np.uint8)  # Assurer que les valeurs restent dans [0, 255] et soient entières
    return image_bruit


def bruit_speckle(image):
    row, col = image.shape
    gauss = np.random.randn(row, col)
    noisy = image + image * gauss
    return noisy.astype(np.uint8)

def SNR(image_or, image_br):
    """
    Entree : image_or : image originale
             image_br : image bruitée
    Sortie : rapport signal beuit en dB (SNR) calculé à partir de l'erreur quadratique moyenne
    """
    # Rogner les images pour qu'elles aient les mêmes dimensions
    min_row = min(image_or.shape[0], image_br.shape[0])
    min_col = min(image_or.shape[1], image_br.shape[1])
    image_or = image_or[:min_row, :min_col]
    image_br = image_br[:min_row, :min_col]

    erreur_quad_moy = np.mean((image_or.astype(np.float32) - image_br.astype(np.float32)) ** 2)
    val_max_pix = 255.0
    SNR_bruit = 20 * math.log10(val_max_pix / (erreur_quad_moy**0.5))
    return SNR_bruit

def zoom(image, fact):
    """
    Entree : image : image sous forme de matrice
             fact : facteur de zoom
    Sortie :  Zoom de l'image par un facteur fact (à partir du coin gauche de l'image)
    """
    img_zoom=[]
    ligne=[]
    li, col=image.shape
    for i in range(li//fact):
        for j in range(col//fact):
            ligne.append(image[i][j])
        img_zoom.append(ligne)
        ligne=[]
    return img_zoom

def intervalle_coef(L):
    """
    Entree : une liste de listes (tableau numpy) 
    Sortie : un tuple contenant le minimum et le maximum golbaux de tous les éléments contenus dans la liste
    """
    k=[]
    m=[]
    for i in range(len(L)):
        k.append(max(L[i]))
        m.append(min(L[i]))
    return (min(m), max(k))

def contraste(img):
  """"
  Entree : image en niveaux de gris
  Sortie : valeur du contraste de l'image
  """
  max= np.max(img)
  min= np.min(img)
  return (max-min)/(max+min)



def hard_thresholding(wavelet_transform, threshold):
    """
    Applique un seuillage dur sur une transformée en ondelettes.

    Paramètres:
    -----------
    wavelet_transform : numpy.ndarray
        La transformée en ondelettes à traiter.
    threshold : float
        Le seuil à appliquer.

    Retourne:
    ---------
    numpy.ndarray
        La transformée en ondelettes après seuillage dur.
    """
    return np.where(np.abs(wavelet_transform) > threshold, wavelet_transform, 0)

def soft_thresholding(coeffs, threshold):
    """
    Entree : coeffs : les quatres listes de la décomposition en ondelettes
             thresgold : seuil de filtrage
    Sortie : les coefficients filtrés par seuillage doux
    """

    # tableau pour stocker les coefficients seuils
    coeffs_thresh = np.zeros_like(coeffs)

    #définition de la fonction de seuillage doux
    coeffs_thresh[(coeffs > threshold)] = coeffs[(coeffs > threshold)] - threshold
    coeffs_thresh[(coeffs < -threshold)] = coeffs[(coeffs < -threshold)] + threshold

    return coeffs_thresh

def debruitage_ht(image, threshold):
    """
    Entree : image : image à débruiter
                threshold : seuil de filtrage 
    Sortie : image débruitée par seuillage dur
    """

    # Décomposer l'image en ondelettes
    cA, (cH, cV, cD) = pywt.dwt2(image, 'db2')

    # Appliquer le hard thresholding sur les coefficients de détail (trois dernieres listes)
    cH_thresh = hard_thresholding(cH, threshold)
    cV_thresh = hard_thresholding(cV, threshold)
    cD_thresh = hard_thresholding(cD, threshold)

    # Reconstruction de l'image avec les coefficients modifiés
    coeffs2 = (cA, (cH_thresh, cV_thresh, cD_thresh))
    img_debruite = pywt.idwt2(coeffs2, 'db2')

    return img_debruite

def debruitage_st(image, threshold):
    """
    Entree : image : image à débruiter
                threshold : seuil de filtrage
    Sortie : image débruitée par seuillage doux
    """

    # Décomposer l'image en ondelettes
    cA, (cH, cV, cD) = pywt.dwt2(image, 'db2')

    # Appliquer le hard thresholding sur les coefficients de détail
    cH_thresh = soft_thresholding(cH, threshold)
    cV_thresh = soft_thresholding(cV, threshold)
    cD_thresh = soft_thresholding(cD, threshold)

    # Reconstruction de l'image avec les coefficients modifiés
    coeffs2 = (cA, (cH_thresh, cV_thresh, cD_thresh))
    img_debruite = pywt.idwt2(coeffs2, 'db2')

    return img_debruite

def triFus_SNR(L, img):
    n=len(L)
    if n<=1:
        return L
    return fusion_SNR(triFus_SNR(L[:n//2],img), triFus_SNR(L[n//2:],img), img)

def fusion_SNR(l1,l2, img):
    if not l1:
        return l2
    if not l2:
        return l1
    if SNR(img, l1[0])>SNR(img, l2[0]):
        return [l1[0]]+fusion_SNR(l1[1:],l2, img)
    return [l2[0]]+fusion_SNR(l2[1:],l1, img)

def analyse(img, var, alpha):
  """"
  Entree : image en niveaux de gris, variance du bruit, alpha est le facteur de contraste qui détermine le contraste minimal souhaité
  Sortie : seuil optimal
  La fonction utilise deux critères : Le pSNR comme critère principal et veille
  à ce que l'image finale ne perde pas trop de contraste
  """
  # Ajouter du bruit à l'image
  img_br = bruit(img, var)
  cA, (cH, cV, cD) = pywt.dwt2(img_br, 'db2')
  intervalle = intervalle_coef(cH)
  img_contraste=contraste(img)

  # Générer les seuils
  ########Choix de la distribution des seuils à étudier 
  thr = np.logspace(0.1*np.sqrt(var), intervalle[1]+0.5*np.sqrt(var), num=20)

  # Appliquer le débruitage par ondelettes avec hard thresholding
  dwt = [debruitage_ht(img_br, t) for t in thr]

  max_snr=0
  thresh=None
  for i in range(len(dwt)):
        current_snr = SNR(dwt[i], img)
        current_contrast = contraste(dwt[i])

        if (abs(current_contrast - img_contraste) <= alpha * img_contraste) and current_snr > max_snr: #vérification du contraste et du snr
            max_snr = current_snr
            thresh = thr[i]

  
  if thresh is None : #dans le cas où il n'y a pas de seuil acceptable compte tenu des critères de contraste et de SNR, on choisit le seuil qui maximise le SNR
    print("pas de seuil convenable")
    maxi=0
    indmax=0
    for i in range(len(dwt)):
        if SNR(dwt[i], img)>maxi:
            maxi=SNR(dwt[i], img)
            indmax=i
    return thr[indmax]
  print("seuil convenable trouvé : ", thresh)
  return thresh

def seuil_visushrink(image,var, wavelet='db2'):
    """ 
    Entree : image : image à débruiter
             var : variance du bruit
             wavelet : ondelette à utiliser
    Sortie : seuil optimal pour la méthode de VisuShrink
    """
    imgbr= bruit(image, var)
    cA, (cH, cV, cD) = pywt.dwt2(imgbr, wavelet)

    N=cD.size
    ecart_typ= np.sqrt(var)
    # Estimer l'écart-type du bruit à partir des coefficients de la première échelle
    threshold = np.sqrt(2 * np.log(N)*ecart_typ)

    return threshold

def affiche(signaux,titre, x):
    """
    Entree :signaux : liste de signaux à afficher
            titre : titre du graphique
            x : 1 pour un seul signal, 0 pour plusieurs signaux
    Sortie : affichage du signal

    Ce programme permet d'économiser du code en évitant de répéter les mêmes instructions pour afficher un ou plusieurs signaux
    """
    if not titre:
        titre= "Signal"
    if x==1:
        plt.plot(signaux)
        plt.title(titre)
        plt.show()
    else:
        for s in signaux: 
            plt.plot(s)
            plt.title(titre)
            plt.show()

def bruit1D(signal, var):
    """
    Entree : signal : signal 1D
             var : variance du bruit gaussien
    Sortie : signal bruité par bruit gaussien de variance var
    """
    sigma = var**0.5
    gauss= np.random.normal(0, sigma, len(signal))
    signal_bruit = signal + gauss  # construction de l'image bruitée par addition des matrices
    
    return signal_bruit

def bruit_f(s):
    """
    Entree : s : un signal
    Sortie : un signal bruité fortement
    """
    sigma = 0.5
    bruit = sigma * np.random.randn(len(s))
    return s + bruit

def ht1D(signal, seuil):
    """
    Entree : signal : signal 1D
             thresh : seuil de filtrage
    Sortie : signal débruité par seuillage dur
    """
    return signal * (np.abs(signal) > seuil)

def ht1D_ca(signal, thresh):
    """
    Entree : signal : signal 1D
             thresh : seuil de filtrage
    Sortie : signal débruité par seuillage dur
    """
    cA,cD= pywt.dwt(signal, 'db2')
    cA_thresh = hard_thresholding(cA, thresh)
    signal_debruite = pywt.idwt(cA_thresh,cD, 'db2')
    return signal_debruite

def st1D(signal, thresh):
    """
    Entree : signal : signal 1D
             thresh : seuil de filtrage
    Sortie : signal débruité par seuillage doux
    """
    cA,cD= pywt.dwt(signal, 'db2')
    cD_thresh = soft_thresholding(cD, thresh)
    signal_debruite = pywt.idwt(cA, cD_thresh, 'db2')
    return signal_debruite

def bandeau(coeffD):
    """
    Entree : coeffD : coefficients de détails
    Sortie : bandeau bas et haut dans un tuple (bandeau bas, bandeau haut)
    """
    return (np.mean(coeffD)-2*np.std(coeffD), np.mean(coeffD)+2*np.std(coeffD))

def erreur_bandeau(signal, bandeau_haut, bandeau_bas):
    """
    Entree : signal : signal 1D
             bandeau_haut : valeur du haut du bandeau
             bandeau_bas : valeur du bas du bandeau
    Sortie : erreur du bandeau
    """
    erreur=0
    for i in signal:
        if i>bandeau_haut or i<bandeau_bas:
            erreur+=1
    return erreur/len(signal)

def erreur_qad_bandeau(signal, bandeau_haut, bandeau_bas):
    """
    Entree : signal : signal 1D
             bandeau_haut : valeur du haut du bandeau
             bandeau_bas : valeur du bas du bandeau
    Sortie : erreur quadratique moyenne du bandeau
    """
    erreur=0
    for i in signal:
        if i>bandeau_haut or i<bandeau_bas:
            erreur+= (i-np.mean(signal))**2
    return erreur/len(signal)

def SNR1D(signal, signal_bruit):
    """
    Entree : signal : signal 1D
             signal_bruit : signal 1D bruité
    Sortie : rapport signal bruit en dB (SNR) calculé à partir de l'erreur quadratique moyenne
    """
    erreur_quad_moy = np.mean((signal.astype(np.float32) - signal_bruit.astype(np.float32)) ** 2)  # calcul de l'erreur quadratique moyenne
    val_max_pix = max(signal)
    SNR_bruit = 20 * math.log10(val_max_pix / (erreur_quad_moy**0.5))  # formule du SNR en dB

    return SNR_bruit 

def analyse1D(sign_br, var):
    """
    Entree : signal : signal 1D
             
    Sortie : seuil optimal pour le débruitage
    """  
    cA, cD= pywt.dwt(sign_br, 'db2') #décomposition en ondelettes
    intervalle= (np.mean(cD)-2*np.std(cD), np.mean(cD)+2*np.std(cD)) #on crée les bornes de l'intervalle avec le critère précédent
    thr = np.logspace(0.1*np.sqrt(var), intervalle[1]+0.5*np.sqrt(var), num=20) #on se donne l'intervalle des seuils à tester  

    dwt = [ht1D(sign_br, t) for t in thr] #on applique le seuillage dur pour chaque coefficient de détail

    max_snr=0
    thresh=None
    for i in range(len(dwt)):
        current_snr = SNR1D(dwt[i], sign_br)
        if current_snr > max_snr:
            max_snr = current_snr
            thresh = thr[i]
    return thresh

def regression_poly_cd(cd, deg):
    """
    Entree : cd_sorted : liste triée des coefficients de détail
    Sortie : tableau des valeurs de la régression linéaire
    """
    cd=sorted(cd)
    x = np.arange(len(cd))
    p = np.polyfit(x, cd, deg)  
    y = np.polyval(p, x)
    return y

def analyse1D_regression(sign_br):
    """
    Entree : signal : signal 1D
             
    Sortie : seuil optimal pour le débruitage
    """  
    cA, cD= pywt.dwt(sign_br, 'db2') #décomposition en ondelettes
    thr = regression_lineaire_cd(np.sort(cD), 10)
    thr = [val for val in thr if val > 0]
    dwt = [ht1D(sign_br, t) for t in thr] #on applique le seuillage dur pour chaque seuil

    max_snr=0
    thresh=None
    for i in range(len(dwt)):
        current_snr = SNR1D(dwt[i], sign_br)
        if current_snr > max_snr:
            max_snr = current_snr
            thresh = thr[i]
    return thresh




    


def perform_wavelet_transform(x, Jmin, direction, options=None):
    """
    Execute une transformation en ondelettes 1D ou 2D avec lissage symetrique
    et traitement des bords symetriques.
    
    Paramètres:
    -----------
    x : numpy.ndarray
        Signal ou image d'entrée
    Jmin : int
        Niveau de décomposition minimum
    direction : int
        1 pour transformée directe, -1 pour transformée inverse
    options : dict, optional
        liste des options:
        - filter: string ou array, coefficients du filtre d'ondelettes
          Peut être 'haar', 'linear'/'5-3', ou '9-7'
        - ti: bool, si vrai, calculer la transformée invariante de translation
        - separable: bool, si vrai, calculer la transformée 2D separable
        - ndims: int, force la dimension (1 ou 2) si elle est mal détectée
    
    Sortie:
    --------
    numpy.ndarray
        Coefficients d'ondelettes ou signal reconstruit
    
    """
        
    if options is None:
        options = {}
    
    # Options par défaut
    h = options.get('filter', '9-7')
    separable = options.get('separable', False)
    
    # Détection de la dimension
    d = options.get('ndims', -1)
    if d < 0:
        if x.ndim == 1 or (x.ndim == 2 and (x.shape[0] == 1 or x.shape[1] == 1)):
            d = 1
        else:
            d = 2
    
    # assurer que x est un tableau numpy
    x = np.array(x, dtype=float)
    
    # Gérer l'orientation du vecteur 1D
    if d == 1 and x.ndim == 2 and x.shape[0] == 1:
        x = x.T
    
    # Convertir les noms de filtres en chaînes en coefficients
    if isinstance(h, str):
        h = h.lower()
        if h == 'haar':
            if d == 1 or not separable:
                return perform_haar_transform(x, Jmin, direction, options)
        elif h in ['linear', '5-3']:
            h = np.array([1/2, 1/4, np.sqrt(2)])
        elif h in ['9-7', '7-9']:
            h = np.array([1.586134342, -0.05298011854, -0.8829110762, 0.4435068522, 1.149604398])
        else:
            raise ValueError(f"Unknown filter: {h}")
    
    # Option translation invariant
    ti = options.get('ti', False)
    
    # Gérer la transformée 2D separable
    if d == 2 and separable:
        options_copy = options.copy()
        options_copy['ti'] = False
        if ti:
            print("Warning: Separable does not work for translation invariant transform")
        
        n = x.shape[0]
        if direction == 1:
            for i in range(n):
                x[:, i] = perform_wavelet_transform(x[:, i], Jmin, direction, options_copy)
            for i in range(n):
                x[i, :] = perform_wavelet_transform(x[i, :], Jmin, direction, options_copy)
        else:
            for i in range(n):
                x[i, :] = perform_wavelet_transform(x[i, :], Jmin, direction, options_copy)
            for i in range(n):
                x[:, i] = perform_wavelet_transform(x[:, i], Jmin, direction, options_copy)
        return x
    
    # Calculer les étapes de lifting
    n = x.shape[0]
    m = (len(h) - 1) // 2
    Jmax = int(math.log2(n)) - 1
    jlist = list(range(Jmax, Jmin-1, -1))
    if direction == -1:
        jlist = jlist[::-1]
    
    if not ti:
        # Transformée subsamplée
        for j in jlist:
            if d == 1:
                x[:2**(j+1)] = lifting_step(x[:2**(j+1)], h, direction)
            else:
                x[:2**(j+1), :2**(j+1)] = lifting_step(x[:2**(j+1), :2**(j+1)], h, direction)
                x[:2**(j+1), :2**(j+1)] = lifting_step(x[:2**(j+1), :2**(j+1)].T, h, direction).T
    else: #i.e if ti is True
        #----préparer les données pour la transformée invariante de translation----
        nJ = Jmax - Jmin + 1
        if direction == 1:
            if d == 1: #pour un signal 1D,
                x = np.repeat(x[:, np.newaxis], nJ+1, axis=1) #on crée une matrice 3D avec nJ+1 colonnes
            elif d == 2: #pour une image 2D,
                x_3d = np.zeros((x.shape[0], x.shape[1], 3*nJ+1)) #on crée une matrice 3D avec 3*nJ+1 colonnes
                x_3d[:,:,0] = x #on copie la matrice d'entrée dans la première colonne de la matrice 3D
                x = x_3d #on affecte la matrice 3D à x
        
        for j in jlist:
            dist = 2**(Jmax-j) #distance entre les coefficients
            if d == 1: #pour un signal 1D,
                if direction == 1: #pour une transformée directe,
                    result = lifting_step_ti(x[:, 0], h, direction, dist) #on applique la transformée invariante de translation
                    x[:, 0] = result[:, 0] #on affecte la première colonne de la matrice 3D à x
                    x[:, j-Jmin+1] = result[:, 1] #on affecte la deuxième colonne de la matrice 3D à x
                else: #pour une transformée inverse,
                    x[:, 0] = lifting_step_ti(x[:, [0, j-Jmin+1]], h, direction, dist) #on applique la transformée inverse
            else:
                dj = 3*(j-Jmin)
                if direction == 1:
                    # Transformée horizontale
                    result = lifting_step_ti(x[:,:,0], h, direction, dist)
                    x[:,:,0] = result[:,:,0]
                    x[:,:,dj+1] = result[:,:,1]
                    
                    # Transformée verticale sur le signal original
                    result = lifting_step_ti(x[:,:,0].T, h, direction, dist)
                    x[:,:,0] = result[:,:,0].T
                    x[:,:,dj+2] = result[:,:,1].T
                    
                    # Transformée verticale sur les détails horizontaux
                    result = lifting_step_ti(x[:,:,dj+1].T, h, direction, dist)
                    x[:,:,dj+1] = result[:,:,0].T
                    x[:,:,dj+3] = result[:,:,1].T
                else:
                    # Inverse l'ordre des opérations
                    x[:,:,dj+1] = x[:,:,dj+1].T
                    x[:,:,dj+3] = x[:,:,dj+3].T
                    x[:,:,dj+1] = lifting_step_ti(x[:,:,[dj+1, dj+3]], h, direction, dist).T
                    
                    x[:,:,0] = x[:,:,0].T
                    x[:,:,dj+2] = x[:,:,dj+2].T
                    x[:,:,0] = lifting_step_ti(x[:,:,[0, dj+2]], h, direction, dist).T
                    
                    x[:,:,0] = lifting_step_ti(x[:,:,[0, dj+1]], h, direction, dist)
        
        if direction == -1:
            if d == 1:
                x = x[:, 0]
            else:
                x = x[:,:,0]
    
    return x


def lifting_step(x, h, direction):
    """
    Effectue une étape de lifting pour la transformée en ondelettes subsamplée.
    
    Paramètres:
    -----------
    x : numpy.ndarray
        Signal ou partie de la transformée
    h : numpy.ndarray
        Coefficients du filtre de lifting
    direction : int
        1 pour transformée directe, -1 pour transformée inverse
    
    Returns:
    --------
    numpy.ndarray
        Résultat après l'étape de lifting
    """
    m = (len(h) - 1) // 2
    
    if direction == 1:
        #séparer les coefs
        d = x[1::2].copy()
        s = x[::2].copy()
        
        # lissage
        for i in range(m):
            # prédiction
            s_extended = np.concatenate([s, s[-1:]])
            d = d - h[2*i] * (s + s_extended[1:])
            
            # amelioration
            d_extended = np.concatenate([d[:1], d])
            s = s + h[2*i+1] * (d + d_extended[:-1])
        
        # scalling
        s = s * h[-1]
        d = d / h[-1]
        
        # Fusionner dans le domaine de la transformée
        return np.concatenate([s, d])
    else:
        # Séparer dans le domaine de la transformée
        half = len(x) // 2
        s = x[:half].copy() / h[-1]
        d = x[half:].copy() * h[-1]
        
        # étapes de lifting inverses
        for i in range(m-1, -1, -1):
            # amelioration inverse
            d_extended = np.concatenate([d[:1], d])
            s = s - h[2*i+1] * (d + d_extended[:-1])
            
            # prédiction inverse
            s_extended = np.concatenate([s, s[-1:]])
            d = d + h[2*i] * (s + s_extended[1:])
        
        # Fusionner
        result = np.zeros(len(x))
        result[::2] = s
        result[1::2] = d
        
        return result

def lifting_step_ti(x, h, direction, dist):
    """
    Effectue une étape de lifting pour la transformée en ondelettes invariante de translation.
    
    Paramètres:
    -----------
    x : numpy.ndarray
        Signal ou partie de la transformée
    h : numpy.ndarray
        Coefficients du filtre de lifting
    direction : int
        1 pour transformée directe, -1 pour transformée inverse
    dist : int
        Paramètre de distance pour la translation
    
    Returns:
    --------
    numpy.ndarray
        Résultat après l'étape de lifting
    """
    m = (len(h) - 1) // 2
    n = x.shape[0]
    
    # Calculer les indices décalés avec les conditions aux limites
    s1 = np.arange(n) + dist
    s2 = np.arange(n) - dist
    
    # Appliquer les conditions aux limites
    s1[s1 >= n] = 2*n - s1[s1 >= n] - 1
    s1[s1 < 0] = -s1[s1 < 0] - 1
    s2[s2 >= n] = 2*n - s2[s2 >= n] - 1
    s2[s2 < 0] = -s2[s2 < 0] - 1
    #cette étape consiste à étendre les coefficients de manière symétrique aux niveaux des bords.
    
    if direction == 1:
        # séparer les coefs
        d = x.copy()
        s = x.copy()
        
        # étapes de lifting
        for i in range(m):
            d = d - h[2*i] * (s[s1] + s[s2])
            s = s + h[2*i+1] * (d[s1] + d[s2])
        
        # scalling
        s = s * h[-1]
        d = d / h[-1]
        
        # Retourner les deux parties de fréquence basse et haute
        if x.ndim == 1:
            return np.column_stack([s, d])
        else:
            result = np.zeros((x.shape[0], x.shape[1], 2))
            result[:,:,0] = s
            result[:,:,1] = d
            return result
    else:
        # Pour la transformée inverse, x doit avoir deux composants
        if x.ndim == 2 and x.shape[1] == 2:
            s = x[:, 0].copy() / h[-1]
            d = x[:, 1].copy() * h[-1]
        else:
            s = x[:,:,0].copy() / h[-1]
            d = x[:,:,1].copy() * h[-1]
        
        # étapes de lifting inverses
        for i in range(m-1, -1, -1):
            s = s - h[2*i+1] * (d[s1] + d[s2])
            d = d + h[2*i] * (s[s1] + s[s2])
        
        # Fusionner
        return (s + d) / 2


def perform_haar_transform(x, Jmin, direction, options=None):
    """
    Perform Haar wavelet transform.
    
    Parameters:
    -----------
    x : numpy.ndarray
        Input signal or image
    Jmin : int
        Minimum decomposition level
    direction : int
        1 for forward transform, -1 for inverse transform
    options : dict, optional
        Dictionary of options
    
    Returns:
    --------
    numpy.ndarray
        Wavelet coefficients or reconstructed signal
    """
    if options is None:
        options = {}
    
    # Detect dimensionality
    d = options.get('ndims', -1)
    if d < 0:
        if x.ndim == 1 or (x.ndim == 2 and (x.shape[0] == 1 or x.shape[1] == 1)):
            d = 1
        else:
            d = 2
    
    # Ensure x is a numpy array
    x = np.array(x, dtype=float)
    
    # Handle 1D vector orientation
    if d == 1 and x.ndim == 2 and x.shape[0] == 1:
        x = x.T
    
    # Translation invariant option
    ti = options.get('ti', False)
    
    # Compute Haar transform
    n = x.shape[0]
    Jmax = int(math.log2(n)) - 1
    jlist = list(range(Jmax, Jmin-1, -1))
    if direction == -1:
        jlist = jlist[::-1]
    
    if not ti:
        # Standard Haar transform
        for j in jlist:
            if d == 1:
                size = 2**(j+1)
                if direction == 1:
                    # Forward transform
                    temp = x[:size].copy()
                    x[:size//2] = (temp[::2] + temp[1::2]) / np.sqrt(2)
                    x[size//2:size] = (temp[::2] - temp[1::2]) / np.sqrt(2)
                else:
                    # Inverse transform
                    temp = x[:size].copy()
                    x[0:size:2] = (temp[:size//2] + temp[size//2:size]) / np.sqrt(2)
                    x[1:size:2] = (temp[:size//2] - temp[size//2:size]) / np.sqrt(2)
            else:
                size = 2**(j+1)
                # Apply to rows
                for i in range(size):
                    if direction == 1:
                        temp = x[i, :size].copy()
                        x[i, :size//2] = (temp[::2] + temp[1::2]) / np.sqrt(2)
                        x[i, size//2:size] = (temp[::2] - temp[1::2]) / np.sqrt(2)
                    else:
                        temp = x[i, :size].copy()
                        x[i, 0:size:2] = (temp[:size//2] + temp[size//2:size]) / np.sqrt(2)
                        x[i, 1:size:2] = (temp[:size//2] - temp[size//2:size]) / np.sqrt(2)
                
                # Apply to columns
                for i in range(size):
                    if direction == 1:
                        temp = x[:size, i].copy()
                        x[:size//2, i] = (temp[::2] + temp[1::2]) / np.sqrt(2)
                        x[size//2:size, i] = (temp[::2] - temp[1::2]) / np.sqrt(2)
                    else:
                        temp = x[:size, i].copy()
                        x[0:size:2, i] = (temp[:size//2] + temp[size//2:size]) / np.sqrt(2)
                        x[1:size:2, i] = (temp[:size//2] - temp[size//2:size]) / np.sqrt(2)
    else:
        # Translation invariant Haar transform
        # This is a simplified implementation
        nJ = Jmax - Jmin + 1
        if direction == 1:
            if d == 1:
                result = np.zeros((n, nJ+1))
                result[:, 0] = x
                
                for j_idx, j in enumerate(jlist):
                    size = 2**(j+1)
                    shift = 2**(Jmax-j)
                    
                    # No shift
                    temp = result[:size, 0].copy()
                    result[:size//2, j_idx+1] = (temp[::2] + temp[1::2]) / np.sqrt(2)
                    result[size//2:size, j_idx+1] = (temp[::2] - temp[1::2]) / np.sqrt(2)
                    
                    # With shift
                    temp = np.roll(result[:size, 0], shift//2)
                    shifted_low = (temp[::2] + temp[1::2]) / np.sqrt(2)
                    shifted_high = (temp[::2] - temp[1::2]) / np.sqrt(2)
                    
                    # Average the results
                    result[:size//2, 0] = (result[:size//2, j_idx+1] + np.roll(shifted_low, -shift//4)[:size//2]) / 2
                
                return result
            else:
                # 2D TI transform is more complex and not fully implemented here
                pass
        else:
            # Inverse TI transform
            if d == 1:
                # Simplified inverse
                result = x[:, 0].copy()
                
                for j_idx, j in enumerate(reversed(jlist)):
                    size = 2**(j+1)
                    low = result[:size//2].copy()
                    high = x[:, j_idx+1][size//2:size].copy()
                    
                    # Reconstruct
                    result[0:size:2] = (low + high) / np.sqrt(2)
                    result[1:size:2] = (low - high) / np.sqrt(2)
                
                return result
            else:
                # 2D inverse TI transform not fully implemented
                pass
    
    return x


    def plot_wavelet(wt, Jmin):
        """
    Visualise les coefficients d'ondelettes.
    
    Parameters:
    -----------
    wt : numpy.ndarray
        Coefficients d'ondelettes
    Jmin : int
        Niveau de décomposition minimum
    """
    plt.plot(wt)
    plt.title("Coefficients d'ondelettes")
    plt.grid(True)

def create_dyadic_grid(wt, Jmin):
    """
    Crée une grille dyadique pour visualiser les coefficients d'ondelettes.
    
    Parameters:
    -----------
    wt : numpy.ndarray
        Coefficients d'ondelettes
    Jmin : int
        Niveau de décomposition minimum
    
    Returns:
    --------
    numpy.ndarray
        Grille dyadique des coefficients
    """
    n = int(np.log2(len(wt)))
    # Extraire les coefficients par niveau
    matrice2 = np.zeros((n, 2**(n-1)))
    for j in range(n):
        matrice2[j, :2**j] = wt[2**j:2**(j+1)]
    
    # Créer la grille dyadique avec redondance
    matrice3 = np.zeros((n, 2**(n-1)))
    for j in range(n):
        for k in range(2**j):
            block_size = int(2**(n-1) / 2**j)
            matrice3[j, k*block_size:(k+1)*block_size] = matrice2[j, k]
    
    return matrice3

def plot_dyadic_grid(matrice, title, filename, cmap='viridis'):
    """
    Affiche une grille dyadique avec des couleurs nettes et des séparations de cases visibles.
    
    Parameters:
    -----------
    matrice : numpy.ndarray
        Matrice contenant la grille dyadique
    title : str
        Titre du graphique
    filename : str
        Nom du fichier pour sauvegarder l'image
    cmap : str, optional
        Colormap à utiliser (par défaut 'viridis')
    """
    n, width = matrice.shape
    
    # Créer une figure avec des dimensions adaptées
    fig, ax = plt.subplots(figsize=(10, 6))
    
    # Utiliser pcolormesh pour un affichage plus net des blocs
    # Créer une grille de coordonnées pour pcolormesh
    x = np.arange(0, width + 1)
    y = np.arange(0, n + 1)
    
    # Afficher la matrice avec pcolormesh pour des blocs nets
    im = ax.pcolormesh(x, y, matrice, cmap=cmap, shading='flat')
    
    # Ajouter une barre de couleur
    cbar = plt.colorbar(im, ax=ax)
    
    # Configurer les axes
    ax.set_yticks(np.arange(0.5, n + 0.5))
    ax.set_yticklabels(range(1, n + 1))
    
    # Ajouter un titre et des étiquettes d'axes
    plt.title(title, fontsize=14)
    plt.xlabel('Position', fontsize=12)
    plt.ylabel('Échelle', fontsize=12)
    
    # Inverser l'axe y pour que l'échelle 1 soit en haut
    ax.invert_yaxis()
    
    # Ajuster les limites des axes
    ax.set_xlim(0, width)
    ax.set_ylim(n, 0)
    
    # Supprimer les graduations mineures
    ax.minorticks_off()
    
    # Enregistrer et afficher
    plt.tight_layout()
    plt.savefig(filename, dpi=150)
    plt.show()