add files

LICENSE

+Apache License 2.0
+
+But please contact me if you want to use this software in a commercial application. 
+
+Note that Apache License 2.0 asks to include a copyright notice if you use this software. 

PB2_mex.cpp

+#include "matrix.h"
+#include "mex.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+/*************************
+ *
+  clc
+  mex PB2_mex.cpp -output PB2
+ *
+ */
+
+double max(double a, double b){
+    if(a>b){
+        return a;
+    }
+    else{
+        return b;
+    }
+}
+
+double min(double a, double b){
+    if(a>b){
+        return b;
+    }
+    else{
+        return a;
+    }
+}
+
+void projection_6(double * entree, double *sortie){
+// Input:3x2 matrix A in double variables a11.. . a32
+// Output : projection \Pi(A) in double variables a11.. . a32
+// Compute D= AˆT A
+
+  double a11= entree[0];
+  double a12= entree[1];
+  double a21= entree[2];
+  double a22= entree[3];
+  double a31= entree[4];
+  double a32= entree[5];
+
+  double d11= a11 * a11 + a21 * a21 + a31 * a31 ;
+  double d12= a12 * a11 + a22 * a21 + a32 * a31 ;
+  double d22= a12 * a12 + a22 * a22 + a32 * a32 ;
+
+
+  // Compute Eigenvalues ( lmax , lmin )
+  //=squares of singular values ( smax , smin )
+
+  double trace = d11 + d22 ;
+  double det = d11 * d22 - d12 * d12 ;
+  double d = sqrt( max( 0.0 , 0.25 * trace * trace - det ) ) ;
+  double lmax = max( 0.0 , 0.5 * trace + d ) ;
+  double lmin = max( 0.0 , 0.5 * trace - d ) ;
+  double smax = sqrt ( lmax ) ;
+  double smin = sqrt ( lmin ) ;
+  // Check whether matrix lies outside K and need  o be projed
+  if ( smax + smin > 1.0 ) {
+    // Compute orthonormalsy st emofEigenvectors , such t h a t
+    // ( v11 , v21 ) belongs to lmax and ( v12 , v22 ) belongs to lmin
+    double v11 , v12 , v21 , v22 ;
+    if ( d12 == 0.0 ) {
+      if ( d11 >= d22 ) { v11 = 1.0 ; v21 = 0.0 ; v12 = 0.0 ; v22 = 1.0 ; }
+      else
+	{ v11 = 0.0 ; v21 = 1.0 ; v12 = 1.0 ; v22 = 0.0 ; }
+    }
+    else {
+      v11 = lmax - d22 ; v21 = d12 ;
+      double l1 = hypot( v11 , v21 ) ;
+      v11 /= l1 ; v21 /= l1 ;
+      v12 = lmin - d22 ; v22 = d12 ;
+      double l2 = hypot( v12 , v22 ) ;
+      v12 /= l2 ; v22 /= l2 ;
+    }
+    // Project ( smax , smin ) toline ( 0 , 1 ) + tau * (1 , -1) , 0<=tau <=1.
+    double tau = 0.5 * ( smax - smin + 1.0 ) ;
+    double s1 = min ( 1.0 , tau ) ;
+    double s2 = 1.0 - s1 ;
+    // Compute \ Sigmaˆ+ * \ Sigma p
+    s1 /= smax ;
+    s2 = ( smin > 0.0 ) ? s2 / smin : 0.0 ;
+    
+    // Compute
+    double t11= s1 * v11 * v11 + s2 * v12 * v12 ;
+    double t12= s1 * v11 * v21 + s2 * v12 * v22 ;
+    double t21= s1 * v21 * v11 + s2 * v22 * v12 ;
+    double t22= s1 * v21 * v21 + s2 * v22 * v22 ;
+    
+    // Result
+    a11 = a11* t11 + a12 *t21 ;
+    a12 = a11* t12 + a12 *t22 ;
+    a21 = a21* t11 + a22 *t21 ;
+    a22 = a21* t12 + a22 *t22 ;
+    a31 = a31* t11 + a32 *t21 ;
+    a32 = a31* t12 + a32 *t22 ;
+
+    sortie[0]= a11;
+    sortie[1]= a12;
+    sortie[2]= a21;
+    sortie[3]= a22;
+    sortie[4]= a31;
+    sortie[5]= a32;
+
+  }
+}
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+    /* Declarations.*/
+    double *Wl;
+    const double *Wr;
+    int height, width, ndim, longueur, i, pix, itterations;
+    const mwSize *dims;
+    double lambda, normi, sigma, tau, theta;
+    
+    /* Verification du nombre d'arguments.*/
+    if(nrhs!=1) {
+        mexErrMsgIdAndTxt("MyToolbox:PW_Total:nrhs","One inputs required.");
+    }
+    if(nlhs!=1) {
+        mexErrMsgIdAndTxt("MyToolbox:PW_Total:nlhs","One output required.");
+    }
+    
+    /*gestion entrees-sorties.*/
+    ndim= mxGetNumberOfDimensions(prhs[0]);
+    dims= mxGetDimensions(prhs[0]);
+    
+    height = dims[0];
+    width = dims[1];
+    int nb_chanel_dir= dims[2];
+    
+    
+    Wr = (const double*)(mxGetPr(prhs[0]));
+    
+    const int n_dim_result= 3;
+    mwSize dims_result[n_dim_result];
+    dims_result[0]= height;
+    dims_result[1]= width;
+    dims_result[2]= nb_chanel_dir;
+    
+    plhs[0]= mxCreateNumericArray(n_dim_result, dims_result, mxDOUBLE_CLASS, mxREAL);
+    Wl= (double*)mxGetPr(plhs[0]);
+    
+    double *x= (double*)malloc(nb_chanel_dir* sizeof(double));
+    double *y= (double*)malloc(nb_chanel_dir* sizeof(double));
+    
+    /* Copie des vecteurs pour traitement*/
+    for(int lin=0; lin< height; lin++){
+        for(int col=0; col< width; col++){
+            for(int w=0; w< nb_chanel_dir; w++){
+                x[w]= Wr[lin + height * col + w * width * height];
+            }
+            
+            projection_6(x,y);
+            
+            for(int w=0; w< nb_chanel_dir; w++){
+                Wl[lin + height * col + w * width * height]= y[w];
+            }
+        }
+    }
+    
+    free(x);
+    free(y);
+    
+    return;
+}
+

divergence_mex.cpp

+/*void divergence(double *gradient_x, double *gradient_y, double * result, int height, int width )
+{
+    int col, lin;
+    double x_1;
+    double x_2;
+    int N;
+    N=width*height;
+    for(lin=0; lin<height; lin++)
+    {
+        for(col=0; col<width; col++)
+        {
+            if(lin==0)
+            {
+                x_1=gradient_x[lin+height*col];
+            }
+            else if (lin==(height-1))
+            {
+                x_1=-gradient_x[lin+height*col-1];
+            }
+            else
+            {
+                x_1=gradient_x[lin+height*col]
+                        -gradient_x[lin+height*col-1];
+            }
+            
+            if (col==0)
+            {
+                x_2=gradient_y[lin+height*col];
+            }
+            else if (col==(width-1))
+            {
+                x_2=-gradient_y[lin+height*(col-1)];
+            }
+            else
+            {
+                x_2=gradient_y[lin+height*col]
+                        -gradient_y[lin+height*(col-1)];
+            }
+            result[lin+height*col]= x_1+x_2;
+        }
+    }
+}
+*/
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "matrix.h"
+#include "mex.h"
+
+/*************************************************************************
+ 
+ AUTEUR : Fabien PIERRE.
+ DATE   : 15 avril 2013.
+ 
+ 
+ Compilation :
+ 
+  % linux
+  clc
+  mex CFLAGS="\$CFLAGS -std=c99" divergence_mex.c -output divergence_mex
+ 
+  % windows
+  clc
+  mex  PG_poly_mex3.c -output PG_poly_mex
+ 
+ 
+ ************************************************************************/
+/*void gradient(double * image_initiale, double *gradient_x, double *gradient_y, int height, int width )
+{
+    int N;
+    int col, lin;
+    N=width*height;
+    for(col=0; col<width; col++)
+    {
+        for(lin=0; lin<height; lin++)
+        {
+            if(lin==(height-1))
+            {
+                gradient_x[lin+col*height]=0;
+            }
+            else
+            {
+                gradient_x[lin+col*height]=
+                        image_initiale[lin+(col)*height+1]
+                        -image_initiale[lin+col*height];
+            }
+            if((col==(width-1)))
+            {
+                gradient_y[lin+col*height]= 0;
+            }
+            else
+            {
+                gradient_y[lin+col*height]=
+                        image_initiale[lin+(col+1)*height]
+                        -image_initiale[lin+col*height];
+            }
+        }
+    }
+}*/
+
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+    /* Declarations.*/
+    double *Wl;
+    const double *Wr;
+    int height, width, ndim, longueur, i, pix, itterations;
+    const mwSize *dims;
+    double lambda, normi, sigma, tau, theta;
+    
+    /* V�rification du nombre d'arguments.*/
+    if(nrhs!=1) {
+        mexErrMsgIdAndTxt("MyToolbox:PW_Total:nrhs","One inputs required.");
+    }
+    if(nlhs!=1) {
+        mexErrMsgIdAndTxt("MyToolbox:PW_Total:nlhs","One output required.");
+    }
+    
+    /*gestion entrees-sorties.*/
+    ndim= mxGetNumberOfDimensions(prhs[0]);
+    dims= mxGetDimensions(prhs[0]);
+    
+    height = dims[0];
+    width = dims[1];
+    int nb_chanel_dir= dims[2];
+    
+    
+    Wr = (const double*)(mxGetPr(prhs[0]));
+    
+    const int n_dim_result= 3;
+    mwSize dims_result[n_dim_result];
+    dims_result[0]= height;
+    dims_result[1]= width;
+    dims_result[2]= nb_chanel_dir/2;
+    
+    
+    
+    plhs[0]= mxCreateNumericArray(n_dim_result, dims_result, mxDOUBLE_CLASS, mxREAL);
+    Wl= (double*)mxGetPr(plhs[0]);
+    
+        /* Initialisation de l'algorithme.*/
+  /*  for(int pix = 0; pix < height*width*nb_chanel_dir ; pix++)
+    {
+        Wl[pix] = Wr[pix];
+    }
+    */
+    int nb_color= nb_chanel_dir/2;       
+    
+    double x_1;
+    double x_2;
+    
+    for(int c=0; c < nb_color; c++){
+        
+    for(int lin=0; lin<height; lin++)
+    {
+        for(int col=0; col<width; col++)
+        {
+            if(lin==0)
+            {
+                x_1=Wr[lin+height*col + (2*c+1)*height*width];
+            }
+            else if (lin==(height-1))
+            {
+                x_1=-Wr[lin+height*col-1+ (2*c+1)*height*width ];
+            }
+            else
+            {
+                x_1=Wr[lin+height*col+ (2*c+1)*height*width]
+                        -Wr[lin+height*col-1+ (2*c+1)*height*width];
+            }
+            
+            if (col==0)
+            {
+                x_2=Wr[lin+height*col+ (2*c)*height*width];
+            }
+            else if (col==(width-1))
+            {
+                x_2=-Wr[lin+height*(col-1)+ (2*c)*height*width];
+            }
+            else
+            {
+                x_2=Wr[lin+height*col+ (2*c)*height*width]
+                        -Wr[lin+height*(col-1)+ (2*c)*height*width];
+            }
+            Wl[lin+height*col + (c)*height*width]= x_1+x_2;
+        }
+    }
+}
+    
+    return;
+}
+

gradient_mex.cpp

+
+#include "matrix.h"
+#include "mex.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+/*************************************************************************
+ 
+ AUTEUR : Fabien PIERRE.
+ DATE   : 15 avril 2013.
+ 
+ 
+ Compilation :
+ 
+  % linux
+  clc
+  mex CFLAGS="\$CFLAGS -std=c99" gradient_mex.c -output gradient_mex
+ 
+  % windows
+  clc
+  mex  PG_poly_mex3.c -output PG_poly_mex
+ 
+ 
+ ************************************************************************/
+void gradient(double * image_initiale, double *gradient_x, double *gradient_y, int height, int width )
+{
+    int N;
+    int col, lin;
+    N=width*height;
+    for(col=0; col<width; col++)
+    {
+        for(lin=0; lin<height; lin++)
+        {
+            if(lin==(height-1))
+            {
+                gradient_x[lin+col*height]=0;
+            }
+            else
+            {
+                gradient_x[lin+col*height]=
+                        image_initiale[lin+(col)*height+1]
+                        -image_initiale[lin+col*height];
+            }
+            if((col==(width-1)))
+            {
+                gradient_y[lin+col*height]= 0;
+            }
+            else
+            {
+                gradient_y[lin+col*height]=
+                        image_initiale[lin+(col+1)*height]
+                        -image_initiale[lin+col*height];
+            }
+        }
+    }
+}
+
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+    /* Declarations.*/
+    double *Wl;
+    const double *Wr;
+    int height, width, ndim, longueur, i, pix, itterations;
+    const mwSize *dims;
+    double lambda, normi, sigma, tau, theta;
+    
+    /* V�rification du nombre d'arguments.*/
+    if(nrhs!=1) {
+        mexErrMsgIdAndTxt("MyToolbox:PW_Total:nrhs","One inputs required.");
+    }
+    if(nlhs!=1) {
+        mexErrMsgIdAndTxt("MyToolbox:PW_Total:nlhs","One output required.");
+    }
+    
+    /*gestion entrees-sorties.*/
+    ndim= mxGetNumberOfDimensions(prhs[0]);
+    dims= mxGetDimensions(prhs[0]);
+    
+    height = dims[0];
+    width = dims[1];
+    int nb_chanel_dir= dims[2];
+    
+    
+    Wr = (const double*)(mxGetPr(prhs[0]));
+    
+    const int n_dim_result= 3;
+    mwSize dims_result[n_dim_result];
+    dims_result[0]= height;
+    dims_result[1]= width;
+    dims_result[2]= 2*nb_chanel_dir;
+    
+    
+    
+    plhs[0]= mxCreateNumericArray(n_dim_result, dims_result, mxDOUBLE_CLASS, mxREAL);
+    Wl= (double*)mxGetPr(plhs[0]);
+    
+        /* Initialisation de l'algorithme.*/
+    /*for(int pix = 0; pix < height*width*nb_chanel_dir ; pix++)
+    {
+        Wl[pix] = Wr[pix];
+    }*/
+    
+    int nb_color= nb_chanel_dir;       
+    
+    for(int c=0; c < nb_color; c++){
+        
+        
+        for(int col=0; col<width; col++)
+        {
+            for(int lin=0; lin<height; lin++)
+            {
+                if(lin==(height-1))
+                {
+                    Wl[lin+col*height + (2*c+1)*height*width]=0;
+                }
+                else
+                {
+                    Wl[lin+col*height + (2*c+1)*height*width]=
+                            Wr[lin+col*height+1 + c*height*width]
+                            -Wr[lin+col*height  + c*height*width];
+                }
+                if((col==(width-1)))
+                {
+                    Wl[lin+col*height + (2*c)*height*width]= 0;
+                }
+                else
+                {
+                    Wl[lin+col*height + (2*c)*height*width]=
+                            Wr[lin+(col+1)*height +c *height*width]
+                            -Wr[lin+col*height + c*height*width];
+                }
+            }
+        }
+        
+       
+        
+        
+    }
+    
+    return;
+}

main.m

+%% Compilation
+
+clc
+mex PB2_mex.cpp -output PB2
+mex divergence_mex.cpp -output divergence
+mex gradient_mex.cpp -output gradient
+
+%% Script de test (equation 4.1)
+
+clc, clear all, close all
+U = double(imread('lena.png'));
+
+figure
+imagesc(uint8(U))
+title('initialisation de U')
+
+[height, width,~] = size(U);
+Z = zeros(height, width, 6); % variable duale
+
+% Variables
+sigma = 1;
+tau = 1;
+ITT = 2500;
+lambda = 0.02;
+
+Donnees = double(U);
+
+for itt = 1:ITT
+    Z = PB2(Z + sigma * gradient(U));
+    U = (U + tau * ( divergence(Z) + lambda  * Donnees)) / ...
+             ( 1 + tau * lambda);
+end
+
+figure
+imagesc(uint8(U))
+title('U convergence')
+axis off
+
+imwrite(uint8(U), 'resultat.png')
+