Merge pull request #132 from vais-ral/nlmult

bilinear interpolation rescaling

3 files changed, 138 insertions, 57 deletions

diff --git a/src/Core/regularisers_CPU/utils.c b/src/Core/regularisers_CPU/utils.c
index 7a4e80b..5bb3a5c 100644
--- a/src/Core/regularisers_CPU/utils.c
+++ b/src/Core/regularisers_CPU/utils.c
@@ -1,21 +1,21 @@
 /*
-This work is part of the Core Imaging Library developed by
-Visual Analytics and Imaging System Group of the Science Technology
-Facilities Council, STFC
-
-Copyright 2017 Daniil Kazanteev
-Copyright 2017 Srikanth Nagella, Edoardo Pasca
-
-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.
-*/
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2017 Daniil Kazanteev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ *
+ * 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 "utils.h"
 #include <math.h>
@@ -23,19 +23,19 @@ limitations under the License.
 /* Copy Image (float) */
 float copyIm(float *A, float *U, long dimX, long dimY, long dimZ)
 {
-	long j;
+    long j;
 #pragma omp parallel for shared(A, U) private(j)
-	for (j = 0; j<dimX*dimY*dimZ; j++)  U[j] = A[j];
-	return *U;
+    for (j = 0; j<dimX*dimY*dimZ; j++)  U[j] = A[j];
+    return *U;
 }
 
 /* Copy Image */
 unsigned char copyIm_unchar(unsigned char *A, unsigned char *U, int dimX, int dimY, int dimZ)
 {
-	int j;
+    int j;
 #pragma omp parallel for shared(A, U) private(j)
-	for (j = 0; j<dimX*dimY*dimZ; j++)  U[j] = A[j];
-	return *U;
+    for (j = 0; j<dimX*dimY*dimZ; j++)  U[j] = A[j];
+    return *U;
 }
 
 /*Roll image symmetrically from top to bottom*/
@@ -59,59 +59,88 @@ float copyIm_roll(float *A, float *U, int dimX, int dimY, int roll_value, int sw
 
 /* function that calculates TV energy
  * type - 1:  2*lambda*min||\nabla u|| + ||u -u0||^2
- * type - 2:  2*lambda*min||\nabla u|| 
+ * type - 2:  2*lambda*min||\nabla u||
  * */
 float TV_energy2D(float *U, float *U0, float *E_val, float lambda, int type, int dimX, int dimY)
 {
-	int i, j, i1, j1, index;
-	float NOMx_2, NOMy_2, E_Grad=0.0f, E_Data=0.0f;
-	
-	/* first calculate \grad U_xy*/	
-        for(j=0; j<dimY; j++) {
-            for(i=0; i<dimX; i++) {
-				index = j*dimX+i;
-                /* boundary conditions */
-                i1 = i + 1; if (i == dimX-1) i1 = i;
-                j1 = j + 1; if (j == dimY-1) j1 = j;
-                
-                /* Forward differences */                
-                NOMx_2 = powf((float)(U[j1*dimX + i] - U[index]),2); /* x+ */
-                NOMy_2 = powf((float)(U[j*dimX + i1] - U[index]),2); /* y+ */
-                E_Grad += 2.0f*lambda*sqrtf((float)(NOMx_2) + (float)(NOMy_2)); /* gradient term energy */
-                E_Data += powf((float)(U[index]-U0[index]),2); /* fidelity term energy */
-			}
-		}
-		if (type == 1) E_val[0] = E_Grad + E_Data;
-		if (type == 2) E_val[0] = E_Grad;
-		return *E_val;
+    int i, j, i1, j1, index;
+    float NOMx_2, NOMy_2, E_Grad=0.0f, E_Data=0.0f;
+    
+    /* first calculate \grad U_xy*/
+    for(j=0; j<dimY; j++) {
+        for(i=0; i<dimX; i++) {
+            index = j*dimX+i;
+            /* boundary conditions */
+            i1 = i + 1; if (i == dimX-1) i1 = i;
+            j1 = j + 1; if (j == dimY-1) j1 = j;
+            
+            /* Forward differences */
+            NOMx_2 = powf((float)(U[j1*dimX + i] - U[index]),2); /* x+ */
+            NOMy_2 = powf((float)(U[j*dimX + i1] - U[index]),2); /* y+ */
+            E_Grad += 2.0f*lambda*sqrtf((float)(NOMx_2) + (float)(NOMy_2)); /* gradient term energy */
+            E_Data += powf((float)(U[index]-U0[index]),2); /* fidelity term energy */
+        }
+    }
+    if (type == 1) E_val[0] = E_Grad + E_Data;
+    if (type == 2) E_val[0] = E_Grad;
+    return *E_val;
 }
 
 float TV_energy3D(float *U, float *U0, float *E_val, float lambda, int type, int dimX, int dimY, int dimZ)
 {
-	long i, j, k, i1, j1, k1, index;
-	float NOMx_2, NOMy_2, NOMz_2, E_Grad=0.0f, E_Data=0.0f;
-	
-	/* first calculate \grad U_xy*/	
+    long i, j, k, i1, j1, k1, index;
+    float NOMx_2, NOMy_2, NOMz_2, E_Grad=0.0f, E_Data=0.0f;
+    
+    /* first calculate \grad U_xy*/
     for(j=0; j<(long)(dimY); j++) {
         for(i=0; i<(long)(dimX); i++) {
             for(k=0; k<(long)(dimZ); k++) {
-				index = (dimX*dimY)*k + j*dimX+i;
+                index = (dimX*dimY)*k + j*dimX+i;
                 /* boundary conditions */
                 i1 = i + 1; if (i == (long)(dimX-1)) i1 = i;
                 j1 = j + 1; if (j == (long)(dimY-1)) j1 = j;
                 k1 = k + 1; if (k == (long)(dimZ-1)) k1 = k;
                 
-                /* Forward differences */                
+                /* Forward differences */
                 NOMx_2 = powf((float)(U[(dimX*dimY)*k + j1*dimX+i] - U[index]),2); /* x+ */
                 NOMy_2 = powf((float)(U[(dimX*dimY)*k + j*dimX+i1] - U[index]),2); /* y+ */
                 NOMz_2 = powf((float)(U[(dimX*dimY)*k1 + j*dimX+i] - U[index]),2); /* z+ */
                 
                 E_Grad += 2.0f*lambda*sqrtf((float)(NOMx_2) + (float)(NOMy_2) + (float)(NOMz_2)); /* gradient term energy */
                 E_Data += (powf((float)(U[index]-U0[index]),2)); /* fidelity term energy */
-			}
-		}
-	}
-		if (type == 1) E_val[0] = E_Grad + E_Data;
-		if (type == 2) E_val[0] = E_Grad;
-		return *E_val;
+            }
+        }
+    }
+    if (type == 1) E_val[0] = E_Grad + E_Data;
+    if (type == 2) E_val[0] = E_Grad;
+    return *E_val;
+}
+
+/* Down-Up scaling of 2D images using bilinear interpolation */
+float Im_scale2D(float *Input, float *Scaled, int w, int h, int w2, int h2)
+{
+    int x, y, index, i, j;
+    float x_ratio = ((float)(w-1))/w2;
+    float y_ratio = ((float)(h-1))/h2;
+    float A, B, C, D, x_diff, y_diff, gray;
+    #pragma omp parallel for shared (Input, Scaled) private(x, y, index, A, B, C, D, x_diff, y_diff, gray)
+    for (j=0;j<w2;j++) {
+        for (i=0;i<h2;i++) {
+            x = (int)(x_ratio * j);
+            y = (int)(y_ratio * i);
+            x_diff = (x_ratio * j) - x;
+            y_diff = (y_ratio * i) - y;
+            index = y*w+x ;
+            
+            A = Input[index];
+            B = Input[index+1];
+            C = Input[index+w];
+            D = Input[index+w+1];
+            
+            gray = (float)(A*(1.0 - x_diff)*(1.0 - y_diff) +  B*(x_diff)*(1.0 - y_diff) +
+                    C*(y_diff)*(1.0 - x_diff) +  D*(x_diff*y_diff));
+
+            Scaled[i*w2+j] = gray;
+        }}
+    return *Scaled;
 }
diff --git a/src/Core/regularisers_CPU/utils.h b/src/Core/regularisers_CPU/utils.h
index cfaf6d7..8f0ba59 100644
--- a/src/Core/regularisers_CPU/utils.h
+++ b/src/Core/regularisers_CPU/utils.h
@@ -29,6 +29,8 @@ CCPI_EXPORT unsigned char copyIm_unchar(unsigned char *A, unsigned char *U, int
 CCPI_EXPORT float copyIm_roll(float *A, float *U, int dimX, int dimY, int roll_value, int switcher);
 CCPI_EXPORT float TV_energy2D(float *U, float *U0, float *E_val, float lambda, int type, int dimX, int dimY);
 CCPI_EXPORT float TV_energy3D(float *U, float *U0, float *E_val, float lambda, int type, int dimX, int dimY, int dimZ);
+CCPI_EXPORT float TV_energy3D(float *U, float *U0, float *E_val, float lambda, int type, int dimX, int dimY, int dimZ);
+CCPI_EXPORT float Im_scale2D(float *Input, float *Scaled, int w, int h, int w2, int h2);
 #ifdef __cplusplus
 }
 #endif
diff --git a/src/Matlab/mex_compile/supp_routines/Imscale2D.c b/src/Matlab/mex_compile/supp_routines/Imscale2D.c
new file mode 100644
index 0000000..4576cce
--- /dev/null
+++ b/src/Matlab/mex_compile/supp_routines/Imscale2D.c
@@ -0,0 +1,50 @@
+/*
+ * 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 2019 Daniil Kazantsev
+ * Copyright 2019 Srikanth Nagella, Edoardo Pasca
+ * Copyright 2019 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 "utils.h"
+
+/**************************************************/
+void mexFunction(
+        int nlhs, mxArray *plhs[],
+        int nrhs, const mxArray *prhs[])
+{
+    int number_of_dims,  X_new, Y_new;
+    mwSize dimX, dimY, dimZ;
+    const mwSize *dim_array;
+    float *A, *B;
+    mwSize dim_array2[2]; /* for scaled 2D 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 */
+    X_new = (int) mxGetScalar(prhs[1]);  /* new size for image */
+    Y_new = (int) mxGetScalar(prhs[2]);  /* new size for image */
+
+    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+    dim_array2[0] = X_new; dim_array2[1] = Y_new;
+    
+    B = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array2, mxSINGLE_CLASS, mxREAL));
+    
+    Im_scale2D(A, B, dimX, dimY, X_new, Y_new);
+ }