fixes all matlab issues

6 files changed, 150 insertions, 17 deletions

diff --git a/src/Matlab/mex_compile/compileCPU_mex_Linux.m b/src/Matlab/mex_compile/compileCPU_mex_Linux.m
index 19fb1a5..2ed7ea6 100644
--- a/src/Matlab/mex_compile/compileCPU_mex_Linux.m
+++ b/src/Matlab/mex_compile/compileCPU_mex_Linux.m
@@ -77,5 +77,5 @@ delete PatchSelect_core* Nonlocal_TV_core*
 delete Diffusion_Inpaint_core* NonlocalMarching_Inpaint_core*
 fprintf('%s \n', '<<<<<<< Regularisers successfully compiled! >>>>>>>');
 
-pathA2 = sprintf(['..' fsep '..' fsep '..' fsep '..' fsep 'demos'], 1i);
+pathA2 = sprintf(['..' fsep '..' fsep '..' fsep '..' fsep 'demos' fsep 'Matlab_demos'], 1i);
 cd(pathA2);
\ No newline at end of file
diff --git a/src/Matlab/mex_compile/compileCPU_mex_WINDOWS.m b/src/Matlab/mex_compile/compileCPU_mex_WINDOWS.m
index 3a9e2af..7d0917f 100644
--- a/src/Matlab/mex_compile/compileCPU_mex_WINDOWS.m
+++ b/src/Matlab/mex_compile/compileCPU_mex_WINDOWS.m
@@ -7,7 +7,11 @@
 % Here I present two ways how software can be compiled, if you have some
 % other suggestions/remarks please contact me at dkazanc@hotmail.com 
 % >>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-
+% cuda related messgaes can be  solved in Linux:
+%export LD_LIBRARY_PATH="$LD_LIBRARY_PATH:/usr/local/cuda/lib64"
+%export CUDA_HOME=/usr/local/cuda
+%export LD_PRELOAD=/home/algol/anaconda3/lib/libstdc++.so.6.0.24
+% >>
 fsep = '/';
 
 pathcopyFrom = sprintf(['..' fsep '..' fsep 'Core' fsep 'regularisers_CPU'], 1i);
diff --git a/src/Matlab/mex_compile/compileGPU_mex.m b/src/Matlab/mex_compile/compileGPU_mex.m
index 7e15233..56fcd38 100644
--- a/src/Matlab/mex_compile/compileGPU_mex.m
+++ b/src/Matlab/mex_compile/compileGPU_mex.m
@@ -66,10 +66,14 @@ fprintf('%s \n', 'Compiling ROF-LLT...');
 mex -g -I/usr/local/cuda-10.0/include -L/usr/local/cuda-10.0/lib64 -lcudart -lcufft -lmwgpu LLT_ROF_GPU.cpp LLT_ROF_GPU_core.o
 movefile('LLT_ROF_GPU.mex*',Pathmove);
 
+fprintf('%s \n', 'Compiling PatchSelect...');
+!/usr/local/cuda/bin/nvcc -O0 -c PatchSelect_GPU_core.cu -Xcompiler -fPIC -I~/SOFT/MATLAB9/extern/include/
+mex -g -I/usr/local/cuda-10.0/include -L/usr/local/cuda-10.0/lib64 -lcudart -lcufft -lmwgpu PatchSelect_GPU.cpp PatchSelect_GPU_core.o
+movefile('PatchSelect_GPU.mex*',Pathmove);
 
 delete TV_ROF_GPU_core* TV_FGP_GPU_core* TV_SB_GPU_core* dTV_FGP_GPU_core* NonlDiff_GPU_core* Diffus_4thO_GPU_core* TGV_GPU_core* LLT_ROF_GPU_core* CCPiDefines.h
-delete PatchSelect_core* Nonlocal_TV_core* shared.h
+delete PatchSelect_GPU_core* Nonlocal_TV_core* shared.h
 fprintf('%s \n', 'All successfully compiled!');
 
-pathA2 = sprintf(['..' fsep '..' fsep '..' fsep '..' fsep 'demos'], 1i);
+pathA2 = sprintf(['..' fsep '..' fsep '..' fsep '..' fsep 'demos' fsep 'Matlab_demos'], 1i);
 cd(pathA2);
diff --git a/src/Matlab/mex_compile/regularisers_CPU/ROF_TV.c b/src/Matlab/mex_compile/regularisers_CPU/ROF_TV.c
index ffe7b91..d6bdaa6 100644
--- a/src/Matlab/mex_compile/regularisers_CPU/ROF_TV.c
+++ b/src/Matlab/mex_compile/regularisers_CPU/ROF_TV.c
@@ -26,7 +26,7 @@
  * 
  * Input Parameters:
  * 1. Noisy image/volume [REQUIRED]
- * 2. lambda - regularization parameter [REQUIRED]
+ * 2. lambda - regularization parameter [REQUIRED] scalar or the same size as 1
  * 3. Number of iterations, for explicit scheme >= 150 is recommended  [REQUIRED]
  * 4. tau - marching step for explicit scheme, ~1 is recommended [REQUIRED]
  * 5. eplsilon: tolerance constant [REQUIRED]
@@ -38,7 +38,7 @@
  * This function is based on the paper by
  * [1] Rudin, Osher, Fatemi, "Nonlinear Total Variation based noise removal algorithms"
  *
- * D. Kazantsev, 2016-18
+ * D. Kazantsev, 2016-19
  */
 
 void mexFunction(
@@ -49,18 +49,29 @@ void mexFunction(
     int number_of_dims, iter_numb;
     mwSize dimX, dimY, dimZ;
     const mwSize *dim_array_i;
-    float *Input, *Output=NULL, lambda, tau, epsil;    
+    float *Input, *Output=NULL, lambda_scalar, tau, epsil;        
     float *infovec=NULL;
+    float *lambda;      
     
     dim_array_i = mxGetDimensions(prhs[0]);
     number_of_dims = mxGetNumberOfDimensions(prhs[0]);
     
     /*Handling Matlab input data*/
     Input  = (float *) mxGetData(prhs[0]);
-    lambda =  (float) mxGetScalar(prhs[1]); /* regularization parameter */
+    /*special check to find the input scalar or an array*/
+    int mrows = mxGetM(prhs[1]);
+    int ncols = mxGetN(prhs[1]); 
+    if (mrows==1 && ncols==1) {        
+        lambda = (float*) calloc (1 ,sizeof(float));
+        lambda_scalar =  (float) mxGetScalar(prhs[1]); /* regularization parameter */        
+        lambda[0] = lambda_scalar;
+    }
+    else {
+        lambda =  (float *) mxGetData(prhs[1]); /* regularization parameter */    
+    }
     iter_numb =  (int) mxGetScalar(prhs[2]); /* iterations number */
     tau =  (float) mxGetScalar(prhs[3]); /* marching step parameter */  
-    epsil = (float) mxGetScalar(prhs[4]); /* tolerance */  
+    epsil = (float) mxGetScalar(prhs[4]); /* tolerance */    
     
     if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
     if(nrhs != 5) mexErrMsgTxt("Four inputs reqired: Image(2D,3D), regularization parameter, iterations number,  marching step constant, tolerance");
@@ -80,6 +91,11 @@ void mexFunction(
     mwSize vecdim[1];
     vecdim[0] = 2;
     infovec = (float*)mxGetPr(plhs[1] = mxCreateNumericArray(1, vecdim, mxSINGLE_CLASS, mxREAL));
-     
-    TV_ROF_CPU_main(Input, Output, infovec, lambda, iter_numb, tau, epsil, dimX, dimY, dimZ);    
+    
+    if (mrows==1 && ncols==1) {
+    TV_ROF_CPU_main(Input, Output, infovec, lambda, 0, iter_numb, tau, epsil, dimX, dimY, dimZ);    
+    free(lambda);
+    }
+    else TV_ROF_CPU_main(Input, Output, infovec, lambda, 1, iter_numb, tau, epsil, dimX, dimY, dimZ);     
+        
 }
diff --git a/src/Matlab/mex_compile/regularisers_GPU/PatchSelect_GPU.cpp b/src/Matlab/mex_compile/regularisers_GPU/PatchSelect_GPU.cpp
new file mode 100644
index 0000000..13319fa
--- /dev/null
+++ b/src/Matlab/mex_compile/regularisers_GPU/PatchSelect_GPU.cpp
@@ -0,0 +1,94 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC and Diamond Light Source Ltd. 
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ * Copyright 2018 Diamond Light Source Ltd. 
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include "matrix.h"
+#include "mex.h"
+#include "PatchSelect_GPU_core.h"
+
+/* CUDA implementation of non-local weight pre-calculation for non-local priors
+ * Weights and associated indices are stored into pre-allocated arrays and passed
+ * to the regulariser
+ *
+ *
+ * Input Parameters:
+ * 1. 2D/3D grayscale image/volume
+ * 2. Searching window (half-size of the main bigger searching window, e.g. 11)
+ * 3. Similarity window (half-size of the patch window, e.g. 2)
+ * 4. The number of neighbours to take (the most prominent after sorting neighbours will be taken)
+ * 5. noise-related parameter to calculate non-local weights
+ *
+ * Output [2D]:
+ * 1. AR_i - indeces of i neighbours
+ * 2. AR_j - indeces of j neighbours
+ * 3. Weights_ij - associated weights
+ *
+ * Output [3D]:
+ * 1. AR_i - indeces of i neighbours
+ * 2. AR_j - indeces of j neighbours
+ * 3. AR_k - indeces of j neighbours
+ * 4. Weights_ijk - associated weights
+ */
+/**************************************************/
+void mexFunction(
+        int nlhs, mxArray *plhs[],
+        int nrhs, const mxArray *prhs[])
+{
+    int number_of_dims,  SearchWindow, SimilarWin, NumNeighb;
+    mwSize dimX, dimY, dimZ;
+    const mwSize *dim_array;
+    unsigned short *H_i=NULL, *H_j=NULL, *H_k=NULL;    
+    float *A, *Weights = NULL, h;
+    mwSize dim_array2[3]; /* for 2D data */
+    mwSize dim_array3[4]; /* for 3D data */
+    
+    dim_array = mxGetDimensions(prhs[0]);
+    number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+    
+    /*Handling Matlab input data*/
+    A  = (float *) mxGetData(prhs[0]); /* a 2D or 3D image/volume */
+    SearchWindow = (int) mxGetScalar(prhs[1]);    /* Large Searching window */
+    SimilarWin = (int) mxGetScalar(prhs[2]);    /* Similarity window (patch-search)*/
+    NumNeighb = (int) mxGetScalar(prhs[3]); /* the total number of neighbours to take */
+    h = (float) mxGetScalar(prhs[4]); /* NLM parameter */
+
+    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+    dim_array2[0] = dimX; dim_array2[1] = dimY; dim_array2[2] = NumNeighb;  /* 2D case */
+    dim_array3[0] = dimX; dim_array3[1] = dimY; dim_array3[2] = dimZ; dim_array3[3] = NumNeighb;  /* 3D case */
+    
+    /****************2D INPUT ***************/
+    if (number_of_dims == 2) {
+        dimZ = 0;               
+        H_i = (unsigned short*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array2, mxUINT16_CLASS, mxREAL));
+        H_j = (unsigned short*)mxGetPr(plhs[1] = mxCreateNumericArray(3, dim_array2, mxUINT16_CLASS, mxREAL));
+        Weights = (float*)mxGetPr(plhs[2] = mxCreateNumericArray(3, dim_array2, mxSINGLE_CLASS, mxREAL));
+        }
+    /****************3D INPUT ***************/
+    if (number_of_dims == 3) { 
+        /*
+        H_i = (unsigned short*)mxGetPr(plhs[0] = mxCreateNumericArray(4, dim_array3, mxUINT16_CLASS, mxREAL));
+        H_j = (unsigned short*)mxGetPr(plhs[1] = mxCreateNumericArray(4, dim_array3, mxUINT16_CLASS, mxREAL));
+        H_k = (unsigned short*)mxGetPr(plhs[2] = mxCreateNumericArray(4, dim_array3, mxUINT16_CLASS, mxREAL));
+        Weights = (float*)mxGetPr(plhs[3] = mxCreateNumericArray(4, dim_array3, mxSINGLE_CLASS, mxREAL));        
+        */
+        mexErrMsgTxt("3D version is not yet supported");
+    }
+    
+    PatchSelect_GPU_main(A, H_i, H_j, Weights, (long)(dimX), (long)(dimY), SearchWindow, SimilarWin, NumNeighb, h);         
+ }
diff --git a/src/Matlab/mex_compile/regularisers_GPU/ROF_TV_GPU.cpp b/src/Matlab/mex_compile/regularisers_GPU/ROF_TV_GPU.cpp
index d9b7e83..2cd4469 100644
--- a/src/Matlab/mex_compile/regularisers_GPU/ROF_TV_GPU.cpp
+++ b/src/Matlab/mex_compile/regularisers_GPU/ROF_TV_GPU.cpp
@@ -25,7 +25,7 @@
  * 
  * Input Parameters:
  * 1. Noisy image/volume [REQUIRED]
- * 2. lambda - regularization parameter [REQUIRED]
+ * 2. lambda - regularization parameter [REQUIRED], scalar or the same size as 1
  * 3. Number of iterations, for explicit scheme >= 150 is recommended  [REQUIRED]
  * 4. tau - marching step for explicit scheme, ~1 is recommended [REQUIRED]
  * 5. eplsilon: tolerance constant [REQUIRED]
@@ -37,7 +37,7 @@
  * This function is based on the paper by
  * [1] Rudin, Osher, Fatemi, "Nonlinear Total Variation based noise removal algorithms"
  *
- * D. Kazantsev, 2016-18
+ * D. Kazantsev, 2016-19
  */
 void mexFunction(
         int nlhs, mxArray *plhs[],
@@ -47,15 +47,26 @@ void mexFunction(
      int number_of_dims, iter_numb;
     mwSize dimX, dimY, dimZ;
     const mwSize *dim_array_i;
-    float *Input, *Output=NULL, lambda, tau, epsil;    
+    float *Input, *Output=NULL, lambda_scalar, tau, epsil;    
     float *infovec=NULL;
-    
+    float *lambda;  
+        
     dim_array_i = mxGetDimensions(prhs[0]);
     number_of_dims = mxGetNumberOfDimensions(prhs[0]);
     
     /*Handling Matlab input data*/
     Input  = (float *) mxGetData(prhs[0]);
-    lambda =  (float) mxGetScalar(prhs[1]); /* regularization parameter */
+    /*special check to find the input scalar or an array*/
+    int mrows = mxGetM(prhs[1]);
+    int ncols = mxGetN(prhs[1]); 
+    if (mrows==1 && ncols==1) {        
+        lambda = (float*) calloc (1 ,sizeof(float));
+        lambda_scalar =  (float) mxGetScalar(prhs[1]); /* regularization parameter */        
+        lambda[0] = lambda_scalar;
+    }
+    else {
+        lambda =  (float *) mxGetData(prhs[1]); /* regularization parameter */    
+    }
     iter_numb =  (int) mxGetScalar(prhs[2]); /* iterations number */
     tau =  (float) mxGetScalar(prhs[3]); /* marching step parameter */  
     epsil = (float) mxGetScalar(prhs[4]); /* tolerance */  
@@ -79,5 +90,9 @@ void mexFunction(
     vecdim[0] = 2;
     infovec = (float*)mxGetPr(plhs[1] = mxCreateNumericArray(1, vecdim, mxSINGLE_CLASS, mxREAL));
     
-    TV_ROF_GPU_main(Input, Output, infovec, lambda, iter_numb, tau, epsil, dimX, dimY, dimZ);    
+    if (mrows==1 && ncols==1) {
+    TV_ROF_GPU_main(Input, Output, infovec, lambda, 0, iter_numb, tau, epsil, dimX, dimY, dimZ);    
+    free(lambda);
+    }
+    else TV_ROF_GPU_main(Input, Output, infovec, lambda, 1, iter_numb, tau, epsil, dimX, dimY, dimZ);       
 }