Merge pull request #34 from dkazanc/numproj

Adding rotation based forward projector
dkazanc · Oct 1, 2018 · d78422d · d78422d
2 parents 0212845 + f06aeeb
commit d78422d
Show file tree

Hide file tree

Showing 13 changed files with 315 additions and 31 deletions.
diff --git a/README.md b/README.md
@@ -32,7 +32,7 @@
  * [MATLAB](www.mathworks.com/products/matlab/) 
  * C compilers: GCC/MinGW/[TDM-GCC](http://tdm-gcc.tdragon.net/)/Visual Studio
 
-### Other (optional) dependencies (reconstruction):
+### Other (optional) dependencies:
  * [ASTRA-toolbox](http://www.astra-toolbox.com/)
  * [TomoPy](http://tomopy.readthedocs.io)
 

diff --git a/functions/TomoP2DSinoNum.c b/functions/TomoP2DSinoNum.c
@@ -0,0 +1,70 @@
+/*
+ * Copyright 2017 Daniil Kazantsev
+ *
+ * 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 "mex.h"
+#include <matrix.h>
+#include <math.h>
+#include <stdlib.h>
+#include <memory.h>
+#include <stdio.h>
+#include <malloc.h>
+#include "omp.h"
+
+#include "TomoP2DSinoNum_core.h"
+#include "utils.h"
+
+
+/* C OMP implementation of the forward projection (The Radon Transform)
+ * by rotating a padded image and summing over columns (rotation-based projector 
+ * for parallel beam)
+ *
+ * Input Parameters:
+ * 1. Phantom to calculate projections from [required]
+ * 2. Detector array size P (in pixels) [required]
+ * 3. Projection angles Theta (in degrees) [required]
+ *
+ * Output:
+ * Sinogram [No. angles] x [No. detectors]
+ *
+ * mex ForwardProjNum.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+ * sinogram = ForwardProjNum(single(G), P, single(angles));
+ */
+
+void mexFunction(
+        int nlhs, mxArray *plhs[],
+        int nrhs, const mxArray *prhs[]) {
+    int dimX, dimY,ThetaLength, DetSize, sys=1;
+    const int  *dim_array, *dimsTh;
+    float *Phantom, *Sinogram, *Theta;    
+
+    /*Handling Matlab input data*/
+    Phantom  = (float *) mxGetData(prhs[0]);
+    DetSize  = (int) mxGetScalar(prhs[1]);    
+    Theta  = (float *) mxGetData(prhs[2]);
+
+    dim_array = mxGetDimensions(prhs[0]);
+    dimsTh = mxGetDimensions(prhs[2]);
+
+    dimX = dim_array[0]; dimY = dim_array[1]; /*Image Dimensions (must be squared) */
+    if( dimX != dimY) mexErrMsgTxt("The input image must be squared");
+    if( DetSize < dimX) mexErrMsgTxt("The detector dimension is smaller than the phantom dimension! ");
+
+    ThetaLength = dimsTh[1]; /* angles dimension */           
+
+    /*Handling Matlab output data*/
+    const mwSize N_dims[2] = {ThetaLength, DetSize}; /*format: Y-angles x X-detectors dim*/
+    Sinogram = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, N_dims, mxSINGLE_CLASS, mxREAL));
+
+    TomoP2DSinoNum_core(Sinogram, Phantom, dimX, DetSize, Theta, ThetaLength, sys);
+}
diff --git a/functions/TomoP2DSinoNum_core.c b/functions/TomoP2DSinoNum_core.c
@@ -0,0 +1,133 @@
+/*
+ * Copyright 2017 Daniil Kazantsev
+ *
+ * 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 <math.h>
+#include <stdlib.h>
+#include <memory.h>
+#include <stdio.h>
+#include "omp.h"
+
+#include "utils.h"
+#include "TomoP2DSinoNum_core.h"
+
+#define M_PI 3.14159265358979323846
+#define EPS 0.000000001
+#define MAXCHAR 1000
+
+#define MAX(x, y) (((x) > (y)) ? (x) : (y))
+#define MIN(x, y) (((x) < (y)) ? (x) : (y))
+
+/* C OMP implementation of the forward projection (The Radon Transform)
+ * by rotating a padded image and summing over columns (rotation-based projector 
+ * for parallel beam)
+ *
+ * Input Parameters:
+ * 1. Phantom to calculate projections from [required]
+ * 2. Detector array size P (in pixels) [required]
+ * 3. Projection angles Theta (in degrees) [required]
+ *
+ * Output:
+ * Sinogram [No. angles] x [No. detectors]
+ */
+
+
+float TomoP2DSinoNum_core(float *Sinogram, float *Phantom, int dimX, int DetSize, float *Theta, int ThetaLength, int sys)
+{    
+	int i, j, k, padXY;
+	float *Phantom_pad=NULL, *B=NULL,  angC, ct, st, sumSin;
+
+    padXY = ceil(0.5f*(DetSize - dimX)); /*Size for padding the phantom*/ 
+
+    /*Perform padding of the phantom to the size [DetSize x DetSize] */
+    Phantom_pad = (float*) calloc(DetSize*DetSize,sizeof(float));  /*allocating space*/
+    padding(Phantom, Phantom_pad, DetSize, dimX, padXY, sys);
+
+    /* setting OMP here */
+    #pragma omp parallel for shared (Phantom_pad, Sinogram, dimX, DetSize, Theta, ThetaLength) private(B, k, j, i, sumSin, angC, ct, st)    
+    for (k=0; k < ThetaLength; k++) {
+
+        B = (float*) calloc(DetSize*DetSize,sizeof(float));
+
+        angC = Theta[k]*M_PI/180.0f;
+        ct = cos(angC + 0.5f*M_PI);
+        st = sin(angC + 0.5f*M_PI);
+
+        BilinearInterpolation(Phantom_pad, B, DetSize, ct, st); /* perform interpolation to rotate image on angle angC */
+
+        for (j=0; j < DetSize; j++) {
+            sumSin = 0.0f;
+            Sinogram[j*ThetaLength + k] = 0.0f;
+            for (i=0; i < DetSize; i++) sumSin += B[i*DetSize+j];
+            Sinogram[j*ThetaLength + k] = sumSin;
+        }
+        free(B);
+    }    
+
+    /*freeing the memory*/
+   free(Phantom_pad);
+   return  *Sinogram;     
+}    
+
+float BilinearInterpolation(float *Phantom_pad, float *B, int DetSize, float ct, float st)
+{
+    int i, j, k, i0, j0, i1, j1;
+    float *xs, H_x, s_min, s_max, stepS, x_rs, y_rs, dhalf, xbar, ybar;
+    dhalf = 0.5f*DetSize;
+
+    xs = (float*)calloc(DetSize,sizeof(float));
+    /*calculate grid*/
+    H_x = 1.0f/DetSize;
+    s_min = -1.0f + H_x;
+    s_max = 1.0f - H_x;
+    stepS = 2.0f/DetSize;
+    for (k=0; k < DetSize; k++) {xs[k] = s_min + (k)*(stepS);}
+
+    for (i=0; i < DetSize; i++) {
+        for (j=0; j < DetSize; j++) {
+
+            x_rs = dhalf*(xs[j]*ct - xs[i]*st + s_max);
+            y_rs = dhalf*(xs[j]*st + xs[i]*ct + s_max);
+
+            i1 = MIN(floor(y_rs), DetSize+2);
+            j1 = MIN(floor(x_rs), DetSize+2);
+            i0 = MAX(i1, 0);
+            j0 = MAX(j1, 0);
+
+            xbar = x_rs - j0;
+            ybar = y_rs - i0;
+
+            i0 = i0+1; j0=j0+1;
+
+            if ((i0 < DetSize) && (j0 < DetSize)) {
+                B[i*DetSize+j] = Phantom_pad[i0*DetSize+j0]*(1.-xbar)*(1.-ybar)+Phantom_pad[(i0+1)*DetSize+j0]*ybar*(1.-xbar)+Phantom_pad[i0*DetSize+(j0+1)]*xbar*(1.-ybar)+Phantom_pad[(i0+1)*DetSize+(j0+1)]*xbar*ybar;
+            }
+        }}
+    free(xs);
+    return *B;
+}
+
+float padding(float *Phantom, float *Phantom_pad, int DetSize, int PhantSize, int padXY, int sys)
+{
+    int i,j;   
+    #pragma omp parallel for shared (Phantom_pad, Phantom) private(i, j)
+    for (i=0; i < DetSize; i++) {
+        for (j=0; j < DetSize; j++) {
+            Phantom_pad[j*DetSize+i] = 0.0f;
+             if (((i >= padXY+1) && (i < DetSize-padXY)) &&  ((j >= padXY) && (j < DetSize-padXY))) {
+                  if (sys == 0) Phantom_pad[j*DetSize+i] = Phantom[(i-padXY)*PhantSize + (j-padXY)];
+                  else Phantom_pad[j*DetSize+i] = Phantom[(j-padXY)*PhantSize + (i-padXY)];
+             }
+        }}
+    return *Phantom_pad;
+}
diff --git a/functions/TomoP2DSinoNum_core.h b/functions/TomoP2DSinoNum_core.h
@@ -0,0 +1,23 @@
+/*
+Copyright 2017 Daniil Kazantsev
+
+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 <math.h>
+#include <stdlib.h>
+#include <memory.h>
+#include <stdio.h>
+#include "omp.h"
+
+float TomoP2DSinoNum_core(float *Sinogram, float *Phantom, int dimX, int DetSize, float *Theta, int ThetaLength, int sys);
+float BilinearInterpolation(float *Phantom_pad, float *B, int DetSize, float ct, float st);
+float padding(float *Phantom, float *Phantom_pad, int DetSize, int PhantSize, int padXY, int sys);
diff --git a/functions/utils.c b/functions/utils.c
@@ -294,7 +294,8 @@ float checkParams3D(int *params_switch, int ModelSelected, char *ModelParameters
             /* work with non-# commented lines */
             if(str[0] != '#') {
                 sscanf(str, "%15s : %21[^;];", tmpstr1, tmpstr2);
-                if (strcmp(tmpstr1,"Model")==0) Model = atoi(tmpstr2);   {
+                if (strcmp(tmpstr1,"Model")==0) {
+                    Model = atoi(tmpstr2);  
                     if ((ModelSelected == Model) && (counter == 0)) {
                         /* check if we have a right model */
                         if (fgets(str, MAXCHAR, fp) != NULL) sscanf(str, "%15s : %21[^;];", tmpstr1, tmpstr2);

diff --git a/matlab/Model2DGeneratorDemo.m b/matlab/Model2DGeneratorDemo.m
@@ -14,7 +14,7 @@
 addpath(pathtoModels);
 addpath('compiled'); addpath('supplem'); 
 
-ModelNo = 4; % Select a model from Phantom2DLibrary.dat
+ModelNo = 1; % Select a model from Phantom2DLibrary.dat
 % Define phantom dimensions
 N = 512; % x-y size (squared image)
 
@@ -28,13 +28,15 @@
 %%
 fprintf('%s \n', 'Generating sinogram analytically and numerically with Matlab (radon)...');
 % generate angles
-angles = linspace(0,180,N); % projection angles
+angles = linspace(0,179.9,round(1.3*N)); % projection angles
 
 % lets use Matlab's radon function
 [F_d,xp] = radon(G,angles); % discrete sinogram
 P = size(F_d,1); %detectors dimension
 F_d = F_d';
 
+% [F_num] = TomoP2DSinoNum(G, P, single(angles)); % numerical sinogram
+
 % generate the 2D analytical parallel beam sinogram
 [F_a] = TomoP2DModelSino(ModelNo, N, P, single(angles), pathTP, 'radon'); 
 

diff --git a/matlab/install/compile_mex_linux.m b/matlab/install/compile_mex_linux.m
@@ -24,6 +24,9 @@
 movefile('TomoP3DModel.mex*',Pathmove);
 mex TomoP3DObject.c TomoP3DModel_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
 movefile('TomoP3DObject.mex*',Pathmove);
+mex TomoP2DSinoNum.c TomoP2DSinoNum_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+movefile('TomoP2DSinoNum.mex*',Pathmove);
+
 fprintf('%s \n', 'All compiled!');
 
 cd(UpPath);

diff --git a/matlab/install/compile_mex_windows.m b/matlab/install/compile_mex_windows.m
@@ -32,6 +32,8 @@
 movefile('TomoP3DModel.mex*',Pathmove);
 mex TomoP3DObject.c TomoP3DModel_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99"
 movefile('TomoP3DObject.mex*',Pathmove);
+mex TomoP2DSinoNum.c TomoP2DSinoNum_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99"
+movefile('TomoP2DSinoNum.mex*',Pathmove);
 
 %%% The second approach to compile using TDM-GCC which follows this
 %%% discussion:
@@ -55,7 +57,9 @@
 % movefile('TomoP3DModel.mex*',Pathmove);
 % mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" TomoP3DObject.c TomoP3DModel_core.c utils.c
 % movefile('TomoP3DObject.mex*',Pathmove);
+% mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" TomoP2DSinoNum.c TomoP2DSinoNum_core.c utils.c
+% movefile('TomoP2DSinoNum.mex*',Pathmove);
 % fprintf('%s \n', 'All compiled!');
 
 cd(UpPath);
-cd matlab
+cd matlab