Soutenez-nous
Publicité
+ Répondre à la discussion
Affichage des résultats 1 à 5 sur 5
  1. #1
    Rédacteur/Modérateur

    Avatar de millie
    Profil pro
    Inscrit en
    juin 2006
    Messages
    6 940
    Détails du profil
    Informations personnelles :
    Localisation : Luxembourg

    Informations forums :
    Inscription : juin 2006
    Messages : 6 940
    Points : 8 763
    Points
    8 763

    Par défaut Filtre de choc de Osher et Rudin

    Bonjour,

    Je vais vous donner l'algorithme pour appliquer un filtre de choc de base.
    C'est à dire le filtre de choc de Osher et Rudin qui consiste à déterminer l'état stationnaire de l'équation suivante :

    dI/dt = - sign( Dnn(I)) * || grad I ||
    (ou Dnn est la dérivée directionnelle selon le gradient et grad est le gradient) (et I l'image).

    Il y a d'autres filtres de chocs (qui font par exemple intervenir un arctang ou des complexes, ou encore une convolution avec une gaussienne). Ces autres méthodes sont souvent de meilleur qualité car elles peuvent être appliquées à des images ayant de bruit.

    Le filtre de choc de base a tendance à ne pas trop aimer le bruit.

    La méthode de base qui est souvent choisi pour résoudre et de calculer :

    Code :
    1
    2
    3
    4
    5
    /* I est l'image d'entrée*/
    I_0 = I
    
    I_(n+1) = I_n - dt * sign( Dnn(I)) * || grad I ||
    Il suffit d'itérer cette équation (une 20 ou 100 aines de fois), mais attention la discrétisation de grad et de Dnn demande quelques précautions.

    Tout d'abord :

    Voici les fonctions permettant de calculer la différentielle en un point dans un canal particulier (je suis désolé d'avoir choisi une notation C++ mais ça me donne moins de boulot à traduire. J'espère en tout cas que vous comprendrez la sémantique des opérateurs :
    Code :
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    95
    96
    97
    98
    99
    100
    101
    102
    103
    104
    105
    106
    107
    108
    109
    110
    111
    112
    113
    114
    115
    116
    117
    118
    119
    120
    121
    122
    123
    124
    125
    126
    127
    128
    129
    130
    131
    132
    133
    134
    135
    136
    137
    138
    139
    140
    141
    142
    143
    144
    145
    146
    147
    148
    149
    150
    151
    152
    153
    154
    155
    156
    157
    158
    159
    160
    161
    162
    163
    /** différentielle par rapport à x**/
    
    float diffDx(const Image& in, int x, int y, int canal)
    {
      int width = in.getWidth();
    
      int px = x-1;
      int nx = x+1;
      if (px<0)
        px=0;
      if (nx>=width)
        nx=width-1;
    
      return((in.getPixel(nx,y,canal) - in.getPixel(px,y,canal))/2.0f);
    }
    
    float diffDxN(const Image & in, int x, int y, int canal)
    {
      int width = in.getWidth();
      if (x+1>=width)
        return 0;
    
      return(in.getPixel(x+1, y, canal) - in.getPixel(x,y,canal));
    }
    
    float diffDxP(const Image & in, int x, int y, int canal)
    {
      if(x-1<0)
        return 0;
      return(in.getPixel(x, y, canal) - in.getPixel(x-1,y,canal));
    }
    
    
    /** différentiel par rapport à y **/
    
    float diffDy(const Image& in, int x, int y, int canal)
    {
      int height = in.getHeight();
    
      int py = y-1;
      int ny = y+1;
    
      if (py<0)
        py=0;
      if (ny>=height)
        ny=height-1;
    
      return (in.getPixel(x,ny,canal) - in.getPixel(x,py,canal))/2.0f;
    }
    
    float diffDyN(const Image & in, int x, int y, int canal)
    {
      int height = in.getHeight();
      if (y+1>=height)
        return 0;
    
      return(in.getPixel(x, y+1, canal) - in.getPixel(x,y,canal));
    }
    
    float diffDyP(const Image & in, int x, int y, int canal)
    {
      if(y-1<0)
        return 0;
      return(in.getPixel(x, y, canal) - in.getPixel(x,y-1,canal));
    }
    
    /** norme du gradient **/
    
    /** version minmod **/
    float diffGradNminmod(const Image& in, int x, int y, int canal)
    {
      float dmoins= diffDyP(in, x, y, canal);
      float dplus= diffDyN(in, x, y, canal);
    
      fy= Math::minmod( dmoins, dplus );
    
      dplus= diffDxP(in, x, y, canal);
      dmoins= diffDxN(in,  x, y, canal);
    
      fx= Math::minmod( dmoins, dplus );
    
      return sqrt(fx² + fy²);
    }
    
    /** gradient simple **/
    
    float diffGradN(const Image& in, int x, int y, int canal)
    {
      return sqrt(diffDx(in, x, y, canal)² + diffDy(in, x, y, canal)²);
    }
    
    
    /** différentielle selon xx, xy et yy **/
    
    
    float diffDxx(const Image& in, int x, int y, int canal)
    {
      int width = in.getWidth();
    
      int px = x-1;
      int nx = x+1;
      if (px<0)
        px=0;
      if (nx>=width)
        nx=width-1;
    return (in.getPixel(nx,y,canal) + in.getPixel(px,y,canal) -2.0f * in.getPixel(x,y, canal));
    }
    
    
    float diffDyy(const Image& in, int x, int y, int canal)
    {
      int height = in.getHeight();
    
      int py = y-1;
      int ny = y+1;
    
      if (py<0)
        py=0;
      if (ny>=height)
        ny=height-1;
    
      return (in.getPixel(x,ny,canal) + in.getPixel(x,py,canal) -2.0f * in.getPixel(x,y, canal));
    }
    
    float diffDxy(const Image& in, int x, int y, int canal)
    {
      int width = in.getWidth();
      int height = in.getHeight();
    
      int px = x-1;
      int nx = x+1;
      int py = y-1;
      int ny = y+1;
      if (px<0)
        px=0;
      if (nx>=width)
        nx=width-1;
      if (py<0)
        py=0;
      if (ny>=height)
        ny=height-1;
    
      return 0.25 * (in.getPixel(nx, ny, canal) + in.getPixel(px, py, canal)
                    - in.getPixel(px, ny, canal) - in.getPixel(nx, py, canal));
    }
    
    
    /** différentielle suivant le gradient **/
    
    float Millie::diffDnn(const Image& in, int x, int y, int canal)
    {
      float normecarre = diffGradN(in, x,y, canal)²;
    
      float dx = diffDx(in, x, y, canal);
      float dy = diffDy(in, x, y, canal);
      float dxx = diffDxx(in, x, y, canal);
      float dyy = diffDyy(in, x, y, canal);
      float dxy = diffDxy(in, x, y, canal);
      float temp = dx² * dxx + dy² * dyy
                   + 2.0 * dx * dy * dxy;
      return (temp/(normecarre+0.001));
    }

    Avec la fonction minmod qui est définie par (il existe d'autres formulations) :

    Code :
    1
    2
    3
    4
    5
    6
    7
         float minmod(float a, float b)
          {
            if(a*b>0)
               return min(abs(a),abs(b));
            else
            return 0;
    }
    Et si vous saviez pas, sign est défini par :
    Code :
    1
    2
    3
    4
    5
    6
    7
    8
      float sign(float a)
          {
            if(a<0)
              return -1.0f;
            if(abs(a)==0.00) /*possibilité de mettre <0.01*/
              return 0.0f;
            return 1.0f;
          }

    Et l'algorithme s'écrit donc :

    Code :
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    /* @param out l'image de sortie
     * @param in l'image d'entrée
     * @param dt discrétisation temporelle (souvent proche de 0.1
     * @param nbiter le nombre d'itération
     *
     * erreur si dt<0 ou si nbiter<=0 ou si le nombre de composantes de out et de 
     * in ne sont pas pareil
     */
    void miniShock(Image& out, const Image& in,  float dt, int nbiter)
    {
      if(out.getNumComponents() != in.getNumComponents())
        Erreur
      if(dt<0)
        Erreur
      if(nbiter<=0)
        Erreur
    
      Image itempo(in.getWidth(), in.getHeight(), in.getNumComponents());
    
      out = in;
      itempo = in;
    
      int n;
      int i,j;
    
      for(n=0; n<nbiter;n++)
      {
        for(int canal = 0; canal< in.getNumComponents(); canal++)
          for(j=0;j<in.getHeight();j++)
            for(i=0;i<in.getWidth();i++)
            {
    
              float dnn = diffDnn(out, i, j, canal);
    
              float second_order = abs(diffDx(out, i, j, canal)) + abs(diffDy(out, i, j, canal));
    
              if(second_order < 1.0f) /*possibilité de mettre ==0*/
                dnn = diffDxx(out, i, j, canal);
    
              float grad = diffGradNminmod(out, i, j, canal);
    
              float value =   dt*(-sign(dnn) *  grad ) + out.getPixel(i,j, canal);
    
              itempo.setPixel(i, j, canal, value);
            }
    
        out = itempo;
      }
    
    
    }
    Voici avec : nbiter = 100 et dt = 0.1


    Lenna :

    Correction :




    Ca, c'est une version avec l'application d'une gaussienne (rayon 2, sigma = 10) avec Inn (je n'ai pas noté l'algorithme) : c'est la méthode d'Alvarez et de Mazorra
    Je ne répondrai à aucune question technique en privé

  2. #2
    Invité régulier
    Profil pro Boo Luck
    Inscrit en
    novembre 2009
    Messages
    6
    Détails du profil
    Informations personnelles :
    Nom : Boo Luck

    Informations forums :
    Inscription : novembre 2009
    Messages : 6
    Points : 7
    Points
    7

    Par défaut Et pour les gros floue

    Voici la version matlab. (et optimise pour)

    Code :
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    function i=miniShock_perf1( i, dt, nbiter)
        if(dt<0)
            warning('dt negatif')
        end
        if(nbiter<=0)
            warning('nbiter negatif')
        end
        tx=size(i,1);
        ty=size(i,2);
        maskup=ones(tx,ty);
        iup=ones(tx+2,ty+2);
        h=waitbar(0,'Init');
        for(n=0:nbiter)
            for(c=1:size(i,3))
            
                iup=[[i(1,1,c),i(1,:,c),i(1,end,c)];[i(:,1,c),i(:,:,c),i(:,end,c)];[i(end,1,c),i(end,:,c),i(end,end,c)]];
    
                ioxoy=iup(2:end-1,2:end-1);
                ipxoy=iup(2:end-1,3:end);
                inxoy=iup(2:end-1,1:end-2);
                ioxpy=iup(3:end,2:end-1);
                ioxny=iup(1:end-2,2:end-1);
                ipxpy=iup(3:end,3:end);
                inxpy=iup(3:end,1:end-2);
                inxny=iup(1:end-2,1:end-2);
                ipxny=iup(1:end-2,3:end);
    
    
                dx =(ipxoy-inxoy)/2;
                dxn=(ipxoy-ioxoy);
                dxp=(ioxoy-inxoy);
                dy =(ioxpy-ioxny)/2;
                dyn=(ioxpy-ioxoy);
                dyp=(ioxoy-ioxny);
                dxy=(ipxpy+inxny-ipxny-inxpy)/4;
                dxx=(ipxoy+inxoy-2*ioxoy);
                dyy=(ioxpy+ioxny-2*ioxoy);
    
    
                gradn=sqrt(dxx.*dxx+dyy.*dyy);
                dnn=(dx.*dx.*dxx+dy.*dy.*dyy+2*dx.*dy.*dxy)./(gradn.*gradn+0.001);
    
                minmodx=((dxn.*dxp)>0).*min(dxn,dxp);
                minmody=((dyn.*dyp)>0).*min(dyn,dyp);
                gradminmod=sqrt(minmodx.*minmodx+minmody.*minmody);
    
    
                second_order = abs(dx) + abs(dy);
                dnn=(second_order<1).*dxx+(second_order>=1).*dnn;
                i(:,:,c) = dt*(-sign(dnn).*gradminmod)+i(:,:,c);
                i(:,:,c) = max(min(i(:,:,c),255*maskup),0*maskup); %seuil
            end
    
            waitbar(0,h,['iter:',num2str(n)])
            imshow(i/max(max(max(i))))
            drawnow;
            pause(0.1)
        end
        close (h)
    end
    Il faut avouer que pour les floues de 1 a 5 pixels ca marche plutot bien, en tout cas mieux que photoshop.

    Par contre mon floue est plus de l'ordre de 400pxl ... oui c'est un jolie floue. Mais vu que l'image fait 18Mpxl (et qu'il y a tres peu de bruit (100ISO) : quasiement que de la quntificqtion JPEG), il doit bien avoir un moyen d'avoir qqc de potable avec 1 ou 2 Mpixl?
    (L'origine du floue est un defaut de Mise au point)

    J'ai cherche sur google les filtres choc, mais j'en ai pas trouve beacoups (et je pense pas que ca puisse s'appliquer å mon cas assez "extreme"). J'ai trouve 2 logiciels payant (specialises dans le defloutage) que j'ai essaye en demo et les resultats sont vraiment navrant (pire que photoshop ou l'algo ci dessus).

    D'ailleur c'est etonnant de voir sur le web autant de gens qui se posent ce probleme et si peu de reponses commercial.

  3. #3
    Rédacteur/Modérateur
    Avatar de pseudocode
    Homme Profil pro Xavier Philippeau
    Architecte système
    Inscrit en
    décembre 2006
    Messages
    9 962
    Détails du profil
    Informations personnelles :
    Nom : Homme Xavier Philippeau
    Âge : 41
    Localisation : France, Hérault (Languedoc Roussillon)

    Informations professionnelles :
    Activité : Architecte système
    Secteur : Industrie

    Informations forums :
    Inscription : décembre 2006
    Messages : 9 962
    Points : 15 080
    Points
    15 080

    Par défaut

    Citation Envoyé par factis Voir le message
    Par contre mon floue est plus de l'ordre de 400pxl ... oui c'est un jolie floue. Mais vu que l'image fait 18Mpxl (et qu'il y a tres peu de bruit (100ISO) : quasiement que de la quntificqtion JPEG), il doit bien avoir un moyen d'avoir qqc de potable avec 1 ou 2 Mpixl?
    (L'origine du floue est un defaut de Mise au point)
    Il faut une approche multi-echelle : construire une pyramide d'image a des tailles différentes, faire le défloutage sur l'image la plus petite, reporter les modifs dans l'image de taille juste supérieure... et recommencer.
    ALGORITHME (n.m.): Méthode complexe de résolution d'un problème simple.

  4. #4
    Invité régulier
    Profil pro Boo Luck
    Inscrit en
    novembre 2009
    Messages
    6
    Détails du profil
    Informations personnelles :
    Nom : Boo Luck

    Informations forums :
    Inscription : novembre 2009
    Messages : 6
    Points : 7
    Points
    7

    Par défaut

    Bonne idée, mais il risque pas d'avoir des gros bloc un peu trop visibles lors de l'assemblage?

    Si on décompose , on traite (donc des contours plus net), et on désassembles, on n'a que des contour plus net: on perd beaucoup de détail.

    Pour avoir une idée du problème:

    (ok, en fait y a l'air d'avoir un peu plus de bruit que prévu, mais y a toujours la masse de pixels)

    Si on souhaite appliquer une déconvolution, quelle est la forme de la PSF? (ça doit être un problème classique les floues de Mise au point?)

  5. #5
    Rédacteur/Modérateur
    Avatar de pseudocode
    Homme Profil pro Xavier Philippeau
    Architecte système
    Inscrit en
    décembre 2006
    Messages
    9 962
    Détails du profil
    Informations personnelles :
    Nom : Homme Xavier Philippeau
    Âge : 41
    Localisation : France, Hérault (Languedoc Roussillon)

    Informations professionnelles :
    Activité : Architecte système
    Secteur : Industrie

    Informations forums :
    Inscription : décembre 2006
    Messages : 9 962
    Points : 15 080
    Points
    15 080

    Par défaut

    Ah oui, c'est quand meme un peu flou. La je pense que c'est un peu trop pour un filtre de choc. Il faut passer par des algos de déconvolution et donc effectivement estimer la PSF. Usuellement la PSF pour du flou est modélisée par une gaussienne, mais il y a aussi des techniques pour avoir une meilleure estimation:

    google OUT OF FOCUS BLUR PSF
    ALGORITHME (n.m.): Méthode complexe de résolution d'un problème simple.

Liens sociaux

Règles de messages

  • Vous ne pouvez pas créer de nouvelles discussions
  • Vous ne pouvez pas envoyer des réponses
  • Vous ne pouvez pas envoyer des pièces jointes
  • Vous ne pouvez pas modifier vos messages
  •