Commits

kritisen committed 8c6331d

renamed openct to openrecon, fanrec to xrec, and fanprj to xprj. Also removed some unused codes. Also moved sample config files out.

Comments (0)

Files changed (24)

 *.cfg
 utilities/_weightForShortScanArchives/*
 log
+*.mex*
+delete/*
+sampleConfig/*
-% for fanprj
-mex fanprj.cpp prj.cpp prj_fanED.cpp prj_fanEA.cpp prj_para.cpp img.cpp img_fltTV.cpp img_fltSTH.cpp cpputils.cpp
+% for xprj
+mex xprj.cpp prj.cpp prj_fanED.cpp prj_fanEA.cpp prj_para.cpp img.cpp img_fltTV.cpp img_fltSTH.cpp cpputils.cpp
 
-% for fanrec
-mex fanrec.cpp prj.cpp prj_fanED.cpp prj_fanEA.cpp prj_para.cpp img.cpp img_fltTV.cpp img_fltSTH.cpp rec.cpp rec_null.cpp rec_backPrj.cpp rec_SART.cpp rec_SIRT.cpp rec_FBP.cpp cpputils.cpp
+% for xrec
+mex xrec.cpp prj.cpp prj_fanED.cpp prj_fanEA.cpp prj_para.cpp img.cpp img_fltTV.cpp img_fltSTH.cpp rec.cpp rec_null.cpp rec_backPrj.cpp rec_SART.cpp rec_SIRT.cpp rec_FBP.cpp cpputils.cpp
 
 % build the utilities
 mex utilities/imgTVGradMin.cpp

build_unix.m

-% for fanprj
-mex CXXFLAGS="\$CXXFLAGS -fopenmp" LDFLAGS="\$LDFLAGS -fopenmp -lrt" -lgomp fanprj.cpp prj.cpp prj_fanED.cpp prj_fanEA.cpp prj_para.cpp img.cpp img_fltTV.cpp img_fltSTH.cpp -lgomp
-
-% for fanrec
-mex CXXFLAGS="\$CXXFLAGS -fopenmp" LDFLAGS="\$LDFLAGS -fopenmp -lrt" -lgomp fanrec.cpp prj.cpp prj_fanED.cpp prj_fanEA.cpp prj_para.cpp img.cpp img_fltTV.cpp img_fltSTH.cpp rec.cpp rec_null.cpp rec_backPrj.cpp rec_SART.cpp rec_SIRT.cpp
 I(idx)=0;
 figure;imshow(I, [1.0 1.05]);title('original phantom image');
 
-% P_DBP = fanprj(I,'config_dbp.cfg');
-% I_DBP = fanrec(P_DBP,'config_dbp.cfg');
+% P_DBP = xprj(I,'config_dbp.cfg');
+% I_DBP = xrec(P_DBP,'config_dbp.cfg');
 % figure;imshow(I_DBP, []);
 % title('direct back-projection');
 
-P_PARA=fanprj(I,'config_para.cfg');
-I_PARA=fanrec(P_PARA,'config_para.cfg');
+P_PARA=xprj(I,'config_para.cfg');
+I_PARA=xrec(P_PARA,'config_para.cfg');
 figure;imshow(I_PARA, [1.0 1.05]);
 title('para-beam reconstuction');
 
-P_fan=fanprj(I,'config_fan.cfg');
-I_fan=fanrec(P_fan,'config_fan.cfg');
+P_fan=xprj(I,'config_fan.cfg');
+I_fan=xrec(P_fan,'config_fan.cfg');
 figure;imshow(I_fan, [1.0 1.05]);
 title('fan-beam reconstuction');
 
 % use default.cfg
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-% P_DEF=fanprj(I);
-% I_DEF=fanrec(P_DEF);
+% P_DEF=xprj(I);
+% I_DEF=xrec(P_DEF);
 % figure;imshow(I_DEF, [1.0 1.05]);title('reconstucted image using default config');
 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 % % ROI problem
-% P_ROI=fanprj(I,'config_SART_TV_ROI.cfg');
-% I_ROI=fanrec(P_ROI,'config_SART_TV_ROI.cfg');
+% P_ROI=xprj(I,'config_SART_TV_ROI.cfg');
+% I_ROI=xrec(P_ROI,'config_SART_TV_ROI.cfg');
 % figure;imshow(I_ROI, [1.0 1.05]);title('reconstucted image for ROI problem');
 
 % % few views problem
-% P_FV=fanprj(I,'config_fewviews_SART.cfg');
-% I_FV=fanrec(P_FV,'config_fewviews_SART.cfg');
+% P_FV=xprj(I,'config_fewviews_SART.cfg');
+% I_FV=xrec(P_FV,'config_fewviews_SART.cfg');
 % figure;imshow(I_FV, [1.0 1.05]);title('reconstucted image for few views problem');
 
 % % Reconstruction with Soft Thresholding Regularization
-% P_STH=fanprj(I,'config_SART_STH_YU2009.cfg');
-% I_STH=fanrec(P_STH,'config_SART_STH_YU2009.cfg');
+% P_STH=xprj(I,'config_SART_STH_YU2009.cfg');
+% I_STH=xrec(P_STH,'config_SART_STH_YU2009.cfg');
 % figure;imshow(I_STH, [1.0 1.05]);title('Reconstruction with Soft Thresholding Regularization');
 
 % % Reconstruction with TV Regularization
-% P_TV = fanprj(I,'config_SART_TV_Sidky2006.cfg');
-% I_TV = fanrec(P_TV,'config_SART_TV_Sidky2006.cfg');
+% P_TV = xprj(I,'config_SART_TV_Sidky2006.cfg');
+% I_TV = xrec(P_TV,'config_SART_TV_Sidky2006.cfg');
 % figure;imshow(I_TV, [1.0 1.05]);title('Reconstruction with TV Regularization');
 
 

fanIO.cpp

-//////////////////////////////////
-//  CT/Micro CT lab
-//  Biomedical Imaging Divison
-//  SBES, Virginia Tech.
-//  Version of 2007.07.24
-//////////////////////////////////
-
-#include <mex.h>
-#include <string.h>
-#include <stdio.h>
-#include <math.h>
-
-int showUsage()
-{
-    mexPrintf("fanIO - Fanbeam Prjection data import and outport function v1.0\n");
-    mexPrintf("\tP = fanIO(I) calculates projection P from image I using default configuration (read from file 'default.cfg').\n");
-    
-    return 0;
-}
-
-void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
-{
-    //no any inputs, print usage
-    if(nrhs==0){
-        showUsage();
-        return;
-    }
-
-    //1 input P = fanprj(I)
-    if(nrhs==1){
-        if( mxGetClassID(prhs[0]) == mxDOUBLE_CLASS ){
-            mexPrintf("Loading default configuration file 'default.cfg'.\n");
-            
-            hImg = Img_create(def, mxGetPr(prhs[0]));
-            if(hImg==NULL){
-                mexPrintf("Error: Load [IMAGE] config failed. Please check the .cfg file.\n");
-                goto EE;
-            }
-            
-            hPrj = Prj_create(def, 0);
-            if(hPrj==NULL){
-                mexPrintf("Error: Load [SCANNING_GEOMETRY] config failed. Please check the .cfg file.\n\n");
-                goto EE;
-            }
-            
-			//calc prj
-            Prj_calcPrj(hPrj, hImg);
-
-            plhs[0] = mxCreateDoubleMatrix(hPrj->pa.detN,hPrj->pa.angN,mxREAL);
-			memcpy(mxGetPr(plhs[0]), hPrj->prj, hPrj->pa.detN*hPrj->pa.angN*sizeof(double));
-            
-            Prj_delete(hPrj);
-            Img_delete(hImg);
-
-            return;
-        }
-        else{
-            mexPrintf("Error: Invalid input argument.\n");
-            goto EE;
-        }
-    }
-    
-    //2 inputs P = fanprj(I,'CFG_FILENAME')
-    if(nrhs==2){
-        if(mxGetClassID(prhs[1]) == mxCHAR_CLASS && mxGetClassID(prhs[0]) == mxDOUBLE_CLASS ){
-            char * c = mxArrayToString(prhs[1]);
-            if(!strstr(c,".cfg")){
-				mexPrintf("Error: Invalid config file.\n");
-				goto EE;
-			}
-
-			hImg = Img_create(c, mxGetPr(prhs[0]));
-			if(hImg==NULL){
-				mexPrintf("Error: Load [IMAGE] config failed. Please check the .cfg file.\n");
-				goto EE;
-			}
-
-			hPrj = Prj_create(c, 0);
-			if(hPrj==NULL){
-				mexPrintf("Error: Load [SCANNING_GEOMETRY] config failed. Please check the .cfg file.\n\n");
-				goto EE;
-			}
- 
-			//calc prj
-            Prj_calcPrj(hPrj, hImg);
-
-            plhs[0] = mxCreateDoubleMatrix(hPrj->pa.detN,hPrj->pa.angN,mxREAL);
-			memcpy(mxGetPr(plhs[0]), hPrj->prj, hPrj->pa.detN*hPrj->pa.angN*sizeof(double));
-            
-            Prj_delete(hPrj);
-            Img_delete(hImg);
-
-            return;
-
-        }
-		else{
-			mexPrintf("Error: Invalid input arguments.\n");
-            goto EE;
-		}
-                
-    }
-
-
-EE:
-    mexPrintf("Request can not be handled.\n\n");
-    showUsage();    
-
-    return;
-}

fanprj.cpp

-//////////////////////////////////
-//  CT/Micro CT lab
-//  Biomedical Imaging Divison
-//  SBES, Virginia Tech.
-//  Version of 2007.07.24
-//////////////////////////////////
-
-#include <mex.h>
-#include <string.h>
-#include <stdio.h>
-#include <math.h>
-
-#include "img.h"
-#include "prj.h"
-
-int showUsage()
-{
-    mexPrintf("fanprj - Fanbeam Prjection Kernel v1.0\n");
-    mexPrintf("\tP = fanprj(I) calculates projection P from image I using default configuration (read from file 'default.cfg').\n");
-    mexPrintf("\tP = fanprj(I,'CFG_FILENAME') calculates projection P from image I using configuration from file 'CFG_FILENAME' (must be a .cfg file).\n");
-    
-    return 0;
-}
-
-void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
-{
-    Img_Handle hImg=NULL;
-    Prj_Handle hPrj=NULL;
-    
-    char defcfg[] = "default.cfg";
-    char *def = defcfg;
-    
-    //no any inputs, print usage
-    if(nrhs==0){
-        showUsage();
-        return;
-    }
-
-    //1 input P = fanprj(I)
-    if(nrhs==1){
-        if( mxGetClassID(prhs[0]) == mxDOUBLE_CLASS ){
-            mexPrintf("Loading default configuration file 'default.cfg'.\n");
-            
-            hImg = Img_create(def, mxGetPr(prhs[0]));
-            if(hImg==NULL){
-                mexPrintf("Error: Load [IMAGE] config failed. Please check the .cfg file.\n");
-                goto EE;
-            }
-            
-            hPrj = Prj_create(def, 0);
-            if(hPrj==NULL){
-                mexPrintf("Error: Load [SCANNING_GEOMETRY] config failed. Please check the .cfg file.\n\n");
-                goto EE;
-            }
-            
-			//calc prj
-            Prj_calcPrj(hPrj, hImg);
-
-            plhs[0] = mxCreateDoubleMatrix(hPrj->pa.detN,hPrj->pa.angN,mxREAL);
-			memcpy(mxGetPr(plhs[0]), hPrj->prj, hPrj->pa.detN*hPrj->pa.angN*sizeof(double));
-            
-            Prj_delete(hPrj);
-            Img_delete(hImg);
-
-            return;
-        }
-        else{
-            mexPrintf("Error: Invalid input argument.\n");
-            goto EE;
-        }
-    }
-    
-    //2 inputs P = fanprj(I,'CFG_FILENAME')
-    if(nrhs==2){
-        if(mxGetClassID(prhs[1]) == mxCHAR_CLASS && mxGetClassID(prhs[0]) == mxDOUBLE_CLASS ){
-            char * c = mxArrayToString(prhs[1]);
-            if(!strstr(c,".cfg")){
-				mexPrintf("Error: Invalid config file.\n");
-				goto EE;
-			}
-
-			hImg = Img_create(c, mxGetPr(prhs[0]));
-			if(hImg==NULL){
-				mexPrintf("Error: Load [IMAGE] config failed. Please check the .cfg file.\n");
-				goto EE;
-			}
-
-			hPrj = Prj_create(c, 0);
-			if(hPrj==NULL){
-				mexPrintf("Error: Load [SCANNING_GEOMETRY] config failed. Please check the .cfg file.\n\n");
-				goto EE;
-			}
- 
-			//calc prj
-            Prj_calcPrj(hPrj, hImg);
-
-            plhs[0] = mxCreateDoubleMatrix(hPrj->pa.detN,hPrj->pa.angN,mxREAL);
-			memcpy(mxGetPr(plhs[0]), hPrj->prj, hPrj->pa.detN*hPrj->pa.angN*sizeof(double));
-            
-            Prj_delete(hPrj);
-            Img_delete(hImg);
-
-            return;
-
-        }
-		else{
-			mexPrintf("Error: Invalid input arguments.\n");
-            goto EE;
-		}
-                
-    }
-
-
-EE:
-    mexPrintf("Request can not be handled.\n\n");
-    showUsage();    
-
-    if(hPrj)
-		Prj_delete(hPrj);
-	if(hImg)
-		Img_delete(hImg);
-    
-    return;
-}

fanrec.cpp

-//////////////////////////////////
-//  CT/Micro CT lab
-//  Biomedical Imaging Divison
-//  SBES, Virginia Tech.
-//  Version of 2007.07.24
-//////////////////////////////////
-
-#include <mex.h>
-#include <string.h>
-#include <stdio.h>
-#include <math.h>
-
-#include "img.h"
-#include "prj.h"
-#include "rec.h"
-
-int showUsage()
-{
-    mexPrintf("fanrec - Fanbeam Reconstruction Kernel v1.0\n");
-    mexPrintf("\tI = fanrec(P) reconstructs image I from projections P using default configuration (read from file 'default.cfg').\n");
-    mexPrintf("\tI = fanrec(P, 'CFG_FILENAME') reconstructs image I using configuration file 'CFG_FILENAME' (must be a .cfg file).\n");
-    mexPrintf("\tI = fanrec(P, I0) same as I=fanrec(P) but with a initial image value I0 given.\n");
-    mexPrintf("\tI = fanrec(P, I0, 'CFG_FILENAME') given both initial image value I0 and configuration file 'CFG_FILENAME' (must be a .cfg file).\n");
-    
-    return 0;
-}
-
-void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
-{
-    Img_Handle hImg=NULL;
-    Prj_Handle hPrj=NULL;
-    Rec_Handle hRec=NULL;
-    
-    char defcfg[] = "default.cfg";
-    char * def = defcfg;
-    int prm_set = 0;
-    
-    //no any inputs, print usage
-    if(nrhs==0){
-        showUsage();
-        return;
-    }
-
-    //1 input
-    if(nrhs==1){
-        if( mxGetClassID(prhs[0]) == mxDOUBLE_CLASS ){
-            mexPrintf("Loading default configuration file 'default.cfg'.\n");
-            
-            hImg = Img_create(def,0);
-            if(hImg==NULL){
-                mexPrintf("Error: Load [IMAGE] config failed. Please check the default.cfg file.\n");
-                goto EE;
-            }
-            
-            hPrj = Prj_create(def, mxGetPr(prhs[0]));
-            if(hPrj==NULL){
-                mexPrintf("Error: Load [SCANNING_GEOMETRY] config failed. Please check the default.cfg file.\n\n");
-                goto EE;
-            }
-            
-            hRec = Rec_create(def);
-            if(hRec==NULL){
-                mexPrintf("Error: Load [RECONSTRUCTION] config failed. Please check the default.cfg file.\n");
-                goto EE;
-            }
-             
-            Rec_recImg(hRec, hPrj, hImg);
-
-			plhs[0] = mxCreateDoubleMatrix(hImg->ia.rowN, hImg->ia.colN, mxREAL);
-			memcpy(mxGetPr(plhs[0]), hImg->img, hImg->ia.rowN*hImg->ia.colN*sizeof(double));
-
-            Rec_delete(hRec);
-            Prj_delete(hPrj);
-            Img_delete(hImg);
-
-            return;
-        }
-        else{
-            mexPrintf("Invalid data format. Please make sure the input image data is in double format.\n");
-            goto EE;
-        }
-    }
-    
-    //2 inputs
-    if(nrhs==2){
-        
-        //I = fanrec(P, 'CFG_FILENAME')
-        if(mxGetClassID(prhs[1]) == mxCHAR_CLASS && mxGetClassID(prhs[0]) == mxDOUBLE_CLASS){
-
-			//load config from file
-            char * c = mxArrayToString(prhs[1]);
-            if(!strstr(c,".cfg")){
-				mexPrintf("Error: invalid configuration file %s.\n", c);
-				goto EE;
-			}
-
-			mexPrintf("Load config from file %s.\n",c);
-
-			hImg = Img_create(c,0);
-			if(hImg==NULL){
-				mexPrintf("Error: Load [IMAGE] config failed. Please check the .cfg file.\n");
-				goto EE;
-			}
-
-			hPrj = Prj_create(c,mxGetPr(prhs[0]));
-			if(hPrj==NULL){
-				mexPrintf("Error: Load [SCANNING_GEOMETRY] config failed. Please check the .cfg file.\n\n");
-				goto EE;
-			}
-
-			hRec = Rec_create(c);
-			if(hRec==NULL){
-				mexPrintf("Error: Load [RECONSTRUCTION] config failed. Please check the .cfg file.\n");
-				goto EE;
-			}
-
-			//fanrec
-            Rec_recImg(hRec, hPrj, hImg);
-
-			plhs[0] = mxCreateDoubleMatrix(hImg->ia.rowN, hImg->ia.colN, mxREAL);
-			memcpy(mxGetPr(plhs[0]), hImg->img, hImg->ia.rowN*hImg->ia.colN*sizeof(double));
-
-			Rec_delete(hRec);
-			Prj_delete(hPrj);
-			Img_delete(hImg);
-
-            return;
-            
-        }
-		//I = fanrec(P, I0)
-        else if(mxGetClassID(prhs[1]) == mxDOUBLE_CLASS && mxGetClassID(prhs[0]) == mxDOUBLE_CLASS){
-
-			mexPrintf("Load default config.\n",def);
-
-			hImg = Img_create(def, mxGetPr(prhs[1]));
-			if(hImg==NULL){
-				mexPrintf("Error: Load [IMAGE] config failed. Please check the .cfg file.\n");
-				goto EE;
-			}
-
-			hPrj = Prj_create(def, mxGetPr(prhs[0]));
-			if(hPrj==NULL){
-				mexPrintf("Error: Load [SCANNING_GEOMETRY] config failed. Please check the .cfg file.\n\n");
-				goto EE;
-			}
-
-			hRec = Rec_create(def);
-			if(hRec==NULL){
-				mexPrintf("Error: Load [RECONSTRUCTION] config failed. Please check the .cfg file.\n");
-				goto EE;
-			}
-
-			//fanrec
-            Rec_recImg(hRec, hPrj, hImg);
-
-			plhs[0] = mxCreateDoubleMatrix(hImg->ia.rowN, hImg->ia.colN, mxREAL);
-			memcpy(mxGetPr(plhs[0]), hImg->img, hImg->ia.rowN*hImg->ia.colN*sizeof(double));
-
-			Rec_delete(hRec);
-			Prj_delete(hPrj);
-			Img_delete(hImg);
-
-            return;
-
-		}
-		else{
-            mexPrintf("Error: Invalid argument.\n");
-            goto EE;
-        }
-    }
-    
-    //3 inputs
-    if(nrhs==3){ 
-
-        //I = fanrec(P, I0, 'CFG_FILENAME')
-        if(mxGetClassID(prhs[0]) == mxDOUBLE_CLASS && mxGetClassID(prhs[1]) == mxDOUBLE_CLASS && mxGetClassID(prhs[2]) == mxCHAR_CLASS){
-
-            char * c = mxArrayToString(prhs[2]);
-            if(!strstr(c,".cfg")){
-				mexPrintf("Error: invalid configuration file %s.\n", c);
-				goto EE;
-			}
-
-			mexPrintf("Load config from file %s.\n",c);mexEvalString("drawnow;");
-
-			hImg = Img_create(c, mxGetPr(prhs[1]));
-			if(hImg==NULL){
-				mexPrintf("Error: Load [IMAGE] config failed. Please check the .cfg file.\n");
-				goto EE;
-			}
-
-			hPrj = Prj_create(c, mxGetPr(prhs[0]));
-			if(hPrj==NULL){
-				mexPrintf("Error: Load [SCANNING_GEOMETRY] config failed. Please check the .cfg file.\n\n");
-				goto EE;
-			}
-
-			hRec = Rec_create(c);
-			if(hRec==NULL){
-				mexPrintf("Error: Load [RECONSTRUCTION] config failed. Please check the .cfg file.\n");
-				goto EE;
-			}
-
-			//fanrec
-            Rec_recImg(hRec, hPrj, hImg);
-
-			plhs[0] = mxCreateDoubleMatrix(hImg->ia.rowN, hImg->ia.colN, mxREAL);
-			memcpy(mxGetPr(plhs[0]), hImg->img, hImg->ia.rowN*hImg->ia.colN*sizeof(double));
-
-			Rec_delete(hRec);
-			Prj_delete(hPrj);
-			Img_delete(hImg);
-
-            return;
-		}
-		else{
-			mexPrintf("Error: invalid argument.\n");
-			goto EE;
-		}
-
-	}
-
-    
-EE:
-    mexPrintf("Request can not be handled.\n\n");
-    showUsage();    
-
-    if(hImg) Img_delete(hImg);
-    if(hPrj) Prj_delete(hPrj);
-    if(hRec) Rec_delete(hRec);
-    
-    return;    
-}
 				//	else if (ang >= L3 && ang <= L4)
 				//		factor = Ec;
 				//	prj[pa.detN*angi+deti] *= factor;
-				//} // **** This is carried out in MATLAB function >>weightForShortScan()<< cefore calling fanrec() 
+				//} // **** This is carried out in MATLAB function >>weightForShortScan()<< cefore calling xrec() 
 				//			This accounts for short scans which do not cover the whole angular range 2pi
 				//				but run through pi+2gammaM	***
 				
 				//	double s = (deti - pa.detN/2);     
 				//	prj[pa.detN*angi+deti] *= (D-(tau*s)/D)/(sqrt(D*D+s*s));
 				//} // **** This is carried out in MATLAB function >>weightFanProjection()<< 
-				//				before calling fanrec() 
+				//				before calling xrec() 
 				//			This function accounts for different FBP formula for fanbeam than parallel beam
 				//				and also accounts for center shift***
 				 
 				//filterProjection(prj + pa.detN*angi, pa.detN);	//filter the projection with the ramp 
 																//filter h(s) / smoothed ramp filter
-				// **** filterProjection() is carried out in MATLAB cefore calling fanrec(),
+				// **** filterProjection() is carried out in MATLAB cefore calling xrec(),
 				// by calling MATLAB code for filtering projections ***
 				
 				//for (int j = 0; j < imgN; j++) //loop through all image pixels

release/1.cfg

-// This is an modified config file. 
-// Usage:
-// For projection simulation:
-//   P = fanprj(I,'1.cfg');	     //where I is the 2D phantom image and P is the projection data
-// For reconstruction:
-//   I_new = fanrec(P,'1.cfg');  //where P is the projection data file and I_new is the reconstructed image
-
-[SCANNING_GEOMETRY]
-Scanning_Type = FAN_ED              //FAN_ED: Fan-beam with detector units arranged in line; FAN_EA: Fan-beam with units on an circular arc whose origin is the X-ray source; 
-XRay_To_Rotation_Center = 400
-Detector_To_Rotation_Center = 400
-Total_Num_Of_Angles = 22 
-First_Angle_Degree = 0 
-Angle_Increment = 16.36
-Length_Of_Detector = 512            //To FAN_ED type it means the total length of detectors, while to FAN_EA type it means total angular range.
-Num_Of_Detector_Units = 256
-Bias_Of_Detector_Center = 0
-Gaussian_Noise_Percentage = 0.0     //Adds Gaussian noise to detectors
-Store_Projection_Matrix = 1         //For iterative reconstruction ONLY. Will use a lot of memory to store projection matrix, but greatly save iteration time, set 0 to disable.
-
-[IMAGE]
-Rows_Of_Image = 200
-Cloumns_Of_Image = 200
-Pixel_Width = 1
-Pixel_Height = 1
-Bias_Of_Row_Center_To_Rotation_Center = 0
-Bias_Of_Column_Center_To_Rotation_Center = 0
-
-[RECONSTRUCTION]
-Reconstruction_Algorithm = SART	    //SART
-Num_Iterations = 20
-Lambda = 1
-Regularization = 1					//whether or not to use image regularization method in iteration. If to use, the [REGUALRAIZATION] section must be provided in this file.
-
-[REGULARIZATION]
-Regularization_Algorithm = TV		//TV(Total Variation) Minimization
-TV_Iterations = 20                  //iterations of each TV minimization loop
-TV_Alpha = 0.2                      //for TV ONLY
-

release/README.m

-% #1. A simple demonstration
-% The following lines creates a phantom image, 
-% and then performs projection and reconstruction.
-% The fanprj and fanrec will use default parameters.
-I=phantom;    % phantom is the Matlab default phantom creation function.
-P=fanprj(I);
-Irec=fanrec(P);
-figure;imshow(I2,[0 1.2]);
-
-
-% #2. Another demonstration showing how to change parameters
-% The following lines creates a phantom image, and then performs 
-% projection and reconstruction using specific parameters refering 
-% to file '1.cfg'. More details can be found in .cfg file. The .cfg file
-% can be open and modified directly using Matlab editor.
-I=phantom(200);
-P=fanprj(I,'1.cfg');  % image size has been changed to fit the phantom
-Irec=fanrec(P,'1.cfg');
-figure;imshow(I2,[0 1.2]);
-
-% #3. One can use template.cfg file to create a new config file for a
-% specific problem and use it when call fanprj and fanrec. There are
-% some config files provided for some well known problem and algorithm.
-% They can be referred while creating a new config file.
-

release/config_SART_STH_YU2009.cfg

-// This config file is an example to solve fewview problem using SART reconstruction method with TV regularization.
-// Usage:
-// For projection simulation:
-//   P = fanprj(I, 'SART_STH_YU2009.cfg');	    //where I is the 2D phantom image and P is the projection data
-// For reconstruction:
-//   I_new = fanrec(P, 'SART_STH_YU2009.cfg');  //where P is the projection data file and I_new is the reconstructed image
-
-
-[SCANNING_GEOMETRY]
-Scanning_Type = FAN_ED              //FAN_ED: Fan-beam with detector units arranged in line; FAN_EA: Fan-beam with units on an circular arc whose origin is the X-ray source; 
-XRay_To_Rotation_Center = 285
-Detector_To_Rotation_Center = 0
-Total_Num_Of_Angles = 21 
-First_Angle_Degree = 0 
-Angle_Increment = 17.1429
-Length_Of_Detector = 200            //To FAN_ED type it means the total length of detectors, while to FAN_EA type it means total angular range.
-Num_Of_Detector_Units = 300
-Bias_Of_Detector_Center = 0
-Gaussian_Noise_Percentage = 0.0     //Adds Gaussian noise to detectors
-Store_Projection_Matrix = 1         //For iterative reconstruction ONLY. Will use a lot of memory to store projection matrix, but greatly save iteration time, set 0 to disable.
-
-[IMAGE]
-Rows_Of_Image = 256
-Cloumns_Of_Image = 256
-Pixel_Width = 0.78125
-Pixel_Height = 0.78125
-Bias_Of_Row_Center_To_Rotation_Center = 0
-Bias_Of_Column_Center_To_Rotation_Center = 0
-
-[RECONSTRUCTION]
-Reconstruction_Algorithm = SART	    //SART
-Num_Iterations = 80
-Lambda = 1
-Regularization = 1					//whether or not to use image regularization method in iteration. If to use, the [REGUALRAIZATION] section must be provided in this file.
-
-[REGULARIZATION]
-Regularization_Algorithm = STH		//STH(Soft Threshold) Minimization
-STH_Image_Gradient_Type = DDT       //DGT: use DGT(Discrete Gradient Transform) to compute image gradient; DDT: use DDT(Discrete Difference Transform) to compute image gradient
-STH_Omega = 0.01                    //for STH method ONLY, the threshold Omega can be set according to the noise level
-

release/config_SART_TV_ROI.cfg

-// This config file is to solve ROI problem using SART reconstruction method with TV regularization.
-// Usage:
-// For projection simulation:
-//   P = fanprj(I, 'config_SART_TV_ROI.cfg');	    //where I is the 2D phantom image and P is the projection data
-// For reconstruction:
-//   I_new = fanrec(P, 'config_SART_TV_ROI.cfg');  //where P is the projection data file and I_new is the reconstructed image
-
-[SCANNING_GEOMETRY]
-Scanning_Type = FAN_ED              //FAN_ED: Fan-beam with detector units arranged in line; FAN_EA: Fan-beam with units on an circular arc whose origin is the X-ray source; 
-XRay_To_Rotation_Center = 400
-Detector_To_Rotation_Center = 400
-Total_Num_Of_Angles = 90 
-First_Angle_Degree = 0 
-Angle_Increment = 2 
-Length_Of_Detector = 256            //To FAN_ED type it means the total length of detectors, while to FAN_EA type it means total angular range.
-Num_Of_Detector_Units = 256
-Bias_Of_Detector_Center = 0
-Gaussian_Noise_Percentage = 0.0     //Adds Gaussian noise to detectors
-Store_Projection_Matrix = 1         //For iterative reconstruction ONLY. Will use a lot of memory to store projection matrix, but greatly save iteration time, set 0 to disable.
-
-[IMAGE]
-Rows_Of_Image = 256
-Cloumns_Of_Image = 256
-Pixel_Width = 1
-Pixel_Height = 1
-Bias_Of_Row_Center_To_Rotation_Center = 0
-Bias_Of_Column_Center_To_Rotation_Center = 0
-
-[RECONSTRUCTION]
-Reconstruction_Algorithm = SART	    //SART
-Num_Iterations = 40
-Lambda = 1
-Regularization = 1					//whether or not to use image regularization method in iteration. If to use, the [REGUALRAIZATION] section must be provided in this file.
-
-[REGULARIZATION]
-Regularization_Algorithm = TV		//TV(Total Variation) Minimization
-TV_Iterations = 20                  //iterations of each TV minimization loop
-TV_Alpha = 0.2                      //for TV ONLY
-

release/config_SART_TV_Sidky2006.cfg

-// This config file is to solve fewview problem using SART reconstruction method with TV regularization.
-// Usage:
-// For projection simulation:
-//   P = fanprj(I, 'config_SART_TV_Sidky2006.cfg');	    //where I is the 2D phantom image and P is the projection data
-// For reconstruction:
-//   I_new = fanrec(P, 'config_SART_TV_Sidky2006.cfg');    //where P is the projection data file and I_new is the reconstructed image
-
-[SCANNING_GEOMETRY]
-Scanning_Type = FAN_ED              //FAN_ED: Fan-beam with detector units arranged in line; FAN_EA: Fan-beam with units on an circular arc whose origin is the X-ray source; 
-XRay_To_Rotation_Center = 400
-Detector_To_Rotation_Center = 400
-Total_Num_Of_Angles = 22 
-First_Angle_Degree = 0 
-Angle_Increment = 16.36
-Length_Of_Detector = 512            //To FAN_ED type it means the total length of detectors, while to FAN_EA type it means total angular range.
-Num_Of_Detector_Units = 256
-Bias_Of_Detector_Center = 0
-Gaussian_Noise_Percentage = 0.0     //Adds Gaussian noise to detectors
-Store_Projection_Matrix = 1         //For iterative reconstruction ONLY. Will use a lot of memory to store projection matrix, but greatly save iteration time, set 0 to disable.
-
-[IMAGE]
-Rows_Of_Image = 256
-Cloumns_Of_Image = 256
-Pixel_Width = 1
-Pixel_Height = 1
-Bias_Of_Row_Center_To_Rotation_Center = 0
-Bias_Of_Column_Center_To_Rotation_Center = 0
-
-[RECONSTRUCTION]
-Reconstruction_Algorithm = SART	    //SART
-Num_Iterations = 20
-Lambda = 1
-Regularization = 1					//whether or not to use image regularization method in iteration. If to use, the [REGUALRAIZATION] section must be provided in this file.
-
-[REGULARIZATION]
-Regularization_Algorithm = TV		//TV(Total Variation) Minimization
-TV_Iterations = 20                  //iterations of each TV minimization loop
-TV_Alpha = 0.2                      //for TV ONLY
-

release/config_fewviews_SART.cfg

-// This config file is an example to solve fewview problem using SART reconstruction method.
-// Usage:
-// For projection simulation:
-//   P = fanprj(I, 'config_fewviews_SART.cfg');	     //where I is the 2D phantom image and P is the projection data
-// For reconstruction:
-//   I_new = fanrec(P, 'condig_fewviews_SART.cfg');  //where P is the projection data file and I_new is the reconstructed image
-
-[SCANNING_GEOMETRY]
-Scanning_Type = FAN_ED              //FAN_ED: Fan-beam with detector units arranged in line; FAN_EA: Fan-beam with units on an circular arc whose origin is the X-ray source; 
-XRay_To_Rotation_Center = 400
-Detector_To_Rotation_Center = 400
-Total_Num_Of_Angles = 22 
-First_Angle_Degree = 0 
-Angle_Increment = 16.36
-Length_Of_Detector = 512            //To FAN_ED type it means the total length of detectors, while to FAN_EA type it means total angular range.
-Num_Of_Detector_Units = 256
-Bias_Of_Detector_Center = 0
-Gaussian_Noise_Percentage = 0.0     //Adds Gaussian noise to detectors
-Store_Projection_Matrix = 1         //For iterative reconstruction ONLY. Will use a lot of memory to store projection matrix, but greatly save iteration time, set 0 to disable.
-
-[IMAGE]
-Rows_Of_Image = 256
-Cloumns_Of_Image = 256
-Pixel_Width = 1
-Pixel_Height = 1
-Bias_Of_Row_Center_To_Rotation_Center = 0
-Bias_Of_Column_Center_To_Rotation_Center = 0
-
-[RECONSTRUCTION]
-Reconstruction_Algorithm = SART	    //SART
-Num_Iterations = 20
-Lambda = 1
-

release/default.cfg

-// This is an default config file. 
-// *********************
-// ** NO MODIFICATION **
-// *********************
-// Usage:
-// For projection simulation:
-//   P = fanprj(I);	     //where I is the 2D phantom image and P is the projection data
-// For reconstruction:
-//   I_new = fanrec(P);  //where P is the projection data file and I_new is the reconstructed image
-
-[SCANNING_GEOMETRY]
-Scanning_Type = FAN_ED              //FAN_ED: Fan-beam with detector units arranged in line; FAN_EA: Fan-beam with units on an circular arc whose origin is the X-ray source; 
-XRay_To_Rotation_Center = 400
-Detector_To_Rotation_Center = 400
-Total_Num_Of_Angles = 22 
-First_Angle_Degree = 0 
-Angle_Increment = 16.36
-Length_Of_Detector = 512            //To FAN_ED type it means the total length of detectors, while to FAN_EA type it means total angular range.
-Num_Of_Detector_Units = 256
-Bias_Of_Detector_Center = 0
-Gaussian_Noise_Percentage = 0.0     //Adds Gaussian noise to detectors
-Store_Projection_Matrix = 1         //For iterative reconstruction ONLY. Will use a lot of memory to store projection matrix, but greatly save iteration time, set 0 to disable.
-
-[IMAGE]
-Rows_Of_Image = 256
-Cloumns_Of_Image = 256
-Pixel_Width = 1
-Pixel_Height = 1
-Bias_Of_Row_Center_To_Rotation_Center = 0
-Bias_Of_Column_Center_To_Rotation_Center = 0
-
-[RECONSTRUCTION]
-Reconstruction_Algorithm = SART	    //SART
-Num_Iterations = 20
-Lambda = 1
-Regularization = 1					//whether or not to use image regularization method in iteration. If to use, the [REGUALRAIZATION] section must be provided in this file.
-
-[REGULARIZATION]
-Regularization_Algorithm = TV		//TV(Total Variation) Minimization
-TV_Iterations = 20                  //iterations of each TV minimization loop
-TV_Alpha = 0.2                      //for TV ONLY
-

release/demo.m

-%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% NOTE: All the examples are for demostration with low iteration time.
-% By chosing a long iteration time one can obtain better image quality.
-% To do this please modify the *.cfg file following the instructions. 
-
-clear all;
-
-I=phantom;
-figure;imshow(I, [0 1]);title('original phantom image');
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% use default.cfg
-P_DEF=fanprj(I);
-I_DEF=fanrec(P_DEF);
-figure;imshow(I_DEF, [0 1]);title('reconstucted image using default config');
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% use other .cfg file
-% ROI problem
-P_ROI=fanprj(I,'config_SART_TV_ROI.cfg');
-I_ROI=fanrec(P_ROI,'config_SART_TV_ROI.cfg');
-figure;imshow(I_ROI, [0 1]);title('reconstucted image for ROI problem');
-
-% few views problem using SART method
-P_FV=fanprj(I,'config_fewviews_SART.cfg');
-I_FV=fanrec(P_FV,'config_fewviews_SART.cfg');
-figure;imshow(I_FV, [0 1]);title('reconstucted image for few views problem');
-
-% Reconstruction with Soft Thresholding Regularization
-P_STH=fanprj(I,'config_SART_STH_YU2009.cfg');
-I_STH=fanrec(P_STH,'config_SART_STH_YU2009.cfg');
-figure;imshow(I_STH, [0 1]);title('Reconstruction with Soft Thresholding Regularization');
-
-% Reconstruction with TV Regularization
-P_TV = fanprj(I,'config_SART_TV_Sidky2006.cfg');
-I_TV = fanrec(P_TV,'config_SART_TV_Sidky2006.cfg');
-figure;imshow(I_TV, [0 1]);title('Reconstruction with TV Regularization');
-

release/fanprj.mexw32

Binary file removed.

release/fanrec.mexw32

Binary file removed.

release/template.cfg

-// This is the template of the config file
-// !!(*&)(!*^@&^$#*^!@*(#&^!)!
-//
-//     TO BE COMPLETED
-//
-// !!!(*!&@(#*&!)@*$^!)*@&#^!
-
-[SCANNING_GEOMETRY]
-XRay_To_Rotation_Center = 400
-Detector_To_Rotation_Center = 400
-Total_Num_Of_Angles = 22 
-First_Angle_Degree = 0 
-Angle_Increment = 16.36
-Detector_Arrangement = LINE         //LINE: detector units are arranged in line; ARC: units are on an circular arc with its origin on the X-ray source.
-Length_Of_Detector = 512            //To LINE arrangement it means the total length of detectors, while to ARC arrangement it means total angular range.
-Num_Of_Detector_Units = 256
-Bias_Of_Detector_Center = 0
-Gaussian_Noise_Percentage = 0.0     //Adds Gaussian noise to detectors
-Store_Projection_Matrix = 1         //For iterative reconstruction ONLY. Will use a lot of memory to store projection matrix, but greatly save iteration time, set 0 to disable.
-
-[IMAGE]
-Rows_Of_Image = 256
-Cloumns_Of_Image = 256
-Pixel_Width = 1
-Pixel_Height = 1
-Bias_Of_Row_Center_To_Rotation_Center = 0
-Bias_Of_Column_Center_To_Rotation_Center = 0
-
-[RECONSTRUCTION]
-Reconstruction_Algorithm = SART	    //SART
-Num_Iterations = 20
-Lambda = 1
-Regularization = 1					//whether or not to use image regularization method in iteration. If to use, the [REGUALRAIZATION] section must be provided in this file.
-
-[REGULARIZATION]
-Regularization_Algorithm = TV		//TV(Total Variation) Minimization
-TV_Iterations = 20                  //iterations of each TV minimization loop
-TV_Alpha = 0.2                      //for TV ONLY
-

testing/testFanBeam.m

 clc;
 clear;
 addpath(genpath('/Users/kritisen/Documents/Coding/utilitiesmatlab/'));    
-addpath(genpath('/Users/kritisen/Documents/Coding/openct/'));
+addpath(genpath('/Users/kritisen/Documents/Coding/openrecon/'));
 % addpath('/Users/kritisen/Documents/Coding/recon/utilsproprietary/');
 
 %% prepare the simulated fanbeam using shepp-logan phantom
 detN = floor(rowN*sqrt(2))+10;
 createPrjCFGFanED(filename, I, scanR, detrR, pxlW, pxlH, ...
     a0, d, length(theta), 'detN', detN);
-P = fanprj(I, filename);
+P = xprj(I, filename);
 
 %% Reconstruction using FBP
 
     'a0', a0, 'dAng', d,  ...
     'rowN', r2, 'colN', c2, ...
     'RecType', 'FBP');
-I10 = fanrec(p3,fname);
+I10 = xrec(p3,fname);
 I10 = I10/2;
 I10 = fliplr(I10');
 I10 = I10/2;
             figure, 
             for i = 1: nSARTIterMax
                 fprintf('iteration %d\t', i);
-                IS = fanrec(P, I0, fname);
-                PS = fanprj(IS, fname);
+                IS = xrec(P, I0, fname);
+                PS = xprj(IS, fname);
                 DelX(i) = mexRMS(IS - I0);
                 RE(i) = mexRMSE(PS, P);
 %                 IS = rot90(IS,2);

testing/testFanBeamCenterShift.m

 clc;
 clear;
 addpath(genpath('/Users/kritisen/Documents/Coding/utilitiesmatlab/'));    
-addpath(genpath('/Users/kritisen/Documents/Coding/openct/'));
+addpath(genpath('/Users/kritisen/Documents/Coding/openrecon/'));
 % addpath('/Users/kritisen/Documents/Coding/recon/utilsproprietary/');
 
 %% prepare the simulated fanbeam using shepp-logan phantom
 detN = floor(rowN*sqrt(2))+10;
 createPrjCFGFanED(filename, I, scanR, detrR, pxlW, pxlH, ...
     a0, d, length(theta), 'detN', detN, 'cs', cs-0.5);
-P = fanprj(I, filename);
+P = xprj(I, filename);
 imshow(P, [])
 
 %% Reconstruction using FBP
     'rowN', r2, 'colN', c2, ...
     'cShift', cs-0.5, ...
     'RecType', 'FBP');
-I10 = fanrec(p3,fname);
+I10 = xrec(p3,fname);
 I10 = I10/2;
 I10 = fliplr(I10');
 I10 = I10/2;
             figure, 
             for i = 1: nSARTIterMax
                 fprintf('iteration %d\t', i);
-                IS = fanrec(P, I0, fname);
-                PS = fanprj(IS, fname);
+                IS = xrec(P, I0, fname);
+                PS = xprj(IS, fname);
                 DelX(i) = mexRMS(IS - I0);
                 RE(i) = mexRMSE(PS, P);
 %                 IS = rot90(IS,2);

testing/testSuites

-openCT-TestSuites
+Notes for openRecon-TestSuites
+[Work in progress]
 
-Metrics
+== Metrics ==
 
 Time
 RMSE
 
 
-Test-cases
+== Test-cases ==
 
 Shepp-logan PARA beam
 
 P = radon(I)
-fanrec('FBP');
-fanrec('SART', 20 iterations, 20 TV iterations)
+xrec('FBP');
+xrec('SART', 20 iterations, 20 TV iterations)
 
+//////////////////////////////////
+//  CT/Micro CT lab
+//  Biomedical Imaging Divison
+//  SBES, Virginia Tech.
+//  Version of 2007.07.24
+//////////////////////////////////
+
+#include <mex.h>
+#include <string.h>
+#include <stdio.h>
+#include <math.h>
+
+#include "img.h"
+#include "prj.h"
+
+int showUsage()
+{
+    mexPrintf("xprj - Fanbeam Prjection Kernel v1.0\n");
+    mexPrintf("\tP = xprj(I) calculates projection P from image I using default configuration (read from file 'default.cfg').\n");
+    mexPrintf("\tP = xprj(I,'CFG_FILENAME') calculates projection P from image I using configuration from file 'CFG_FILENAME' (must be a .cfg file).\n");
+    
+    return 0;
+}
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+    Img_Handle hImg=NULL;
+    Prj_Handle hPrj=NULL;
+    
+    char defcfg[] = "default.cfg";
+    char *def = defcfg;
+    
+    //no any inputs, print usage
+    if(nrhs==0){
+        showUsage();
+        return;
+    }
+
+    //1 input P = xprj(I)
+    if(nrhs==1){
+        if( mxGetClassID(prhs[0]) == mxDOUBLE_CLASS ){
+            mexPrintf("Loading default configuration file 'default.cfg'.\n");
+            
+            hImg = Img_create(def, mxGetPr(prhs[0]));
+            if(hImg==NULL){
+                mexPrintf("Error: Load [IMAGE] config failed. Please check the .cfg file.\n");
+                goto EE;
+            }
+            
+            hPrj = Prj_create(def, 0);
+            if(hPrj==NULL){
+                mexPrintf("Error: Load [SCANNING_GEOMETRY] config failed. Please check the .cfg file.\n\n");
+                goto EE;
+            }
+            
+			//calc prj
+            Prj_calcPrj(hPrj, hImg);
+
+            plhs[0] = mxCreateDoubleMatrix(hPrj->pa.detN,hPrj->pa.angN,mxREAL);
+			memcpy(mxGetPr(plhs[0]), hPrj->prj, hPrj->pa.detN*hPrj->pa.angN*sizeof(double));
+            
+            Prj_delete(hPrj);
+            Img_delete(hImg);
+
+            return;
+        }
+        else{
+            mexPrintf("Error: Invalid input argument.\n");
+            goto EE;
+        }
+    }
+    
+    //2 inputs P = xprj(I,'CFG_FILENAME')
+    if(nrhs==2){
+        if(mxGetClassID(prhs[1]) == mxCHAR_CLASS && mxGetClassID(prhs[0]) == mxDOUBLE_CLASS ){
+            char * c = mxArrayToString(prhs[1]);
+            if(!strstr(c,".cfg")){
+				mexPrintf("Error: Invalid config file.\n");
+				goto EE;
+			}
+
+			hImg = Img_create(c, mxGetPr(prhs[0]));
+			if(hImg==NULL){
+				mexPrintf("Error: Load [IMAGE] config failed. Please check the .cfg file.\n");
+				goto EE;
+			}
+
+			hPrj = Prj_create(c, 0);
+			if(hPrj==NULL){
+				mexPrintf("Error: Load [SCANNING_GEOMETRY] config failed. Please check the .cfg file.\n\n");
+				goto EE;
+			}
+ 
+			//calc prj
+            Prj_calcPrj(hPrj, hImg);
+
+            plhs[0] = mxCreateDoubleMatrix(hPrj->pa.detN,hPrj->pa.angN,mxREAL);
+			memcpy(mxGetPr(plhs[0]), hPrj->prj, hPrj->pa.detN*hPrj->pa.angN*sizeof(double));
+            
+            Prj_delete(hPrj);
+            Img_delete(hImg);
+
+            return;
+
+        }
+		else{
+			mexPrintf("Error: Invalid input arguments.\n");
+            goto EE;
+		}
+                
+    }
+
+
+EE:
+    mexPrintf("Request can not be handled.\n\n");
+    showUsage();    
+
+    if(hPrj)
+		Prj_delete(hPrj);
+	if(hImg)
+		Img_delete(hImg);
+    
+    return;
+}
+//////////////////////////////////
+//  CT/Micro CT lab
+//  Biomedical Imaging Divison
+//  SBES, Virginia Tech.
+//  Version of 2007.07.24
+//////////////////////////////////
+
+#include <mex.h>
+#include <string.h>
+#include <stdio.h>
+#include <math.h>
+
+#include "img.h"
+#include "prj.h"
+#include "rec.h"
+
+int showUsage()
+{
+    mexPrintf("xrec - Fanbeam Reconstruction Kernel v1.0\n");
+    mexPrintf("\tI = xrec(P) reconstructs image I from projections P using default configuration (read from file 'default.cfg').\n");
+    mexPrintf("\tI = xrec(P, 'CFG_FILENAME') reconstructs image I using configuration file 'CFG_FILENAME' (must be a .cfg file).\n");
+    mexPrintf("\tI = xrec(P, I0) same as I=xrec(P) but with a initial image value I0 given.\n");
+    mexPrintf("\tI = xrec(P, I0, 'CFG_FILENAME') given both initial image value I0 and configuration file 'CFG_FILENAME' (must be a .cfg file).\n");
+    
+    return 0;
+}
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+    Img_Handle hImg=NULL;
+    Prj_Handle hPrj=NULL;
+    Rec_Handle hRec=NULL;
+    
+    char defcfg[] = "default.cfg";
+    char * def = defcfg;
+    int prm_set = 0;
+    
+    //no any inputs, print usage
+    if(nrhs==0){
+        showUsage();
+        return;
+    }
+
+    //1 input
+    if(nrhs==1){
+        if( mxGetClassID(prhs[0]) == mxDOUBLE_CLASS ){
+            mexPrintf("Loading default configuration file 'default.cfg'.\n");
+            
+            hImg = Img_create(def,0);
+            if(hImg==NULL){
+                mexPrintf("Error: Load [IMAGE] config failed. Please check the default.cfg file.\n");
+                goto EE;
+            }
+            
+            hPrj = Prj_create(def, mxGetPr(prhs[0]));
+            if(hPrj==NULL){
+                mexPrintf("Error: Load [SCANNING_GEOMETRY] config failed. Please check the default.cfg file.\n\n");
+                goto EE;
+            }
+            
+            hRec = Rec_create(def);
+            if(hRec==NULL){
+                mexPrintf("Error: Load [RECONSTRUCTION] config failed. Please check the default.cfg file.\n");
+                goto EE;
+            }
+             
+            Rec_recImg(hRec, hPrj, hImg);
+
+			plhs[0] = mxCreateDoubleMatrix(hImg->ia.rowN, hImg->ia.colN, mxREAL);
+			memcpy(mxGetPr(plhs[0]), hImg->img, hImg->ia.rowN*hImg->ia.colN*sizeof(double));
+
+            Rec_delete(hRec);
+            Prj_delete(hPrj);
+            Img_delete(hImg);
+
+            return;
+        }
+        else{
+            mexPrintf("Invalid data format. Please make sure the input image data is in double format.\n");
+            goto EE;
+        }
+    }
+    
+    //2 inputs
+    if(nrhs==2){
+        
+        //I = xrec(P, 'CFG_FILENAME')
+        if(mxGetClassID(prhs[1]) == mxCHAR_CLASS && mxGetClassID(prhs[0]) == mxDOUBLE_CLASS){
+
+			//load config from file
+            char * c = mxArrayToString(prhs[1]);
+            if(!strstr(c,".cfg")){
+				mexPrintf("Error: invalid configuration file %s.\n", c);
+				goto EE;
+			}
+
+			mexPrintf("Load config from file %s.\n",c);
+
+			hImg = Img_create(c,0);
+			if(hImg==NULL){
+				mexPrintf("Error: Load [IMAGE] config failed. Please check the .cfg file.\n");
+				goto EE;
+			}
+
+			hPrj = Prj_create(c,mxGetPr(prhs[0]));
+			if(hPrj==NULL){
+				mexPrintf("Error: Load [SCANNING_GEOMETRY] config failed. Please check the .cfg file.\n\n");
+				goto EE;
+			}
+
+			hRec = Rec_create(c);
+			if(hRec==NULL){
+				mexPrintf("Error: Load [RECONSTRUCTION] config failed. Please check the .cfg file.\n");
+				goto EE;
+			}
+
+			//xrec
+            Rec_recImg(hRec, hPrj, hImg);
+
+			plhs[0] = mxCreateDoubleMatrix(hImg->ia.rowN, hImg->ia.colN, mxREAL);
+			memcpy(mxGetPr(plhs[0]), hImg->img, hImg->ia.rowN*hImg->ia.colN*sizeof(double));
+
+			Rec_delete(hRec);
+			Prj_delete(hPrj);
+			Img_delete(hImg);
+
+            return;
+            
+        }
+		//I = xrec(P, I0)
+        else if(mxGetClassID(prhs[1]) == mxDOUBLE_CLASS && mxGetClassID(prhs[0]) == mxDOUBLE_CLASS){
+
+			mexPrintf("Load default config.\n",def);
+
+			hImg = Img_create(def, mxGetPr(prhs[1]));
+			if(hImg==NULL){
+				mexPrintf("Error: Load [IMAGE] config failed. Please check the .cfg file.\n");
+				goto EE;
+			}
+
+			hPrj = Prj_create(def, mxGetPr(prhs[0]));
+			if(hPrj==NULL){
+				mexPrintf("Error: Load [SCANNING_GEOMETRY] config failed. Please check the .cfg file.\n\n");
+				goto EE;
+			}
+
+			hRec = Rec_create(def);
+			if(hRec==NULL){
+				mexPrintf("Error: Load [RECONSTRUCTION] config failed. Please check the .cfg file.\n");
+				goto EE;
+			}
+
+			//xrec
+            Rec_recImg(hRec, hPrj, hImg);
+
+			plhs[0] = mxCreateDoubleMatrix(hImg->ia.rowN, hImg->ia.colN, mxREAL);
+			memcpy(mxGetPr(plhs[0]), hImg->img, hImg->ia.rowN*hImg->ia.colN*sizeof(double));
+
+			Rec_delete(hRec);
+			Prj_delete(hPrj);
+			Img_delete(hImg);
+
+            return;
+
+		}
+		else{
+            mexPrintf("Error: Invalid argument.\n");
+            goto EE;
+        }
+    }
+    
+    //3 inputs
+    if(nrhs==3){ 
+
+        //I = xrec(P, I0, 'CFG_FILENAME')
+        if(mxGetClassID(prhs[0]) == mxDOUBLE_CLASS && mxGetClassID(prhs[1]) == mxDOUBLE_CLASS && mxGetClassID(prhs[2]) == mxCHAR_CLASS){
+
+            char * c = mxArrayToString(prhs[2]);
+            if(!strstr(c,".cfg")){
+				mexPrintf("Error: invalid configuration file %s.\n", c);
+				goto EE;
+			}
+
+			mexPrintf("Load config from file %s.\n",c);mexEvalString("drawnow;");
+
+			hImg = Img_create(c, mxGetPr(prhs[1]));
+			if(hImg==NULL){
+				mexPrintf("Error: Load [IMAGE] config failed. Please check the .cfg file.\n");
+				goto EE;
+			}
+
+			hPrj = Prj_create(c, mxGetPr(prhs[0]));
+			if(hPrj==NULL){
+				mexPrintf("Error: Load [SCANNING_GEOMETRY] config failed. Please check the .cfg file.\n\n");
+				goto EE;
+			}
+
+			hRec = Rec_create(c);
+			if(hRec==NULL){
+				mexPrintf("Error: Load [RECONSTRUCTION] config failed. Please check the .cfg file.\n");
+				goto EE;
+			}
+
+			//xrec
+            Rec_recImg(hRec, hPrj, hImg);
+
+			plhs[0] = mxCreateDoubleMatrix(hImg->ia.rowN, hImg->ia.colN, mxREAL);
+			memcpy(mxGetPr(plhs[0]), hImg->img, hImg->ia.rowN*hImg->ia.colN*sizeof(double));
+
+			Rec_delete(hRec);
+			Prj_delete(hPrj);
+			Img_delete(hImg);
+
+            return;
+		}
+		else{
+			mexPrintf("Error: invalid argument.\n");
+			goto EE;
+		}
+
+	}
+
+    
+EE:
+    mexPrintf("Request can not be handled.\n\n");
+    showUsage();    
+
+    if(hImg) Img_delete(hImg);
+    if(hPrj) Prj_delete(hPrj);
+    if(hRec) Rec_delete(hRec);
+    
+    return;    
+}