diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..159cde9 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +/Output/ +*.m~ diff --git a/Lib/bin/DoCESTScanAtGPU.mexa64 b/Lib/bin/DoCESTScanAtGPU.mexa64 index 40771e3..325ffb5 100755 Binary files a/Lib/bin/DoCESTScanAtGPU.mexa64 and b/Lib/bin/DoCESTScanAtGPU.mexa64 differ diff --git a/Lib/bin/DoCalSARAtCPU.mexmaci64 b/Lib/bin/DoCalSARAtCPU.mexmaci64 new file mode 100755 index 0000000..7f45a33 Binary files /dev/null and b/Lib/bin/DoCalSARAtCPU.mexmaci64 differ diff --git a/Lib/bin/DoGMScanAtGPU.mexa64 b/Lib/bin/DoGMScanAtGPU.mexa64 index 19ceef2..60bbae5 100755 Binary files a/Lib/bin/DoGMScanAtGPU.mexa64 and b/Lib/bin/DoGMScanAtGPU.mexa64 differ diff --git a/Lib/bin/DoGridding.mexmaci64 b/Lib/bin/DoGridding.mexmaci64 new file mode 100755 index 0000000..847bfe7 Binary files /dev/null and b/Lib/bin/DoGridding.mexmaci64 differ diff --git a/Lib/bin/DoMEScanAtGPU.mexa64 b/Lib/bin/DoMEScanAtGPU.mexa64 index 9c8eef6..0deb1cf 100755 Binary files a/Lib/bin/DoMEScanAtGPU.mexa64 and b/Lib/bin/DoMEScanAtGPU.mexa64 differ diff --git a/Lib/bin/DoMTScanAtGPU.mexa64 b/Lib/bin/DoMTScanAtGPU.mexa64 index c304d4e..72e762b 100755 Binary files a/Lib/bin/DoMTScanAtGPU.mexa64 and b/Lib/bin/DoMTScanAtGPU.mexa64 differ diff --git a/Lib/bin/DoSARAverageAtCPU.mexmaci64 b/Lib/bin/DoSARAverageAtCPU.mexmaci64 new file mode 100755 index 0000000..8f56445 Binary files /dev/null and b/Lib/bin/DoSARAverageAtCPU.mexmaci64 differ diff --git a/Lib/bin/DoScanAtCPU.mexmaci64 b/Lib/bin/DoScanAtCPU.mexmaci64 new file mode 100755 index 0000000..7731736 Binary files /dev/null and b/Lib/bin/DoScanAtCPU.mexmaci64 differ diff --git a/Lib/bin/DoScanAtGPU.mexa64 b/Lib/bin/DoScanAtGPU.mexa64 index c969aa6..fbcfcb6 100755 Binary files a/Lib/bin/DoScanAtGPU.mexa64 and b/Lib/bin/DoScanAtGPU.mexa64 differ diff --git a/Lib/src/.gitignore b/Lib/src/.gitignore new file mode 100644 index 0000000..84c048a --- /dev/null +++ b/Lib/src/.gitignore @@ -0,0 +1 @@ +/build/ diff --git a/Lib/src/CMakeLists.txt b/Lib/src/CMakeLists.txt old mode 100644 new mode 100755 index 8dd7321..8044d88 --- a/Lib/src/CMakeLists.txt +++ b/Lib/src/CMakeLists.txt @@ -2,6 +2,33 @@ cmake_minimum_required(VERSION 2.8) project(MRiLab) +# Ubuntu Linux (local user install) +SET(ENV{IPP_ROOT} /home.old/arwillis/intel) +SET(ENV{MATLAB_ROOT} /usr/local/bin/matlab/R2018a) +# On Mac OS X +#SET(ENV{IPP_ROOT} /opt/intel) +#SET(ENV{MATLAB_ROOT} /Applications/MATLAB_R2016b.app) + +SET(CMAKE_VERBOSE_MAKEFILE ON) +#MESSAGE( STATUS "WIN32: " ${WIN32} ) +#MESSAGE( STATUS "UNIX: " ${UNIX} ) +#MESSAGE( STATUS "APPLE: " ${APPLE} ) + +if (APPLE) +SET(LLVM_DIR /usr/local/opt/llvm/share/cmake/modules/) +SET(CMAKE_SKIP_BUILD_RPATH FALSE) +SET(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib") +SET(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE) +SET(MACOSX_RPATH 1) +SET(CC /usr/local/opt/llvm/bin/clang) +SET(CXX /usr/local/opt/llvm/bin/clang++) +SET(CMAKE_CXX_COMPILER g++-7) +endif(APPLE) + + +#set(CC gcc) +#set(CXX g++) + set(CMAKE_INSTALL_PREFIX ${CMAKE_SOURCE_DIR}) set(CMAKE_MODULE_PATH ${CMAKE_SOURCE_DIR}/cmake) set(CMAKE_SIZEOF_VOID_P 8) # use 64-bit by default @@ -13,7 +40,14 @@ add_definitions(/DMX_COMPAT_32) find_package(Matlab REQUIRED) find_package(IPP) find_package(Framewave) -find_package(OpenMP) + +if (APPLE) + set(LDFLAGS -L/usr/local/opt/llvm/lib/) + find_package(LLVM) +else(APPLE) + find_package(OpenMP) +endif(APPLE) + find_package(CUDA) find_package(VTK) find_package(Ismrmrd) diff --git a/Lib/src/cmake/FindIPP.cmake b/Lib/src/cmake/FindIPP.cmake index 5494feb..38966f9 100644 --- a/Lib/src/cmake/FindIPP.cmake +++ b/Lib/src/cmake/FindIPP.cmake @@ -20,6 +20,7 @@ ELSE( "$ENV{IPP_ROOT}" STREQUAL "" ) PATH_SUFFIXES ipp ipp/include) INCLUDE_DIRECTORIES(${IPP_INCLUDE_DIR}) + SET(IPP_LIB_PATH $ENV{IPP_ROOT}/lib/intel64/) FIND_LIBRARY( IPP_S_LIBRARY NAMES "ipps_l" @@ -49,7 +50,7 @@ ELSE( "$ENV{IPP_ROOT}" STREQUAL "" ) INCLUDE_DIRECTORIES(${IPP_INCLUDE_DIR}) SET(IPP_LIB_PATH $ENV{IPP_ROOT}/lib) - + FIND_LIBRARY( IPP_SE_LIBRARY NAMES ippsemergedem64t PATHS ${IPP_LIB_PATH}) @@ -115,49 +116,50 @@ ELSE( "$ENV{IPP_ROOT}" STREQUAL "" ) ) MESSAGE (STATUS "IPP_ROOT (IPP 7): $ENV{IPP_ROOT}") - ELSE ($ENV{IPP_ROOT} MATCHES .*composer.*) #IPP 6.X + ELSE ($ENV{IPP_ROOT} MATCHES .*composer.*) #IPP 2018 UNIX/Mac OS X FIND_PATH(IPP_INCLUDE_DIR ipp.h - PATHS $ENV{IPP_ROOT}/include $ENV{IPP_ROOT}/../include /usr/include /usr/local/include + PATHS $ENV{IPP_ROOT}/ipp/include $ENV{IPP_ROOT}/../include /usr/include /usr/local/include PATH_SUFFIXES ipp ipp/include) INCLUDE_DIRECTORIES(${IPP_INCLUDE_DIR}) - - SET(IPP_LIB_PATH $ENV{IPP_ROOT}/lib) - - FIND_LIBRARY( IPP_SE_LIBRARY - NAMES ippsemergedem64t - PATHS ${IPP_LIB_PATH}) + if (APPLE) + # Mac OS X relative path to IPP libraries + SET(IPP_LIB_PATH $ENV{IPP_ROOT}/ipp/lib) + else(APPLE) + # Linux/Unix relative path to 64-bit IPP libraries + SET(IPP_LIB_PATH $ENV{IPP_ROOT}/ipp/lib/intel64) + endif(APPLE) FIND_LIBRARY( IPP_S_LIBRARY - NAMES ippsmergedem64t_t - PATHS ${IPP_LIB_PATH}) - - FIND_LIBRARY( IPP_VME_LIBRARY - NAMES ippvmemergedem64t + NAMES "ipps" PATHS ${IPP_LIB_PATH}) FIND_LIBRARY( IPP_VM_LIBRARY - NAMES ippvmmergedem64t_t + NAMES "ippvm" PATHS ${IPP_LIB_PATH}) - FIND_LIBRARY( IPP_IOMP5_LIBRARY - NAMES "iomp5" - PATHS $ENV{IPP_ROOT}/sharedlib) - FIND_LIBRARY( IPP_CORE_LIBRARY - NAMES ippcoreem64t_t + NAMES "ippcore" PATHS ${IPP_LIB_PATH}) SET(IPP_LIBRARIES - ${IPP_SE_LIBRARY} ${IPP_S_LIBRARY} - ${IPP_VME_LIBRARY} ${IPP_VM_LIBRARY} - ${IPP_IOMP5_LIBRARY} ${IPP_CORE_LIBRARY} ) + # iomp5 library only exists for Mac OS X + if (APPLE) + FIND_LIBRARY( IPP_IOMP5_LIBRARY + NAMES "iomp5" + PATHS $ENV{IPP_ROOT}/lib) + SET(IPP_LIBRARIES + ${IPP_LIBRARIES} + ${IPP_IOMP5_LIBRARY} + ) + endif(APPLE) + MESSAGE (STATUS "IPP_ROOT (IPP 6): $ENV{IPP_ROOT}") ENDIF ($ENV{IPP_ROOT} MATCHES .*composer.*) @@ -170,4 +172,4 @@ IF(IPP_INCLUDE_DIR AND IPP_LIBRARIES) ENDIF(IPP_INCLUDE_DIR AND IPP_LIBRARIES) INCLUDE( "FindPackageHandleStandardArgs" ) -FIND_PACKAGE_HANDLE_STANDARD_ARGS("IPP" DEFAULT_MSG IPP_LIBRARIES IPP_INCLUDE_DIR) \ No newline at end of file +FIND_PACKAGE_HANDLE_STANDARD_ARGS("IPP" DEFAULT_MSG IPP_LIBRARIES IPP_INCLUDE_DIR) diff --git a/Lib/src/cmake/FindMatlab.cmake b/Lib/src/cmake/FindMatlab.cmake index c88238b..de9532e 100644 --- a/Lib/src/cmake/FindMatlab.cmake +++ b/Lib/src/cmake/FindMatlab.cmake @@ -23,12 +23,14 @@ ELSE("$ENV{MATLAB_ROOT}" STREQUAL "" ) FIND_LIBRARY( MATLAB_MEX_LIBRARY NAMES libmex mex PATHS $ENV{MATLAB_ROOT}/bin $ENV{MATLAB_ROOT}/extern/lib - PATH_SUFFIXES glnxa64 glnx86 win64/microsoft win32/microsoft) + PATH_SUFFIXES glnxa64 glnx86 win64/microsoft win32/microsoft maci64 maci + NO_DEFAULT_PATH) FIND_LIBRARY( MATLAB_MX_LIBRARY NAMES libmx mx - PATHS $ENV{MATLAB_ROOT}/bin $ENV{MATLAB_ROOT}/extern/lib - PATH_SUFFIXES glnxa64 glnx86 win64/microsoft win32/microsoft) + PATHS $ENV{MATLAB_ROOT}/bin $ENV{MATLAB_ROOT}/extern/lib + PATH_SUFFIXES glnxa64 glnx86 win64/microsoft win32/microsoft maci64 maci + NO_DEFAULT_PATH) MESSAGE (STATUS "MATLAB_ROOT: $ENV{MATLAB_ROOT}") diff --git a/Lib/src/engine/CMakeLists.txt b/Lib/src/engine/CMakeLists.txt old mode 100644 new mode 100755 diff --git a/Lib/src/engine/CPUEngine/CMakeLists.txt b/Lib/src/engine/CPUEngine/CMakeLists.txt old mode 100644 new mode 100755 index f3327c4..0a2a0c8 --- a/Lib/src/engine/CPUEngine/CMakeLists.txt +++ b/Lib/src/engine/CPUEngine/CMakeLists.txt @@ -2,7 +2,7 @@ INCLUDE_DIRECTORIES(${MATLAB_INCLUDE_DIR} ${IPP_INCLUDE_DIR} ${FW_INCLUDE_DIR}) SET(CPP_FILES DoScanAtCPU DoCalSARAtCPU) - +SET(OpenMP_CXX_FLAGS -fopenmp) foreach(CPP_FILE ${CPP_FILES}) SET(CPP_FILE_NAME ${CPP_FILE}.cpp) add_library(${CPP_FILE} SHARED ${CPP_FILE_NAME} ${CMAKE_SOURCE_DIR}/Matlabdef.def) @@ -14,12 +14,18 @@ foreach(CPP_FILE ${CPP_FILES}) else(CMAKE_CL_64) SET_TARGET_PROPERTIES(${CPP_FILE} PROPERTIES SUFFIX .mexw32 COMPILE_FLAGS ${OpenMP_CXX_FLAGS} LINK_FLAGS ${OpenMP_CXX_FLAGS}) endif(CMAKE_CL_64) - else(WIN32) + elseif(UNIX AND NOT APPLE) if (CMAKE_SIZEOF_VOID_P MATCHES "8") SET_TARGET_PROPERTIES(${CPP_FILE} PROPERTIES SUFFIX .mexa64 PREFIX "" COMPILE_FLAGS ${OpenMP_CXX_FLAGS} LINK_FLAGS ${OpenMP_CXX_FLAGS}) else(CMAKE_SIZEOF_VOID_P MATCHES "8") SET_TARGET_PROPERTIES(${CPP_FILE} PROPERTIES SUFFIX .mexglx PREFIX "" COMPILE_FLAGS ${OpenMP_CXX_FLAGS} LINK_FLAGS ${OpenMP_CXX_FLAGS}) endif (CMAKE_SIZEOF_VOID_P MATCHES "8") + elseif(APPLE) + if (CMAKE_SIZEOF_VOID_P MATCHES "8") + SET_TARGET_PROPERTIES(${CPP_FILE} PROPERTIES SUFFIX .mexmaci64 PREFIX "" COMPILE_FLAGS ${OpenMP_CXX_FLAGS} LINK_FLAGS ${OpenMP_CXX_FLAGS}) + else(CMAKE_SIZEOF_VOID_P MATCHES "8") + SET_TARGET_PROPERTIES(${CPP_FILE} PROPERTIES SUFFIX .mexmaci PREFIX "" COMPILE_FLAGS ${OpenMP_CXX_FLAGS} LINK_FLAGS ${OpenMP_CXX_FLAGS}) + endif (CMAKE_SIZEOF_VOID_P MATCHES "8") endif(WIN32) install(TARGETS ${CPP_FILE} DESTINATION ../bin) diff --git a/Lib/src/engine/CPUEngine/DoCalSARAtCPU.cpp b/Lib/src/engine/CPUEngine/DoCalSARAtCPU.cpp index d8b9fd8..0fa5bf6 100644 --- a/Lib/src/engine/CPUEngine/DoCalSARAtCPU.cpp +++ b/Lib/src/engine/CPUEngine/DoCalSARAtCPU.cpp @@ -63,7 +63,8 @@ and multi-threading (OpenMP) written by Fang Liu (leoliuf@gmail.com). /* includes CPU kernel */ #include "SARKernel.h" -extern "C" bool mxUnshareArray(mxArray *array_ptr, bool noDeepCopy); +//extern "C" bool mxUnshareArray(mxArray *array_ptr, bool noDeepCopy); +extern "C" int mxUnshareArray(mxArray *array_ptr, int noDeepCopy); /* MEX entry function */ void mexFunction(int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[]) diff --git a/Lib/src/engine/CPUEngine/DoScanAtCPU.cpp b/Lib/src/engine/CPUEngine/DoScanAtCPU.cpp index b1c1f28..a6593fa 100644 --- a/Lib/src/engine/CPUEngine/DoScanAtCPU.cpp +++ b/Lib/src/engine/CPUEngine/DoScanAtCPU.cpp @@ -64,7 +64,8 @@ and multi-threading (OpenMP) written by Fang Liu (leoliuf@gmail.com). /* includes CPU kernel */ #include "BlochKernel.h" -extern "C" bool mxUnshareArray(mxArray *array_ptr, bool noDeepCopy); +//extern "C" bool mxUnshareArray(mxArray *array_ptr, bool noDeepCopy); +extern "C" int mxUnshareArray(mxArray *array_ptr, int noDeepCopy); /* MEX entry function */ void mexFunction(int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[]) diff --git a/Lib/src/engine/GPUEngine/CMakeLists.txt b/Lib/src/engine/GPUEngine/CMakeLists.txt index 4e67050..4513b76 100644 --- a/Lib/src/engine/GPUEngine/CMakeLists.txt +++ b/Lib/src/engine/GPUEngine/CMakeLists.txt @@ -1,18 +1,16 @@ IF (IPP_FOUND) # use ipp # compile kernels for shader models 20 30 35 50, support double-precision floating-point operation - set(CUDA_NVCC_FLAGS -Xcompiler -fPIC -use_fast_math -gencode=arch=compute_20,code="sm_20,compute_20" - -gencode=arch=compute_30,code="sm_30,compute_30" - -gencode=arch=compute_35,code="sm_35,compute_35" - -gencode=arch=compute_50,code="sm_50,compute_50" - --ptxas-options=-v -DMATLAB_MEX_FILE -DIPP) + set(CUDA_NVCC_FLAGS -Xcompiler -fPIC -gencode=arch=compute_30,code="sm_30,compute_30" + -gencode=arch=compute_35,code="sm_35,compute_35" + -gencode=arch=compute_50,code="sm_50,compute_50" + --ptxas-options=-v -DMATLAB_MEX_FILE -DIPP) else (IPP_FOUND) # use framewave # compile kernels for shader models 20 21 30 35 50, support double-precision floating-point operation - set(CUDA_NVCC_FLAGS -Xcompiler -fPIC -use_fast_math -gencode=arch=compute_20,code="sm_20,compute_20" - -gencode=arch=compute_30,code="sm_30,compute_30" - -gencode=arch=compute_35,code="sm_35,compute_35" - -gencode=arch=compute_50,code="sm_50,compute_50" - --ptxas-options=-v -DMATLAB_MEX_FILE -DFW) + set(CUDA_NVCC_FLAGS -Xcompiler -fPIC -gencode=arch=compute_30,code="sm_30,compute_30" + -gencode=arch=compute_35,code="sm_35,compute_35" + -gencode=arch=compute_50,code="sm_50,compute_50" + --ptxas-options=-v -DMATLAB_MEX_FILE -DFW) endif (IPP_FOUND) set(CUDA_NVCC_FLAGS ${CUDA_NVCC_FLAGS} -Xcompiler ${OpenMP_CXX_FLAGS}) @@ -38,11 +36,17 @@ foreach(CU_FILE ${CU_FILES}) else(CMAKE_CL_64) SET_TARGET_PROPERTIES(${CU_FILE} PROPERTIES SUFFIX .mexw32) endif(CMAKE_CL_64) - else(WIN32) + elseif(UNIX AND NOT APPLE) if (CMAKE_SIZEOF_VOID_P MATCHES "8") - SET_TARGET_PROPERTIES(${CU_FILE} PROPERTIES SUFFIX .mexa64 PREFIX "") + SET_TARGET_PROPERTIES(${CU_FILE} PROPERTIES SUFFIX .mexa64 PREFIX "") else(CMAKE_SIZEOF_VOID_P MATCHES "8") - SET_TARGET_PROPERTIES(${CU_FILE} PROPERTIES SUFFIX .mexglx PREFIX "") + SET_TARGET_PROPERTIES(${CU_FILE} PROPERTIES SUFFIX .mexglx PREFIX "") + endif (CMAKE_SIZEOF_VOID_P MATCHES "8") + elseif(APPLE) + if (CMAKE_SIZEOF_VOID_P MATCHES "8") + SET_TARGET_PROPERTIES(${CU_FILE} PROPERTIES SUFFIX .mexmaci64 PREFIX "") + else(CMAKE_SIZEOF_VOID_P MATCHES "8") + SET_TARGET_PROPERTIES(${CU_FILE} PROPERTIES SUFFIX .mexmaci PREFIX "") endif (CMAKE_SIZEOF_VOID_P MATCHES "8") endif(WIN32) diff --git a/Lib/src/engine/GPUEngine/DoCESTScanAtGPU.cu b/Lib/src/engine/GPUEngine/DoCESTScanAtGPU.cu index 0bc7c1f..025911f 100644 --- a/Lib/src/engine/GPUEngine/DoCESTScanAtGPU.cu +++ b/Lib/src/engine/GPUEngine/DoCESTScanAtGPU.cu @@ -70,7 +70,8 @@ /* includes CUDA kernel */ #include "BlochKernelCEST.cuh" -extern "C" bool mxUnshareArray(mxArray *array_ptr, bool noDeepCopy); +//extern "C" bool mxUnshareArray(mxArray *array_ptr, bool noDeepCopy); +extern "C" int mxUnshareArray(mxArray *array_ptr, int noDeepCopy); /* MEX entry function */ void mexFunction(int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[]) diff --git a/Lib/src/engine/GPUEngine/DoGMScanAtGPU.cu b/Lib/src/engine/GPUEngine/DoGMScanAtGPU.cu index 4e593b9..6689df3 100644 --- a/Lib/src/engine/GPUEngine/DoGMScanAtGPU.cu +++ b/Lib/src/engine/GPUEngine/DoGMScanAtGPU.cu @@ -70,7 +70,8 @@ /* includes CUDA kernel */ #include "BlochKernelGM.cuh" -extern "C" bool mxUnshareArray(mxArray *array_ptr, bool noDeepCopy); +//extern "C" bool mxUnshareArray(mxArray *array_ptr, bool noDeepCopy); +extern "C" int mxUnshareArray(mxArray *array_ptr, int noDeepCopy); /* MEX entry function */ void mexFunction(int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[]) diff --git a/Lib/src/engine/GPUEngine/DoMEScanAtGPU.cu b/Lib/src/engine/GPUEngine/DoMEScanAtGPU.cu index 618b40c..14bd25d 100644 --- a/Lib/src/engine/GPUEngine/DoMEScanAtGPU.cu +++ b/Lib/src/engine/GPUEngine/DoMEScanAtGPU.cu @@ -70,7 +70,8 @@ /* includes CUDA kernel */ #include "BlochKernelME.cuh" -extern "C" bool mxUnshareArray(mxArray *array_ptr, bool noDeepCopy); +//extern "C" bool mxUnshareArray(mxArray *array_ptr, bool noDeepCopy); +extern "C" int mxUnshareArray(mxArray *array_ptr, int noDeepCopy); /* MEX entry function */ void mexFunction(int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[]) diff --git a/Lib/src/engine/GPUEngine/DoMTScanAtGPU.cu b/Lib/src/engine/GPUEngine/DoMTScanAtGPU.cu index aa6209b..d36ec11 100644 --- a/Lib/src/engine/GPUEngine/DoMTScanAtGPU.cu +++ b/Lib/src/engine/GPUEngine/DoMTScanAtGPU.cu @@ -70,7 +70,8 @@ /* includes CUDA kernel */ #include "BlochKernelMT.cuh" -extern "C" bool mxUnshareArray(mxArray *array_ptr, bool noDeepCopy); +//extern "C" bool mxUnshareArray(mxArray *array_ptr, bool noDeepCopy); +extern "C" int mxUnshareArray(mxArray *array_ptr, int noDeepCopy); /* MEX entry function */ void mexFunction(int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[]) diff --git a/Lib/src/engine/GPUEngine/DoScanAtGPU.cu b/Lib/src/engine/GPUEngine/DoScanAtGPU.cu index e614bea..2d5db0c 100644 --- a/Lib/src/engine/GPUEngine/DoScanAtGPU.cu +++ b/Lib/src/engine/GPUEngine/DoScanAtGPU.cu @@ -70,7 +70,8 @@ /* includes CUDA kernel */ #include "BlochKernel.cuh" -extern "C" bool mxUnshareArray(mxArray *array_ptr, bool noDeepCopy); +//extern "C" bool mxUnshareArray(mxArray *array_ptr, bool noDeepCopy); +extern "C" int mxUnshareArray(mxArray *array_ptr, int noDeepCopy); /* MEX entry function */ void mexFunction(int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[]) diff --git a/Lib/src/interface/DoMatToHDF5/CMakeLists.txt b/Lib/src/interface/DoMatToHDF5/CMakeLists.txt index 5c8a805..d56519e 100644 --- a/Lib/src/interface/DoMatToHDF5/CMakeLists.txt +++ b/Lib/src/interface/DoMatToHDF5/CMakeLists.txt @@ -11,12 +11,18 @@ if(WIN32) else(CMAKE_CL_64) SET_TARGET_PROPERTIES(DoMatToHDF5 PROPERTIES SUFFIX .mexw32) endif(CMAKE_CL_64) -else(WIN32) +elseif(UNIX AND NOT APPLE) if (CMAKE_SIZEOF_VOID_P MATCHES "8") - SET_TARGET_PROPERTIES(DoMatToHDF5 PROPERTIES SUFFIX .mexa64 PREFIX "") + SET_TARGET_PROPERTIES(DoMatToHDF5 PROPERTIES SUFFIX .mexa64 PREFIX "") else(CMAKE_SIZEOF_VOID_P MATCHES "8") SET_TARGET_PROPERTIES(DoMatToHDF5 PROPERTIES SUFFIX .mexglx PREFIX "") endif (CMAKE_SIZEOF_VOID_P MATCHES "8") +else(APPLE) + if (CMAKE_SIZEOF_VOID_P MATCHES "8") + SET_TARGET_PROPERTIES(DoMatToHDF5 PROPERTIES SUFFIX .mexmaci64 PREFIX "") + else(CMAKE_SIZEOF_VOID_P MATCHES "8") + SET_TARGET_PROPERTIES(DoMatToHDF5 PROPERTIES SUFFIX .mexmaci PREFIX "") + endif (CMAKE_SIZEOF_VOID_P MATCHES "8") endif(WIN32) install(TARGETS DoMatToHDF5 DESTINATION ../bin) \ No newline at end of file diff --git a/Lib/src/recon/DoGridding/CMakeLists.txt b/Lib/src/recon/DoGridding/CMakeLists.txt index a4368f7..23dae83 100644 --- a/Lib/src/recon/DoGridding/CMakeLists.txt +++ b/Lib/src/recon/DoGridding/CMakeLists.txt @@ -10,12 +10,18 @@ if(WIN32) else(CMAKE_CL_64) SET_TARGET_PROPERTIES(DoGridding PROPERTIES SUFFIX .mexw32) endif(CMAKE_CL_64) -else(WIN32) - if (CMAKE_SIZEOF_VOID_P MATCHES "8") - SET_TARGET_PROPERTIES(DoGridding PROPERTIES SUFFIX .mexa64 PREFIX "") - else(CMAKE_SIZEOF_VOID_P MATCHES "8") - SET_TARGET_PROPERTIES(DoGridding PROPERTIES SUFFIX .mexglx PREFIX "") - endif (CMAKE_SIZEOF_VOID_P MATCHES "8") +elseif(UNIX AND NOT APPLE) + if (CMAKE_SIZEOF_VOID_P MATCHES "8") + SET_TARGET_PROPERTIES(DoGridding PROPERTIES SUFFIX .mexa64 PREFIX "") + else(CMAKE_SIZEOF_VOID_P MATCHES "8") + SET_TARGET_PROPERTIES(DoGridding PROPERTIES SUFFIX .mexglx PREFIX "") + endif (CMAKE_SIZEOF_VOID_P MATCHES "8") +elseif(APPLE) + if (CMAKE_SIZEOF_VOID_P MATCHES "8") + SET_TARGET_PROPERTIES(DoGridding PROPERTIES SUFFIX .mexmaci64 PREFIX "") + else(CMAKE_SIZEOF_VOID_P MATCHES "8") + SET_TARGET_PROPERTIES(DoGridding PROPERTIES SUFFIX .mexmaci PREFIX "") + endif (CMAKE_SIZEOF_VOID_P MATCHES "8") endif(WIN32) install(TARGETS DoGridding DESTINATION ../bin) diff --git a/Lib/src/renderer/CMakeLists.txt b/Lib/src/renderer/CMakeLists.txt old mode 100644 new mode 100755 index aee193a..bd9a29a --- a/Lib/src/renderer/CMakeLists.txt +++ b/Lib/src/renderer/CMakeLists.txt @@ -1,5 +1,5 @@ -include(${VTK_USE_FILE}) +#include(${VTK_USE_FILE}) if (VTK_FOUND) add_subdirectory(DoKSpaceTrajVTK) diff --git a/Lib/src/renderer/DoKSpaceTrajVTK/CMakeLists.txt b/Lib/src/renderer/DoKSpaceTrajVTK/CMakeLists.txt index 042f0f6..778bc9d 100644 --- a/Lib/src/renderer/DoKSpaceTrajVTK/CMakeLists.txt +++ b/Lib/src/renderer/DoKSpaceTrajVTK/CMakeLists.txt @@ -1,6 +1,7 @@ add_library(DoKSpaceTrajVTK SHARED DoKSpaceTrajVTK.cpp ${CMAKE_SOURCE_DIR}/Matlabdef.def) target_link_libraries(DoKSpaceTrajVTK ${MATLAB_LIBRARIES}) +include_directories(SYSTEM ${VTK_INCLUDE_DIRS}) if(WIN32) if (CMAKE_CL_64) @@ -8,11 +9,17 @@ if(WIN32) else(CMAKE_CL_64) SET_TARGET_PROPERTIES(DoKSpaceTrajVTK PROPERTIES SUFFIX .mexw32) endif(CMAKE_CL_64) -else(WIN32) +elseif(UNIX AND NOT APPLE) if (CMAKE_SIZEOF_VOID_P MATCHES "8") - SET_TARGET_PROPERTIES(DoKSpaceTrajVTK PROPERTIES SUFFIX .mexa64 PREFIX "") + SET_TARGET_PROPERTIES(DoKSpaceTrajVTK PROPERTIES SUFFIX .mexa64 PREFIX "") else(CMAKE_SIZEOF_VOID_P MATCHES "8") - SET_TARGET_PROPERTIES(DoKSpaceTrajVTK PROPERTIES SUFFIX .mexglx PREFIX "") + SET_TARGET_PROPERTIES(DoKSpaceTrajVTK PROPERTIES SUFFIX .mexglx PREFIX "") + endif (CMAKE_SIZEOF_VOID_P MATCHES "8") +elseif(APPLE) + if (CMAKE_SIZEOF_VOID_P MATCHES "8") + SET_TARGET_PROPERTIES(DoKSpaceTrajVTK PROPERTIES SUFFIX .mexmaci64 PREFIX "") + else(CMAKE_SIZEOF_VOID_P MATCHES "8") + SET_TARGET_PROPERTIES(DoKSpaceTrajVTK PROPERTIES SUFFIX .mexmaci PREFIX "") endif (CMAKE_SIZEOF_VOID_P MATCHES "8") endif(WIN32) diff --git a/Lib/src/utilities/DoSARAverageAtCPU/CMakeLists.txt b/Lib/src/utilities/DoSARAverageAtCPU/CMakeLists.txt index 9cdb276..c5fd638 100644 --- a/Lib/src/utilities/DoSARAverageAtCPU/CMakeLists.txt +++ b/Lib/src/utilities/DoSARAverageAtCPU/CMakeLists.txt @@ -10,12 +10,18 @@ if(WIN32) else(CMAKE_CL_64) SET_TARGET_PROPERTIES(DoSARAverageAtCPU PROPERTIES SUFFIX .mexw32) endif(CMAKE_CL_64) -else(WIN32) +elseif(UNIX AND NOT APPLE) if (CMAKE_SIZEOF_VOID_P MATCHES "8") SET_TARGET_PROPERTIES(DoSARAverageAtCPU PROPERTIES SUFFIX .mexa64 PREFIX "") else(CMAKE_SIZEOF_VOID_P MATCHES "8") SET_TARGET_PROPERTIES(DoSARAverageAtCPU PROPERTIES SUFFIX .mexglx PREFIX "") endif (CMAKE_SIZEOF_VOID_P MATCHES "8") +elseif(APPLE) + if (CMAKE_SIZEOF_VOID_P MATCHES "8") + SET_TARGET_PROPERTIES(DoSARAverageAtCPU PROPERTIES SUFFIX .mexmaci64 PREFIX "") + else(CMAKE_SIZEOF_VOID_P MATCHES "8") + SET_TARGET_PROPERTIES(DoSARAverageAtCPU PROPERTIES SUFFIX .mexmaci PREFIX "") + endif (CMAKE_SIZEOF_VOID_P MATCHES "8") endif(WIN32) install(TARGETS DoSARAverageAtCPU DESTINATION ../bin) diff --git a/MacNotes.txt b/MacNotes.txt new file mode 100644 index 0000000..30335e9 --- /dev/null +++ b/MacNotes.txt @@ -0,0 +1,112 @@ +To compile for Mac OS X +======================== +Download the IPP libraries from Intel (https://software.intel.com/en-us/intel-ipp) +Install basic Xcode command line compilation tools +console:> xcode-select install +Install Home-brew (https://brew.sh) +Install gcc, make and llvm using home-brew +brew install cmake gcc +brew install llvm --with-clang +Change a bunch of CMakeLists.txt ... changes noted in the commit + + +Console output of a MAC OS X build: +=========== +wifi_facstf-10-22-39-94:~ arwillis$ cd /Users/arwillis/Documents/ +wifi_facstf-10-22-39-94:Documents arwillis$ cd MRiLab/Lib/src/ +wifi_facstf-10-22-39-94:src arwillis$ mkdir build +wifi_facstf-10-22-39-94:src arwillis$ cd build/ +wifi_facstf-10-22-39-94:build arwillis$ cmake .. +-- The C compiler identification is AppleClang 9.1.0.9020039 +-- The CXX compiler identification is AppleClang 9.1.0.9020039 +-- Check for working C compiler: /Library/Developer/CommandLineTools/usr/bin/cc +-- Check for working C compiler: /Library/Developer/CommandLineTools/usr/bin/cc -- works +-- Detecting C compiler ABI info +-- Detecting C compiler ABI info - done +-- Detecting C compile features +-- Detecting C compile features - done +-- Check for working CXX compiler: /Library/Developer/CommandLineTools/usr/bin/c++ +-- Check for working CXX compiler: /Library/Developer/CommandLineTools/usr/bin/c++ -- works +-- Detecting CXX compiler ABI info +-- Detecting CXX compiler ABI info - done +-- Detecting CXX compile features +-- Detecting CXX compile features - done +-- MATLAB_ROOT: /Applications/MATLAB_R2016b.app +-- Matlab mex will be used +-- IPP_ROOT (IPP 6): /opt/intel +-- Found IPP: /opt/intel/ipp/lib/libipps.dylib;/opt/intel/ipp/lib/libippvm.dylib;/opt/intel/lib/libiomp5.dylib;/opt/intel/ipp/lib/libippcore.dylib +-- FRAMEWAVE_ROOT environment variable not set. +-- In Linux this can be done in your user .bashrc file by appending the corresponding line, e.g: +-- export FRAMEWAVE_ROOT=/usr/local/framewave/build +-- In Windows this can be done by adding system variable, e.g: +-- FRAMEWAVE_ROOT=C:\Program Files\FRAMEWAVE_1.3.1_SRC\Framewave\build +-- Could NOT find Framewave (missing: FW_LIBRARIES FW_INCLUDE_DIR) +-- IPP and/or Framewave Found, CPU engine will be compiled +-- IPP will be used +CUDA not found. GPU engine will not be compiled +VTK not found. MRiLab DoKSpaceTrajVTK will not be compiled. +ISMRMRD and/or Boost not found. MRiLab to Gadgetron interface will not be compiled +-- Configuring done +CMake Warning (dev): + Policy CMP0042 is not set: MACOSX_RPATH is enabled by default. Run "cmake + --help-policy CMP0042" for policy details. Use the cmake_policy command to + set the policy and suppress this warning. + + MACOSX_RPATH is not specified for the following targets: + + DoCalSARAtCPU + DoGridding + DoSARAverageAtCPU + DoScanAtCPU + +This warning is for project developers. Use -Wno-dev to suppress it. + +-- Generating done +-- Build files have been written to: /Users/arwillis/Documents/MRiLab/Lib/src/build +wifi_facstf-10-22-39-94:build arwillis$ make +/usr/local/Cellar/cmake/3.11.0/bin/cmake -H/Users/arwillis/Documents/MRiLab/Lib/src -B/Users/arwillis/Documents/MRiLab/Lib/src/build --check-build-system CMakeFiles/Makefile.cmake 0 +/usr/local/Cellar/cmake/3.11.0/bin/cmake -E cmake_progress_start /Users/arwillis/Documents/MRiLab/Lib/src/build/CMakeFiles /Users/arwillis/Documents/MRiLab/Lib/src/build/CMakeFiles/progress.marks +/Library/Developer/CommandLineTools/usr/bin/make -f CMakeFiles/Makefile2 all +/Library/Developer/CommandLineTools/usr/bin/make -f engine/CPUEngine/CMakeFiles/DoCalSARAtCPU.dir/build.make engine/CPUEngine/CMakeFiles/DoCalSARAtCPU.dir/depend +cd /Users/arwillis/Documents/MRiLab/Lib/src/build && /usr/local/Cellar/cmake/3.11.0/bin/cmake -E cmake_depends "Unix Makefiles" /Users/arwillis/Documents/MRiLab/Lib/src /Users/arwillis/Documents/MRiLab/Lib/src/engine/CPUEngine /Users/arwillis/Documents/MRiLab/Lib/src/build /Users/arwillis/Documents/MRiLab/Lib/src/build/engine/CPUEngine /Users/arwillis/Documents/MRiLab/Lib/src/build/engine/CPUEngine/CMakeFiles/DoCalSARAtCPU.dir/DependInfo.cmake --color= +Scanning dependencies of target DoCalSARAtCPU +/Library/Developer/CommandLineTools/usr/bin/make -f engine/CPUEngine/CMakeFiles/DoCalSARAtCPU.dir/build.make engine/CPUEngine/CMakeFiles/DoCalSARAtCPU.dir/build +[ 12%] Building CXX object engine/CPUEngine/CMakeFiles/DoCalSARAtCPU.dir/DoCalSARAtCPU.cpp.o +cd /Users/arwillis/Documents/MRiLab/Lib/src/build/engine/CPUEngine && g++-7 -DDoCalSARAtCPU_EXPORTS -DIPP -DMATLAB_MEX_FILE -DMX_COMPAT_32 -I/Applications/MATLAB_R2016b.app/extern/include -I/opt/intel/ipp/include -O3 -DNDEBUG -fPIC -fopenmp -o CMakeFiles/DoCalSARAtCPU.dir/DoCalSARAtCPU.cpp.o -c /Users/arwillis/Documents/MRiLab/Lib/src/engine/CPUEngine/DoCalSARAtCPU.cpp +[ 25%] Linking CXX shared library DoCalSARAtCPU.mexmaci64 +cd /Users/arwillis/Documents/MRiLab/Lib/src/build/engine/CPUEngine && /usr/local/Cellar/cmake/3.11.0/bin/cmake -E cmake_link_script CMakeFiles/DoCalSARAtCPU.dir/link.txt --verbose=1 +g++-7 -O3 -DNDEBUG -dynamiclib -Wl,-headerpad_max_install_names -fopenmp -o DoCalSARAtCPU.mexmaci64 -install_name /Users/arwillis/Documents/MRiLab/Lib/src/build/engine/CPUEngine/DoCalSARAtCPU.mexmaci64 CMakeFiles/DoCalSARAtCPU.dir/DoCalSARAtCPU.cpp.o -Wl,-rpath,/Applications/MATLAB_R2016b.app/bin/maci64 -Wl,-rpath,/opt/intel/ipp/lib -Wl,-rpath,/opt/intel/lib /Applications/MATLAB_R2016b.app/bin/maci64/libmex.dylib /Applications/MATLAB_R2016b.app/bin/maci64/libmx.dylib /opt/intel/ipp/lib/libipps.dylib /opt/intel/ipp/lib/libippvm.dylib /opt/intel/lib/libiomp5.dylib /opt/intel/ipp/lib/libippcore.dylib +[ 25%] Built target DoCalSARAtCPU +/Library/Developer/CommandLineTools/usr/bin/make -f engine/CPUEngine/CMakeFiles/DoScanAtCPU.dir/build.make engine/CPUEngine/CMakeFiles/DoScanAtCPU.dir/depend +cd /Users/arwillis/Documents/MRiLab/Lib/src/build && /usr/local/Cellar/cmake/3.11.0/bin/cmake -E cmake_depends "Unix Makefiles" /Users/arwillis/Documents/MRiLab/Lib/src /Users/arwillis/Documents/MRiLab/Lib/src/engine/CPUEngine /Users/arwillis/Documents/MRiLab/Lib/src/build /Users/arwillis/Documents/MRiLab/Lib/src/build/engine/CPUEngine /Users/arwillis/Documents/MRiLab/Lib/src/build/engine/CPUEngine/CMakeFiles/DoScanAtCPU.dir/DependInfo.cmake --color= +Scanning dependencies of target DoScanAtCPU +/Library/Developer/CommandLineTools/usr/bin/make -f engine/CPUEngine/CMakeFiles/DoScanAtCPU.dir/build.make engine/CPUEngine/CMakeFiles/DoScanAtCPU.dir/build +[ 37%] Building CXX object engine/CPUEngine/CMakeFiles/DoScanAtCPU.dir/DoScanAtCPU.cpp.o +cd /Users/arwillis/Documents/MRiLab/Lib/src/build/engine/CPUEngine && g++-7 -DDoScanAtCPU_EXPORTS -DIPP -DMATLAB_MEX_FILE -DMX_COMPAT_32 -I/Applications/MATLAB_R2016b.app/extern/include -I/opt/intel/ipp/include -O3 -DNDEBUG -fPIC -fopenmp -o CMakeFiles/DoScanAtCPU.dir/DoScanAtCPU.cpp.o -c /Users/arwillis/Documents/MRiLab/Lib/src/engine/CPUEngine/DoScanAtCPU.cpp +[ 50%] Linking CXX shared library DoScanAtCPU.mexmaci64 +cd /Users/arwillis/Documents/MRiLab/Lib/src/build/engine/CPUEngine && /usr/local/Cellar/cmake/3.11.0/bin/cmake -E cmake_link_script CMakeFiles/DoScanAtCPU.dir/link.txt --verbose=1 +g++-7 -O3 -DNDEBUG -dynamiclib -Wl,-headerpad_max_install_names -fopenmp -o DoScanAtCPU.mexmaci64 -install_name /Users/arwillis/Documents/MRiLab/Lib/src/build/engine/CPUEngine/DoScanAtCPU.mexmaci64 CMakeFiles/DoScanAtCPU.dir/DoScanAtCPU.cpp.o -Wl,-rpath,/Applications/MATLAB_R2016b.app/bin/maci64 -Wl,-rpath,/opt/intel/ipp/lib -Wl,-rpath,/opt/intel/lib /Applications/MATLAB_R2016b.app/bin/maci64/libmex.dylib /Applications/MATLAB_R2016b.app/bin/maci64/libmx.dylib /opt/intel/ipp/lib/libipps.dylib /opt/intel/ipp/lib/libippvm.dylib /opt/intel/lib/libiomp5.dylib /opt/intel/ipp/lib/libippcore.dylib +[ 50%] Built target DoScanAtCPU +/Library/Developer/CommandLineTools/usr/bin/make -f recon/DoGridding/CMakeFiles/DoGridding.dir/build.make recon/DoGridding/CMakeFiles/DoGridding.dir/depend +cd /Users/arwillis/Documents/MRiLab/Lib/src/build && /usr/local/Cellar/cmake/3.11.0/bin/cmake -E cmake_depends "Unix Makefiles" /Users/arwillis/Documents/MRiLab/Lib/src /Users/arwillis/Documents/MRiLab/Lib/src/recon/DoGridding /Users/arwillis/Documents/MRiLab/Lib/src/build /Users/arwillis/Documents/MRiLab/Lib/src/build/recon/DoGridding /Users/arwillis/Documents/MRiLab/Lib/src/build/recon/DoGridding/CMakeFiles/DoGridding.dir/DependInfo.cmake --color= +Scanning dependencies of target DoGridding +/Library/Developer/CommandLineTools/usr/bin/make -f recon/DoGridding/CMakeFiles/DoGridding.dir/build.make recon/DoGridding/CMakeFiles/DoGridding.dir/build +[ 62%] Building CXX object recon/DoGridding/CMakeFiles/DoGridding.dir/DoGridding.cpp.o +cd /Users/arwillis/Documents/MRiLab/Lib/src/build/recon/DoGridding && g++-7 -DDoGridding_EXPORTS -DMATLAB_MEX_FILE -DMX_COMPAT_32 -I/Applications/MATLAB_R2016b.app/extern/include -I/opt/intel/ipp/include -O3 -DNDEBUG -fPIC -o CMakeFiles/DoGridding.dir/DoGridding.cpp.o -c /Users/arwillis/Documents/MRiLab/Lib/src/recon/DoGridding/DoGridding.cpp +[ 75%] Linking CXX shared library DoGridding.mexa64 +cd /Users/arwillis/Documents/MRiLab/Lib/src/build/recon/DoGridding && /usr/local/Cellar/cmake/3.11.0/bin/cmake -E cmake_link_script CMakeFiles/DoGridding.dir/link.txt --verbose=1 +g++-7 -O3 -DNDEBUG -dynamiclib -Wl,-headerpad_max_install_names -o DoGridding.mexa64 -install_name /Users/arwillis/Documents/MRiLab/Lib/src/build/recon/DoGridding/DoGridding.mexa64 CMakeFiles/DoGridding.dir/DoGridding.cpp.o -Wl,-rpath,/Applications/MATLAB_R2016b.app/bin/maci64 /Applications/MATLAB_R2016b.app/bin/maci64/libmex.dylib /Applications/MATLAB_R2016b.app/bin/maci64/libmx.dylib +[ 75%] Built target DoGridding +/Library/Developer/CommandLineTools/usr/bin/make -f utilities/DoSARAverageAtCPU/CMakeFiles/DoSARAverageAtCPU.dir/build.make utilities/DoSARAverageAtCPU/CMakeFiles/DoSARAverageAtCPU.dir/depend +cd /Users/arwillis/Documents/MRiLab/Lib/src/build && /usr/local/Cellar/cmake/3.11.0/bin/cmake -E cmake_depends "Unix Makefiles" /Users/arwillis/Documents/MRiLab/Lib/src /Users/arwillis/Documents/MRiLab/Lib/src/utilities/DoSARAverageAtCPU /Users/arwillis/Documents/MRiLab/Lib/src/build /Users/arwillis/Documents/MRiLab/Lib/src/build/utilities/DoSARAverageAtCPU /Users/arwillis/Documents/MRiLab/Lib/src/build/utilities/DoSARAverageAtCPU/CMakeFiles/DoSARAverageAtCPU.dir/DependInfo.cmake --color= +Scanning dependencies of target DoSARAverageAtCPU +/Library/Developer/CommandLineTools/usr/bin/make -f utilities/DoSARAverageAtCPU/CMakeFiles/DoSARAverageAtCPU.dir/build.make utilities/DoSARAverageAtCPU/CMakeFiles/DoSARAverageAtCPU.dir/build +[ 87%] Building CXX object utilities/DoSARAverageAtCPU/CMakeFiles/DoSARAverageAtCPU.dir/DoSARAverageAtCPU.cpp.o +cd /Users/arwillis/Documents/MRiLab/Lib/src/build/utilities/DoSARAverageAtCPU && g++-7 -DDoSARAverageAtCPU_EXPORTS -DMATLAB_MEX_FILE -DMX_COMPAT_32 -I/Applications/MATLAB_R2016b.app/extern/include -I/opt/intel/ipp/include -O3 -DNDEBUG -fPIC -o CMakeFiles/DoSARAverageAtCPU.dir/DoSARAverageAtCPU.cpp.o -c /Users/arwillis/Documents/MRiLab/Lib/src/utilities/DoSARAverageAtCPU/DoSARAverageAtCPU.cpp +[100%] Linking CXX shared library DoSARAverageAtCPU.mexa64 +cd /Users/arwillis/Documents/MRiLab/Lib/src/build/utilities/DoSARAverageAtCPU && /usr/local/Cellar/cmake/3.11.0/bin/cmake -E cmake_link_script CMakeFiles/DoSARAverageAtCPU.dir/link.txt --verbose=1 +g++-7 -O3 -DNDEBUG -dynamiclib -Wl,-headerpad_max_install_names -o DoSARAverageAtCPU.mexa64 -install_name /Users/arwillis/Documents/MRiLab/Lib/src/build/utilities/DoSARAverageAtCPU/DoSARAverageAtCPU.mexa64 CMakeFiles/DoSARAverageAtCPU.dir/DoSARAverageAtCPU.cpp.o -Wl,-rpath,/Applications/MATLAB_R2016b.app/bin/maci64 /Applications/MATLAB_R2016b.app/bin/maci64/libmex.dylib /Applications/MATLAB_R2016b.app/bin/maci64/libmx.dylib +[100%] Built target DoSARAverageAtCPU +/usr/local/Cellar/cmake/3.11.0/bin/cmake -E cmake_progress_start /Users/arwillis/Documents/MRiLab/Lib/src/build/CMakeFiles 0 +wifi_facstf-10-22-39-94:build arwillis$ cp engine/CPUEngine/DoScanAtCPU.mexmaci64 engine/CPUEngine/DoCalSARAtCPU.mexmaci64 ../../bin/ +wifi_facstf-10-22-39-94:build arwillis$ diff --git a/Src/Engine/MATLAB/BlochKernelNormalCPU.m b/Src/Engine/MATLAB/BlochKernelNormalCPU.m new file mode 100644 index 0000000..3de0e07 --- /dev/null +++ b/Src/Engine/MATLAB/BlochKernelNormalCPU.m @@ -0,0 +1,530 @@ +function [Mx,My,Mz] = BlochKernelNormalCPU(Gyro, CS, SpinNum, Rho, T1, T2, Mz, My, Mx, ... + dB0, dWRnd, Gzgrid, Gygrid, Gxgrid, TxCoilmg, TxCoilpe, ... + dt, rfAmp, rfPhase, rfFreq, GzAmp, GyAmp, GxAmp, ... + Typei, SpinMxNum, SpinMxSliceNum, TxCoilNum) +%/* variables for dealing multi-Tx */ +rfAmpSum = 0; +%rfAmpBuf; +%rfPhaseBuf; +%rfFreqBuf; + +%/* Calculate dW */ +buffer = zeros(1, SpinMxNum, 'single'); +buffer1 = zeros(1, SpinMxNum, 'single'); +buffer2 = zeros(1, SpinMxNum, 'single'); +buffer3 = zeros(1, SpinMxNum, 'single'); +buffer4 = zeros(1, SpinMxNum, 'single'); +buffer5 = zeros(1, SpinMxNum, 'single'); +buffer6 = zeros(1, SpinMxNum, 'single'); +SinAlpha = zeros(1, SpinMxNum, 'single'); +SinBeta = zeros(1, SpinMxNum, 'single'); +CosAlpha = zeros(1, SpinMxNum, 'single'); +CosBeta = zeros(1, SpinMxNum, 'single'); +CosPhi = zeros(1, SpinMxNum, 'single'); +SinPhi = zeros(1, SpinMxNum, 'single'); +minusCosPhi = zeros(1, SpinMxNum, 'single'); +minusSinPhi = zeros(1, SpinMxNum, 'single'); +dW = zeros(1, SpinMxNum, 'single'); +alpha = zeros(1, SpinMxNum, 'single'); +beta = zeros(1, SpinMxNum, 'single'); + +% ippsMulC_32f(dB0, Gyro, dW, SpinMxNum); /* add dw caused by main field imhomogeneity */ +dW = dB0 * Gyro; % /* add dw caused by main field imhomogeneity */ + +% ippsAdd_32f_I(dWRnd, dW, SpinMxNum); /* add dWRnd for simulating T2* effect */ +dW = dW + dWRnd; % /* add dWRnd for simulating T2* effect */ +% ippsAddC_32f_I((float)(CS[Typei]) * 2 * PI, dW, SpinMxNum); /* add dw caused by chemical shift */ +dW = dW + CS(Typei+1) * 2 * pi; +% ippsMulC_32f(Gzgrid, (GzAmp) * (Gyro), buffer, SpinMxNum); +buffer = Gzgrid * (GzAmp * Gyro); +% ippsAdd_32f_I(buffer, dW, SpinMxNum); /* add dw caused by Gz for space encoding */ +dW = dW + buffer; % /* add dw caused by Gz for space encoding */ +% ippsMulC_32f(Gygrid, (GyAmp) * (Gyro), buffer, SpinMxNum); +buffer = Gygrid * (GyAmp * Gyro); +% ippsAdd_32f_I(buffer, dW, SpinMxNum); /* add dw caused by Gy for space encoding */ +dW = dW + buffer; % /* add dw caused by Gy for space encoding */ +% ippsMulC_32f(Gxgrid, (GxAmp) * (Gyro), buffer, SpinMxNum); +buffer = Gxgrid * (GxAmp * Gyro); +% ippsAdd_32f_I(buffer, dW, SpinMxNum); /* add dw caused by Gx for space encoding */ +dW = dW + buffer; % /* add dw caused by Gx for space encoding */ + +% /* Spin process */ +for i=0:(TxCoilNum-1) + rfAmpSum=rfAmpSum + abs(rfAmp(i+1)); +end + +if (rfAmpSum ~= 0) + if (TxCoilNum == 1) + rfAmpBuf = rfAmp; + rfPhaseBuf = rfPhase; + rfFreqBuf = rfFreq; % /* note rfFreq is defined as fB0-frf */ + if (rfAmpBuf<0) + rfAmpBuf=abs(rfAmpBuf); + rfPhaseBuf=rfPhaseBuf+pi; + end + + % ippsMulC_32f(TxCoilmg, rfAmpBuf, buffer1, SpinMxNum); /* deal with B1+ magnitude */ + buffer1 = TxCoilmg .* rfAmpBuf; + % ippsAddC_32f(TxCoilpe, rfPhaseBuf, buffer2, SpinMxNum); /* deal with B1+ phase*/ + buffer2 = TxCoilpe + rfPhaseBuf; + % ippsAddC_32f_I(rfFreqBuf * 2 * PI, dW, SpinMxNum); /* add dw caused by off-resonance rf effect */ + dW = dW + rfFreqBuf * 2 * pi; + % ippsMulC_32f_I(Gyro, buffer1, SpinMxNum); + buffer1 = buffer1 * Gyro; + % ippsMulC_32f_I(-1, buffer2, SpinMxNum); + buffer2 = -1 * buffer2; + else + % ippsZero_32f(buffer5, SpinMxNum); + buffer5(:) = 0; + % ippsZero_32f(buffer6, SpinMxNum); + buffer6(:) = 0; + for i=0:(TxCoilNum-1) % /* multi-Tx, sum all (B1+ * rf) */ + rfAmpBuf = rfAmp(i+1); + rfPhaseBuf = rfPhase(i+1); + rfFreqBuf = rfFreq(i+1); % /* note rfFreq is defined as fB0-frf */ + if (rfAmpBuf ~=0 ) + if (rfAmpBuf<0) + rfAmpBuf=abs(rfAmpBuf); + rfPhaseBuf=rfPhaseBuf+pi; + end + + % ippsMulC_32f(TxCoilmg + i * SpinMxSliceNum * SpinMxNum, rfAmpBuf, buffer1, SpinMxNum); /* deal with B1+ magnitude */ + buffer1 = TxCoilmg(i * SpinMxSliceNum * SpinMxNum + 1) .* rfAmpBuf; + % ippsAddC_32f(TxCoilpe + i * SpinMxSliceNum * SpinMxNum, rfPhaseBuf, buffer2, SpinMxNum); /* deal with B1+ phase*/ + buffer2 = TxCoilpe(i * SpinMxSliceNum * SpinMxNum + 1) .* rfPhaseBuf; + % #ifdef IPP + % ippsPolarToCart_32f(buffer1, buffer2, buffer3, buffer4, SpinMxNum); + buffer3 = buffer1 .* cos(buffer2); + buffer4 = buffer1 .* sin(buffer2); + % ippsAdd_32f_I(buffer3, buffer5, SpinMxNum); /* add real from multi-Tx */ + buffer5 = buffer5 + buffer3; + % ippsAdd_32f_I(buffer4, buffer6, SpinMxNum); /* add imag from multi-Tx */ + buffer6 = buffer6 + buffer4; + % ippsAddC_32f_I(rfFreqBuf * 2 * PI, dW, SpinMxNum); /* add dw caused by off-resonance rf effect */ + dW = dW + rfFreqBuf * 2 * pi; + end + end + + % ippsCartToPolar_32f(buffer5, buffer6, buffer1, buffer2, SpinMxNum); /* buffer1 is mag, buffer2 is phase*/ + buffer1 = sqrt(buffer5.^2 + buffer6.^2); + buffer2 = atan2(buffer6, buffer5); + % #endif + + % #ifdef FW + % fwsCos_32f_A24 (buffer2, buffer3, SpinMxNum); + % fwsSin_32f_A24 (buffer2, buffer4, SpinMxNum); + % fwsMul_32f_I(buffer1, buffer3, SpinMxNum); + % fwsMul_32f_I(buffer1, buffer4, SpinMxNum); + % fwsAdd_32f_I(buffer3, buffer5, SpinMxNum); /* add real from multi-Tx */ + % fwsAdd_32f_I(buffer4, buffer6, SpinMxNum); /* add imag from multi-Tx */ + % fwsAddC_32f_I(rfFreqBuf * 2 * PI, dW, SpinMxNum); /* add dw caused by off-resonance rf effect */ + % end + % end + % fwsAtan2_32f_A24(buffer6, buffer5, buffer2, SpinMxNum); /* buffer2 is phase */ + % fwsSqr_32f_I(buffer5, SpinMxNum); + % fwsSqr_32f_I(buffer6, SpinMxNum); + % fwsAdd_32f_I(buffer5, buffer6, SpinMxNum); + % fwsSqrt_32f(buffer6, buffer1, SpinMxNum); /* buffer1 is mag */ + + % #endif + + % ippsMulC_32f_I(Gyro, buffer1, SpinMxNum); + buffer1 = Gyro * buffer1; + % ippsMulC_32f_I(-1, buffer2, SpinMxNum); + buffer2 = -1 * buffer2; + end + + % /* Calculate beta */ + % ippsDiv_32f(buffer1, dW, buffer, SpinMxNum); + buffer = dW ./ buffer1; + % ippsArctan_32f(buffer, beta, SpinMxNum); + beta = atan(buffer); + % /* Calculate alpha */ + % ippsSqr_32f(dW, buffer, SpinMxNum); + buffer = dW.^2; + % ippsSqr_32f_I(buffer1, SpinMxNum); + buffer1 = buffer1.^2; + % ippsAdd_32f_I(buffer1, buffer, SpinMxNum); + buffer = buffer + buffer1; + % ippsSqrt_32f_I(buffer, SpinMxNum); + buffer = sqrt(buffer); + % ippsMulC_32f(buffer, dt, alpha, SpinMxNum); + alpha = buffer * dt; + % /* Trigonometry */ + % ippsSin_32f_A24 (alpha, SinAlpha, SpinMxNum); + SinAlpha = sin(alpha); + % ippsSin_32f_A24 (beta, SinBeta, SpinMxNum); + SinBeta = sin(beta); + % ippsCos_32f_A24 (alpha, CosAlpha, SpinMxNum); + CosAlpha = cos(alpha); + % ippsCos_32f_A24 (beta, CosBeta, SpinMxNum); + CosBeta = cos(beta); + % ippsSin_32f_A24 (buffer2, SinPhi, SpinMxNum); + SinPhi = sin(buffer2); + % ippsMulC_32f(SinPhi, -1, minusSinPhi, SpinMxNum); + minusSinPhi = -1 * SinPhi; + % ippsCos_32f_A24 (buffer2, CosPhi, SpinMxNum); + CosPhi = cos(buffer2); + % ippsMulC_32f(CosPhi, -1, minusCosPhi, SpinMxNum); + minusCosPhi = -1 * CosPhi; + + % /* Calculate Mx */ + + % ippsMul_32f(SinBeta, CosAlpha, buffer1, SpinMxNum); + buffer1 = SinBeta .* CosAlpha; + % ippsMul_32f_I(minusCosPhi, buffer1, SpinMxNum); + buffer1 = buffer1 .* minusCosPhi; + % ippsMul_32f(SinAlpha, SinPhi, buffer3, SpinMxNum); + buffer3 = SinAlpha .* SinPhi; + % ippsAdd_32f_I(buffer3, buffer1, SpinMxNum); + buffer1 = buffer1 + buffer3; + % ippsMul_32f_I(SinBeta, buffer1, SpinMxNum); + buffer1 = buffer1 .* SinBeta; + % ippsMulC_32f_I(-1, buffer1, SpinMxNum); + buffer1 = -1 * buffer1; + % ippsSqr_32f(CosBeta, buffer3, SpinMxNum); + buffer3 = CosBeta.^2; + % ippsMul_32f_I(CosPhi, buffer3, SpinMxNum); + buffer3 = buffer3 .* CosPhi; + % ippsAdd_32f_I(buffer3, buffer1, SpinMxNum); + buffer1 = buffer1 + buffer3; + % ippsCopy_32f(buffer1, buffer3, SpinMxNum); + buffer3 = buffer1; + + % ippsMul_32f_I(CosPhi, buffer1, SpinMxNum); + buffer1 = buffer1 .* CosPhi; + % ippsMul_32f_I(SinPhi, buffer3, SpinMxNum); + buffer3 = buffer3 .* SinPhi; + + % ippsMul_32f(SinAlpha, SinBeta, buffer2, SpinMxNum); + buffer2 = SinAlpha .* SinBeta; + % ippsMul_32f_I(CosPhi, buffer2, SpinMxNum); + buffer2 = buffer2 .* CosPhi; + % ippsMul_32f(CosAlpha,SinPhi,buffer4, SpinMxNum); + buffer4 = CosAlpha .* SinPhi; + % ippsAdd_32f_I(buffer4, buffer2, SpinMxNum); + buffer2 = buffer2 + buffer4; + % ippsCopy_32f(buffer2, buffer4, SpinMxNum); + buffer4 = buffer2; + + % ippsMul_32f_I(SinPhi,buffer2, SpinMxNum); + buffer2 = buffer2 .* SinPhi; + % ippsMul_32f_I(minusCosPhi, buffer4, SpinMxNum); + buffer4 = buffer4 .* minusCosPhi; + + + % ippsAdd_32f_I(buffer1, buffer2, SpinMxNum); + buffer2 = buffer2 + buffer1; + % ippsMul_32f_I(Mx, buffer2, SpinMxNum); + buffer2 = buffer2 .* Mx; + + % ippsAdd_32f_I(buffer3, buffer4, SpinMxNum); + buffer4 = buffer4 + buffer3; + % ippsMul_32f_I(My, buffer4, SpinMxNum); + buffer4 = buffer4 .* My; + % ippsMulC_32f_I(-1, buffer4, SpinMxNum); + buffer4 = -1 * buffer4; + + % ippsAdd_32f(buffer2, buffer4, buffer, SpinMxNum); + buffer = buffer2 + buffer4; + + % ippsMul_32f(SinBeta, CosBeta, buffer1, SpinMxNum); + buffer1 = SinBeta .* CosBeta; + % ippsMul_32f_I(CosPhi, buffer1, SpinMxNum); + buffer1 = buffer1 .* CosPhi; + + % ippsMul_32f(SinBeta, CosAlpha, buffer2, SpinMxNum); + buffer2 = SinBeta .* CosAlpha; + % ippsMul_32f_I(minusCosPhi, buffer2, SpinMxNum); + buffer2 = buffer2 .* minusCosPhi; + % ippsMul_32f(SinAlpha, SinPhi, buffer3, SpinMxNum); + buffer3 = SinAlpha .* SinPhi; + % ippsAdd_32f_I(buffer3, buffer2, SpinMxNum); + buffer2 = buffer2 + buffer3; + % ippsMul_32f_I(CosBeta, buffer2, SpinMxNum); + buffer2 = buffer2 .* CosBeta; + + % ippsAdd_32f_I(buffer2, buffer1, SpinMxNum); + buffer1 = buffer1 + buffer2; + % ippsMul_32f_I(Mz, buffer1, SpinMxNum); + buffer1 = buffer1 .* Mz; + + % ippsAdd_32f(buffer, buffer1, alpha, SpinMxNum); /* use the space of alpha as buffer for storing Mx */ + alpha = buffer + buffer1; + + % /* + % Mx*(cos(PHI)*(cos(BETA)^2*cos(PHI) - sin(BETA)*(sin(ALPHA)*sin(PHI) - cos(ALPHA)*cos(PHI)*sin(BETA))) + sin(PHI)*(cos(ALPHA)*sin(PHI) + cos(PHI)*sin(ALPHA)*sin(BETA))) + % - My*(sin(PHI)*(cos(BETA)^2*cos(PHI) - sin(BETA)*(sin(ALPHA)*sin(PHI) - cos(ALPHA)*cos(PHI)*sin(BETA))) - cos(PHI)*(cos(ALPHA)*sin(PHI) + cos(PHI)*sin(ALPHA)*sin(BETA))) + % + Mz*(cos(BETA)*(sin(ALPHA)*sin(PHI) - cos(ALPHA)*cos(PHI)*sin(BETA)) + cos(BETA)*cos(PHI)*sin(BETA)) + % */ + + % /* Calculate My */ + + + % ippsMul_32f(CosAlpha, SinBeta, buffer1, SpinMxNum); + buffer1 = CosAlpha .* SinBeta; + % ippsMul_32f_I(SinPhi,buffer1, SpinMxNum); + buffer1 = buffer1 .* SinPhi; + % ippsMul_32f(SinAlpha,CosPhi, buffer3, SpinMxNum); + buffer3 = SinAlpha .* CosPhi; + % ippsAdd_32f_I(buffer3, buffer1, SpinMxNum); + buffer1 = buffer1 + buffer3; + % ippsMul_32f_I(SinBeta, buffer1, SpinMxNum); + buffer1 = buffer1 .* SinBeta; + % ippsSqr_32f(CosBeta, buffer3, SpinMxNum); + buffer3 = CosBeta.^2; + % ippsMul_32f_I(SinPhi, buffer3, SpinMxNum); + buffer3 = buffer3 .* SinPhi; + % ippsAdd_32f_I(buffer3, buffer1, SpinMxNum); + buffer1 = buffer1 + buffer3; + % ippsCopy_32f(buffer1, buffer3, SpinMxNum); + buffer3 = buffer1; + + % ippsMul_32f_I(SinPhi, buffer1, SpinMxNum); + buffer1 = buffer1 .* SinPhi; + % ippsMul_32f_I(CosPhi, buffer3, SpinMxNum); + buffer3 = buffer3 .* CosPhi; + + % ippsMul_32f(SinAlpha, SinBeta, buffer2, SpinMxNum); + buffer2 = SinAlpha .* SinBeta; + % ippsMul_32f_I(minusSinPhi, buffer2, SpinMxNum); + buffer2 = buffer2 .* minusSinPhi; + % ippsMul_32f(CosAlpha, CosPhi, buffer4, SpinMxNum); + buffer4 = CosAlpha .* CosPhi; + % ippsAdd_32f_I(buffer4, buffer2, SpinMxNum); + buffer2 = buffer2 + buffer4; + % ippsCopy_32f(buffer2, buffer4, SpinMxNum); + buffer4 = buffer2; + + % ippsMul_32f_I(CosPhi,buffer2, SpinMxNum); + buffer2 = buffer2 .* CosPhi; + % ippsMul_32f_I(minusSinPhi, buffer4, SpinMxNum); + buffer4 = buffer4 .* minusSinPhi; + + % ippsAdd_32f_I(buffer1, buffer2, SpinMxNum); + buffer2 = buffer2 + buffer1; + % ippsMul_32f_I(My, buffer2, SpinMxNum); + buffer2 = buffer2 .* My; + + % ippsAdd_32f_I(buffer3, buffer4, SpinMxNum); + buffer4 = buffer4 + buffer3; + % ippsMul_32f_I(Mx, buffer4, SpinMxNum); + buffer4 = buffer4 .* Mx; + % ippsMulC_32f_I(-1, buffer4, SpinMxNum); + buffer4 = -1 * buffer4; + + % ippsAdd_32f(buffer2, buffer4, buffer, SpinMxNum); + buffer = buffer2 + buffer4; + + % ippsMul_32f(CosBeta, SinBeta, buffer1, SpinMxNum); + buffer1 = CosBeta .* SinBeta; + % ippsMul_32f_I(minusSinPhi, buffer1, SpinMxNum); + buffer1 = buffer1 .* minusSinPhi; + + % ippsMul_32f(CosAlpha, SinBeta, buffer2, SpinMxNum); + buffer2 = CosAlpha .* SinBeta; + % ippsMul_32f_I(SinPhi, buffer2, SpinMxNum); + buffer2 = buffer2 .* SinPhi; + + % ippsMul_32f(SinAlpha, CosPhi, buffer3, SpinMxNum); + buffer3 = SinAlpha .* CosPhi; + + % ippsAdd_32f_I(buffer3, buffer2, SpinMxNum); + buffer2 = buffer2 + buffer3; + % ippsMul_32f_I(CosBeta, buffer2, SpinMxNum); + buffer2 = buffer2 .* CosBeta; + + % ippsAdd_32f_I(buffer2, buffer1, SpinMxNum); + buffer1 = buffer1 + buffer2; + % ippsMul_32f_I(Mz, buffer1, SpinMxNum); + buffer1 = buffer1 .* Mz; + + % ippsAdd_32f(buffer, buffer1, beta, SpinMxNum); /* use the space of beta as buffer for storing My */ + beta = buffer + buffer1; + + % /* + % My*(sin(PHI)*(sin(PHI)*cos(BETA)^2 + sin(BETA)*(cos(PHI)*sin(ALPHA) + cos(ALPHA)*sin(BETA)*sin(PHI))) + cos(PHI)*(cos(ALPHA)*cos(PHI) - sin(ALPHA)*sin(BETA)*sin(PHI))) + % - Mx*(cos(PHI)*(sin(PHI)*cos(BETA)^2 + sin(BETA)*(cos(PHI)*sin(ALPHA) + cos(ALPHA)*sin(BETA)*sin(PHI))) - sin(PHI)*(cos(ALPHA)*cos(PHI) - sin(ALPHA)*sin(BETA)*sin(PHI))) + % + Mz*(cos(BETA)*(cos(PHI)*sin(ALPHA) + cos(ALPHA)*sin(BETA)*sin(PHI)) - cos(BETA)*sin(BETA)*sin(PHI)) + % */ + + % /*Calculate Mz */ + + % ippsMul_32f(CosAlpha, CosBeta, buffer1, SpinMxNum); + buffer1 = CosAlpha .* CosBeta; + % ippsMul_32f_I(SinBeta, buffer1, SpinMxNum); + buffer1 = buffer1 .* SinBeta; + % ippsMulC_32f_I(-1, buffer1, SpinMxNum); + buffer1 = -1 * buffer1; + + % ippsMul_32f(CosBeta, SinBeta, buffer3, SpinMxNum); + buffer3 = CosBeta .* SinBeta; + % ippsAdd_32f_I(buffer3, buffer1, SpinMxNum); + buffer1 = buffer1 + buffer3; + % ippsCopy_32f(buffer1, buffer3, SpinMxNum); + buffer3 = buffer1; + + % ippsMul_32f_I(CosPhi, buffer1, SpinMxNum); + buffer1 = buffer1 .* CosPhi; + % ippsMul_32f_I(SinPhi, buffer3, SpinMxNum); + buffer3 = buffer3 .* SinPhi; + + % ippsMul_32f(CosBeta, SinAlpha, buffer2, SpinMxNum); + buffer2 = CosBeta .* SinAlpha; + % ippsMul_32f_I(minusSinPhi, buffer2, SpinMxNum); + buffer2 = buffer2 .* minusSinPhi; + + % ippsMul_32f(CosBeta, SinAlpha, buffer4, SpinMxNum); + buffer4 = CosBeta .* SinAlpha; + % ippsMul_32f_I(CosPhi, buffer4, SpinMxNum); + buffer4 = buffer4 .* CosPhi; + + % ippsAdd_32f_I(buffer1, buffer2, SpinMxNum); + buffer2 = buffer2 + buffer1; + % ippsMul_32f_I(Mx, buffer2, SpinMxNum); + buffer2 = buffer2 .* Mx; + + % ippsAdd_32f_I(buffer3, buffer4, SpinMxNum); + buffer4 = buffer4 + buffer3; + % ippsMul_32f_I(My, buffer4, SpinMxNum); + buffer4 = buffer4 .* My; + % ippsMulC_32f_I(-1, buffer4, SpinMxNum); + buffer4 = -1 * buffer4; + + % ippsAdd_32f(buffer2, buffer4, buffer, SpinMxNum); + buffer = buffer2 + buffer4; + + % ippsSqr_32f(SinBeta, buffer1, SpinMxNum); + buffer1 = SinBeta.^2; + + % ippsSqr_32f(CosBeta, buffer2, SpinMxNum); + buffer2 = CosBeta.^2; + % ippsMul_32f_I(CosAlpha, buffer2, SpinMxNum); + buffer2 = buffer2 .* CosAlpha; + + % ippsAdd_32f_I(buffer2, buffer1, SpinMxNum); + buffer1 = buffer1 + buffer2; + % ippsMul_32f_I(Mz, buffer1, SpinMxNum); + buffer1 = buffer1 .* Mz; + + % ippsAdd_32f(buffer, buffer1, dW, SpinMxNum); /* use the space of dW as buffer for storing Mz */ + dW = buffer + buffer1; + + % /* + % Mx*(cos(PHI)*(cos(BETA)*sin(BETA) - cos(ALPHA)*cos(BETA)*sin(BETA)) - cos(BETA)*sin(ALPHA)*sin(PHI)) + % - My*(sin(PHI)*(cos(BETA)*sin(BETA) - cos(ALPHA)*cos(BETA)*sin(BETA)) + cos(BETA)*cos(PHI)*sin(ALPHA)) + % + Mz*(cos(ALPHA)*cos(BETA)^2 + sin(BETA)^2) + % */ + + % ippsCopy_32f(alpha, Mx, SpinMxNum); + Mx = alpha; + % ippsCopy_32f(beta, My, SpinMxNum); + My = beta; + % ippsCopy_32f(dW, Mz, SpinMxNum); + Mz = dW; + +else + % /* Calculate alpha */ + % ippsMulC_32f(dW, dt, alpha, SpinMxNum); + alpha = dW * dt; + % ippsSin_32f_A24 (alpha, SinAlpha, SpinMxNum); + SinAlpha = sin(alpha); + % ippsCos_32f_A24 (alpha, CosAlpha, SpinMxNum); + CosAlpha = cos(alpha); + + % /* Calculate Mx */ + % ippsMul_32f(Mx, CosAlpha, buffer1, SpinMxNum); + buffer1 = Mx .* CosAlpha; + % ippsMul_32f(My, SinAlpha, buffer2, SpinMxNum); + buffer2 = My .* SinAlpha; + % ippsAdd_32f(buffer1, buffer2, alpha, SpinMxNum); /* use the space of alpha as buffer for storing Mx */ + alpha = buffer1 + buffer2; + + % /* Calculate My */ + % ippsMul_32f(My, CosAlpha, buffer3, SpinMxNum); + buffer3 = My .* CosAlpha; + % ippsMul_32f(Mx, SinAlpha, buffer4, SpinMxNum); + buffer4 = Mx .* SinAlpha; + % ippsMulC_32f_I(-1, buffer4, SpinMxNum); + buffer4 = -1 * buffer4; + % ippsAdd_32f(buffer3, buffer4, beta, SpinMxNum); /* use the space of beta as buffer for storing My */ + beta = buffer3 + buffer4; + + % /* + % Mx*cos(ALPHA) + My*sin(ALPHA) + % + % My*cos(ALPHA) - Mx*sin(ALPHA) + % + % Mz + % */ + % ippsCopy_32f(alpha, Mx, SpinMxNum); + Mx = alpha; + % ippsCopy_32f(beta, My, SpinMxNum); + My = beta; +end + +% /* Relax */ + +% ippsInv_32f_A24(T2, buffer1, SpinMxNum); +buffer1 = 1./T2; +% ippsMulC_32f_I(dt, buffer1, SpinMxNum); +buffer1 = buffer1 * dt; +% ippsMulC_32f_I(-1, buffer1, SpinMxNum); +buffer1 = -1 * buffer1; +% ippsThreshold_LT_32f_I(buffer1, SpinMxNum, -200);/* Prevent -1.#IND error (i.e. overflow caused by following exp) */ +buffer1(buffer1<-200)=-200; +% ippsExp_32f_I(buffer1, SpinMxNum); +buffer1 = exp(buffer1); +% ippsInv_32f_A24(T1, buffer2, SpinMxNum); +buffer2 = 1./T1; +% ippsMulC_32f_I(dt, buffer2, SpinMxNum); +buffer2 = buffer2 * dt; +% ippsMulC_32f_I(-1, buffer2, SpinMxNum); +buffer2 = -1 * buffer2; +% ippsThreshold_LT_32f_I(buffer2, SpinMxNum, -200);/* Prevent -1.#IND error (i.e. overflow caused by following exp) */ +buffer2(buffer2<-200)=-200; +% ippsExp_32f_I(buffer2, SpinMxNum); +buffer2 = exp(buffer2); +% ippsMul_32f_I(buffer1, Mx, SpinMxNum); +Mx = buffer1 .* Mx; +% ippsMul_32f_I(buffer1, My, SpinMxNum); +My = buffer1 .* My; +% ippsMulC_32f(buffer2, -1, buffer3, SpinMxNum); +buffer3 = -1 * buffer2; +% ippsAddC_32f_I(1, buffer3, SpinMxNum); +buffer3 = buffer3 + 1; +% ippsDivC_32f(Rho, *SpinNum, buffer4, SpinMxNum); +buffer4 = Rho ./ single(SpinNum); +% ippsMul_32f_I(buffer4, buffer3, SpinMxNum); +buffer3 = buffer4 .* buffer3; +% ippsMul_32f(Mz, buffer2, buffer4 ,SpinMxNum); +buffer4 = Mz .* buffer2; +% ippsAdd_32f(buffer3, buffer4, Mz, SpinMxNum); +Mz = buffer3 + buffer4; +% /* +% Mx/exp(dt/T2) +% My/exp(dt/T2) +% Mz/exp(dt/T1) - (Rho*(1/exp(dt/T1) - 1))/SpinNum +% */ + + +% /* Free memory */ +% ippsFree(buffer); +% ippsFree(buffer1); +% ippsFree(buffer2); +% ippsFree(buffer3); +% ippsFree(buffer4); +% ippsFree(buffer5); +% ippsFree(buffer6); +% ippsFree(SinAlpha); +% ippsFree(SinBeta); +% ippsFree(SinPhi); +% ippsFree(minusSinPhi); +% ippsFree(CosAlpha); +% ippsFree(CosBeta); +% ippsFree(CosPhi); +% ippsFree(minusCosPhi); +% ippsFree(dW); +% ippsFree(alpha); +% ippsFree(beta); +end diff --git a/Src/Engine/MATLAB/DoScanAtCPU_MATLAB.m b/Src/Engine/MATLAB/DoScanAtCPU_MATLAB.m new file mode 100644 index 0000000..42070b7 --- /dev/null +++ b/Src/Engine/MATLAB/DoScanAtCPU_MATLAB.m @@ -0,0 +1,582 @@ +function DoScanAtCPU_MATLAB +global VSeq; +global VObj; +global VCtl; +global VMag; +global VCoi; +global VVar; +global VSig; +% const mwSize *SpinMxDims; +% float *MzBase, *MyBase, *MxBase, *Mz, *My, *Mx, *RhoBase, *T1Base, *T2Base, *Rho, *T1, *T2; + +%/* loop control */ +i=int32(0); j=int32(0); s=int32(0); Signali=int32(0); Spini=int32(0); Typei=int32(0); RxCoili=int32(0); TxCoili=int32(0); +Slicei=int32(0); +%mwSize MaxStep, MaxutsStep, MaxrfStep, MaxGzStep, MaxGyStep, MaxGxStep, *SpinNum, *TypeNum, *TxCoilNum, *RxCoilNum, *SignalNum; +flag = zeros(1,6,'double'); +% /* IPP or FW buffer */ +% Ipp32f *buffer1, *buffer2, *buffer3, *buffer4, *buffer; +% /* Function status */ +% int ExtCall; +% /* assign pointers */ +Gyro = VObj.Gyro; +MzBase = VObj.Mz; +MyBase = VObj.My; +MxBase = VObj.Mx; +RhoBase = VObj.Rho; +T1Base = VObj.T1; +T2Base = VObj.T2; +SpinNum = VObj.SpinNum; +TypeNum = VObj.TypeNum; + +dB0Base = VMag.dB0; +dWRndBase = VMag.dWRnd; +GzgridBase = VMag.Gzgrid; +GygridBase = VMag.Gygrid; +GxgridBase = VMag.Gxgrid; + +TxCoilmgBase = VCoi.TxCoilmg; +TxCoilpeBase = VCoi.TxCoilpe; +RxCoilx = VCoi.RxCoilx; +RxCoily = VCoi.RxCoily; +TxCoilNum = VCoi.TxCoilNum; +RxCoilNum = VCoi.RxCoilNum; +TxCoilDefault = VCoi.TxCoilDefault; +RxCoilDefault = VCoi.RxCoilDefault; + +CS = VCtl.CS; +TRNum = VCtl.TRNum; +RunMode = VCtl.RunMode; +MaxThreadNum = VCtl.MaxThreadNum; + +utsLine = VSeq.utsLine; +tsLine = VSeq.tsLine; +rfAmpLine = VSeq.rfAmpLine; +rfPhaseLine = VSeq.rfPhaseLine; +rfFreqLine = VSeq.rfFreqLine; +rfCoilLine = VSeq.rfCoilLine; +GzAmpLine = VSeq.GzAmpLine; +GyAmpLine = VSeq.GyAmpLine; +GxAmpLine = VSeq.GxAmpLine; +ADCLine = VSeq.ADCLine; +ExtLine = VSeq.ExtLine; +flagsLine = VSeq.flagsLine; + +MaxStep = length(VSeq.tsLine); +MaxutsStep = length(VSeq.utsLine); +MaxrfStep = length(VSeq.rfAmpLine); +MaxGzStep = length(VSeq.GzAmpLine); +MaxGyStep = length(VSeq.GyAmpLine); +MaxGxStep = length(VSeq.GxAmpLine); + +t = VVar.t; +dt = VVar.dt; +rfAmp = VVar.rfAmp; +rfPhase = VVar.rfPhase; +rfFreq = VVar.rfFreq; +rfCoil = VVar.rfCoil; +rfRef = VVar.rfRef; +GzAmp = VVar.GzAmp; +GyAmp = VVar.GyAmp; +GxAmp = VVar.GxAmp; +ADC = VVar.ADC; +Ext = VVar.Ext; +KzTmp = VVar.Kz; +KyTmp = VVar.Ky; +KxTmp = VVar.Kx; +utsi = VVar.utsi; +rfi = VVar.rfi; +Gzi = VVar.Gzi; +Gyi = VVar.Gyi; +Gxi = VVar.Gxi; +ADCi = VVar.ADCi; +Exti = VVar.Exti; +TRCount = VVar.TRCount; + +Sy = VSig.Sy; +Sx = VSig.Sx; +Kz = VSig.Kz; +Ky = VSig.Ky; +Kx = VSig.Kx; +MzsBase = VSig.Mz; +MysBase = VSig.My; +MxsBase = VSig.Mx; +Muts = VSig.Muts; +SignalNum = VSig.SignalNum; + +% /* get dimensions of spin matrix */ +SpinMxDimNum = length(size(VObj.Mz));%mxGetNumberOfDimensions(mxGetField(mexGetVariablePtr("global", "VObj"), 0, "Mz")); +% SpinMxDims = (mwSize*) mxCalloc(SpinMxDimNum, sizeof(mwSize)); +SpinMxDims = size(VObj.Mz);%mxGetDimensions(mxGetField(mexGetVariablePtr("global", "VObj"), 0, "Mz")); +SpinMxNum = prod(SpinMxDims(1:2));%SpinMxDims[0] * SpinMxDims[1]; + + + +if (SpinMxDimNum == 2) + SpinMxSliceNum = 1; +else + SpinMxSliceNum = SpinMxDims(3); +end + +% /* assign spins to multi-threads */ +% if (SpinMxSliceNum < *MaxThreadNum) +% ThreadNum = SpinMxSliceNum; +% else +% ThreadNum = *MaxThreadNum; /* Full CPU load */ +% /* Number of parallel threads that the code will try to use (integer). +% Set to 1 if you want to use a single core, set >1 for multithreading. +% This number should be less than the #cpus per machine x #cores/cpu x #thread/core. +% i.e. a machine which two quad-core cpu's could have upto 8 threads. +% Note that if there are fewer images than threads... +% then the code will automatically turn down the number of threads +% (since the extra ones do nothing and waste resources) */ +% +% /* set buffer */ +%buffer1 = zeros(SpinMxSliceNum,SpinMxNum,'single'); +%buffer2 = zeros(SpinMxSliceNum,SpinMxNum,'single'); +%buffer3 = zeros(SpinMxSliceNum,SpinMxNum,'single'); +%buffer4 = zeros(SpinMxSliceNum,SpinMxNum,'single'); +try + buffer1 = zeros(1, SpinMxSliceNum*SpinMxNum, 'single'); + buffer2 = zeros(1, SpinMxSliceNum*SpinMxNum, 'single'); + buffer3 = zeros(1, SpinMxSliceNum*SpinMxNum, 'single'); + buffer4 = zeros(1, SpinMxSliceNum*SpinMxNum, 'single'); + buffer = zeros(1,1,'single'); + % + % /* start simulator execution loop */ + fprintf("TR Counts: %d of %d\n", 1, TRNum); + while (i < MaxStep-1) + % /* check MR sequence pulse flag */ + flag(1:6)=0; + if (tsLine(i+1) ~= tsLine(i+2)) + flag(1)= flag(1) + flagsLine(1,i+1); + flag(2)= flag(2) + flagsLine(2,i+1); + flag(3)= flag(3) + flagsLine(3,i+1); + flag(4)= flag(4) + flagsLine(4,i+1); + flag(5)= flag(5) + flagsLine(5,i+1); + flag(6)= flag(6) + flagsLine(6,i+1); + i = i + 1; + else + flag(1)= flag(1) + flagsLine(1,i+1); + flag(2)= flag(2) + flagsLine(2,i+1); + flag(3)= flag(3) + flagsLine(3,i+1); + flag(4)= flag(4) + flagsLine(4,i+1); + flag(5)= flag(5) + flagsLine(5,i+1); + flag(6)= flag(6) + flagsLine(6,i+1); + + while (tsLine(i+1)==tsLine(i+2)) + flag(1)= flag(1) + flagsLine(1,i+2); + flag(2)= flag(2) + flagsLine(2,i+2); + flag(3)= flag(3) + flagsLine(3,i+2); + flag(4)= flag(4) + flagsLine(4,i+2); + flag(5)= flag(5) + flagsLine(5,i+2); + flag(6)= flag(6) + flagsLine(6,i+2); + i = i + 1; + if (i==MaxStep-1) + break; + end + end + i = i + 1; + end + + % /* update pulse status */ + % *t = *(utsLine + *utsi); + t = utsLine(utsi+1); + VVar.t = t; + % *dt = *(utsLine + (int)fmin(*utsi+1, MaxutsStep-1))-*(utsLine + *utsi); + dt = utsLine(min(utsi+2, MaxutsStep))-utsLine(utsi+1); + VVar.dt = dt; + utsi = min(utsi+1, MaxutsStep); + VVar.utsi = utsi; + + if (flag(1) >= 1) % /* update rfAmp, rfPhase, rfFreq, rfCoil for multiple rf lines*/ + for j=0:(flag(1)-1) + rfCoil = rfCoilLine(rfi+1); + VVar.rfCoil = rfCoil; + TxCoili = int32(rfCoil); + s = rfi + 1; + while (s < MaxrfStep) + if (rfCoil == rfCoilLine(s+1)) + if (abs(rfAmpLine(rfi+1)) <= abs(rfAmpLine(s+1))) + rfAmp(TxCoili)= rfAmpLine(rfi+1); + else + rfAmp(TxCoili)= rfAmpLine(s+1); + end + if (abs(rfPhaseLine(rfi+1)) <= abs(rfPhaseLine(s+1))) + rfPhase(TxCoili)= rfPhaseLine(rfi+1); + else + rfPhase(TxCoili)= rfPhaseLine(s+1); + end + if (abs(rfFreqLine(rfi+1)) <= abs(rfFreqLine(s+1))) + rfFreq(TxCoili)= rfFreqLine(rfi+1); + else + rfFreq(TxCoili)= rfFreqLine(s+1); + end + VVar.rfAmp(TxCoili) = rfAmp(TxCoili); + VVar.rfFreq(TxCoili) = rfFreq(TxCoili); + VVar.rfPhase(TxCoili) = rfPhase(TxCoili); + break; + end + s = s + 1; + end + + rfi = rfi + 1; + VVar.rfi = rfi; + end + end + + if (flag(2)==1)% /*update GzAmp */ + if (abs(GzAmpLine(Gzi+1)) <= abs(GzAmpLine(int32(min(Gzi+2, MaxGzStep))))) + GzAmp = GzAmpLine(Gzi+1); + else + GzAmp = GzAmpLine(Gzi+2); + end + VVar.GzAmp = GzAmp; + Gzi = Gzi + 1; + VVar.Gzi = Gzi; + end + + if (flag(3)==1 )% /*update GyAmp */ + if (abs(GyAmpLine(Gyi+1)) <= abs(GyAmpLine(int32(min(Gyi+2, MaxGyStep))))) + GyAmp = GyAmpLine(Gyi+1); + else + GyAmp = GyAmpLine(Gyi+2); + end + VVar.GyAmp = GyAmp; + Gyi = Gyi + 1; + VVar.Gyi = Gyi; + end + + if (flag(4)==1 )% /*update GxAmp */ + if (abs(GxAmpLine(Gxi+1)) <= abs(GxAmpLine(int32(min(Gxi+2, MaxGxStep))))) + GxAmp = GxAmpLine(Gxi+1); + else + GxAmp = GxAmpLine(Gxi+2); + end + VVar.GxAmp = GxAmp; + Gxi = Gxi + 1; + VVar.Gxi = Gxi; + end + + ADC = 0; %/* prevent ADC overflow*/ + VVar.ADC = ADC; + if (flag(5)==1)% /* update ADC */ + ADC = ADCLine(ADCi+1); + VVar.ADC = ADC; + ADCi = ADCi + 1; + VVar.ADCi = ADCi; + end + + if (flag(6)==1)% /* update Ext */ + Ext = ExtLine(Exti+1); + VVar.Ext = Ext; + %/* execute extended process */ + if (Ext ~= 0) + %ExtCall = mexEvalString("DoExtPlugin"); + %ExtCall = DoExtPlugin; + %disp("DoExtPlugin: ERROR NOT IMPLEMENTED"); + ExtCall = false; + if (ExtCall) + disp("Extended process encounters ERROR!"); + return; + end + + % /* update pointers, avoid pointer change between Matlab and Mex call */ + % MzBase = (float*) mxGetData(mxGetField(mexGetVariablePtr("global", "VObj"), 0, "Mz")); + % MyBase = (float*) mxGetData(mxGetField(mexGetVariablePtr("global", "VObj"), 0, "My")); + % MxBase = (float*) mxGetData(mxGetField(mexGetVariablePtr("global", "VObj"), 0, "Mx")); + % RhoBase = (float*) mxGetData(mxGetField(mexGetVariablePtr("global", "VObj"), 0, "Rho")); + % T1Base = (float*) mxGetData(mxGetField(mexGetVariablePtr("global", "VObj"), 0, "T1")); + % T2Base = (float*) mxGetData(mxGetField(mexGetVariablePtr("global", "VObj"), 0, "T2")); + % dWRndBase = (float*) mxGetData(mxGetField(mexGetVariablePtr("global", "VMag"), 0, "dWRnd")); + % dB0Base = (float*) mxGetData(mxGetField(mexGetVariablePtr("global", "VMag"), 0, "dB0")); + % GzgridBase = (float*) mxGetData(mxGetField(mexGetVariablePtr("global", "VMag"), 0, "Gzgrid")); + % GygridBase = (float*) mxGetData(mxGetField(mexGetVariablePtr("global", "VMag"), 0, "Gygrid")); + % GxgridBase = (float*) mxGetData(mxGetField(mexGetVariablePtr("global", "VMag"), 0, "Gxgrid")); + MzBase = VObj.Mz; + MyBase = VObj.My; + MxBase = VObj.Mx; + RhoBase = VObj.Rho; + T1Base = VObj.T1; + T2Base = VObj.T2; + dB0Base = VMag.dB0; + dWRndBase = VMag.dWRnd; + GzgridBase = VMag.Gzgrid; + GygridBase = VMag.Gygrid; + GxgridBase = VMag.Gxgrid; + + % TxCoilmgBase = (float*) mxGetData(mxGetField(mexGetVariablePtr("global", "VCoi"), 0, "TxCoilmg")); + % TxCoilpeBase = (float*) mxGetData(mxGetField(mexGetVariablePtr("global", "VCoi"), 0, "TxCoilpe")); + % RxCoilx = (float*) mxGetData(mxGetField(mexGetVariablePtr("global", "VCoi"), 0, "RxCoilx")); + % RxCoily = (float*) mxGetData(mxGetField(mexGetVariablePtr("global", "VCoi"), 0, "RxCoily")); + TxCoilmgBase = VCoi.TxCoilmg; + TxCoilpeBase = VCoi.TxCoilpe; + RxCoilx = VCoi.RxCoilx; + RxCoily = VCoi.RxCoily; + + % t = (double*) mxGetData(mxGetField(mexGetVariablePtr("global", "VVar"), 0, "t")); + % dt = (double*) mxGetData(mxGetField(mexGetVariablePtr("global", "VVar"), 0, "dt")); + % rfAmp = (double*) mxGetData(mxGetField(mexGetVariablePtr("global", "VVar"), 0, "rfAmp")); + % rfPhase = (double*) mxGetData(mxGetField(mexGetVariablePtr("global", "VVar"), 0, "rfPhase")); + % rfFreq = (double*) mxGetData(mxGetField(mexGetVariablePtr("global", "VVar"), 0, "rfFreq")); + % rfCoil = (double*) mxGetData(mxGetField(mexGetVariablePtr("global", "VVar"), 0, "rfCoil")); + % rfRef = (double*) mxGetData(mxGetField(mexGetVariablePtr("global", "VVar"), 0, "rfRef")); + % GzAmp = (double*) mxGetData(mxGetField(mexGetVariablePtr("global", "VVar"), 0, "GzAmp")); + % GyAmp = (double*) mxGetData(mxGetField(mexGetVariablePtr("global", "VVar"), 0, "GyAmp")); + % GxAmp = (double*) mxGetData(mxGetField(mexGetVariablePtr("global", "VVar"), 0, "GxAmp")); + % ADC = (double*) mxGetData(mxGetField(mexGetVariablePtr("global", "VVar"), 0, "ADC")); + % Ext = (double*) mxGetData(mxGetField(mexGetVariablePtr("global", "VVar"), 0, "Ext")); + % KzTmp = (double*) mxGetData(mxGetField(mexGetVariablePtr("global", "VVar"), 0, "Kz")); + % KyTmp = (double*) mxGetData(mxGetField(mexGetVariablePtr("global", "VVar"), 0, "Ky")); + % KxTmp = (double*) mxGetData(mxGetField(mexGetVariablePtr("global", "VVar"), 0, "Kx")); + % utsi = (int*) mxGetData(mxGetField(mexGetVariablePtr("global", "VVar"), 0, "utsi")); + % rfi = (int*) mxGetData(mxGetField(mexGetVariablePtr("global", "VVar"), 0, "rfi")); + % Gzi = (int*) mxGetData(mxGetField(mexGetVariablePtr("global", "VVar"), 0, "Gzi")); + % Gyi = (int*) mxGetData(mxGetField(mexGetVariablePtr("global", "VVar"), 0, "Gyi")); + % Gxi = (int*) mxGetData(mxGetField(mexGetVariablePtr("global", "VVar"), 0, "Gxi")); + % ADCi = (int*) mxGetData(mxGetField(mexGetVariablePtr("global", "VVar"), 0, "ADCi")); + % Exti = (int*) mxGetData(mxGetField(mexGetVariablePtr("global", "VVar"), 0, "Exti")); + % TRCount = (int*) mxGetData(mxGetField(mexGetVariablePtr("global", "VVar"), 0, "TRCount")); + t = VVar.t; + dt = VVar.dt; + rfAmp = VVar.rfAmp; + rfPhase = VVar.rfPhase; + rfFreq = VVar.rfFreq; + rfCoil = VVar.rfCoil; + rfRef = VVar.rfRef; + GzAmp = VVar.GzAmp; + GyAmp = VVar.GyAmp; + GxAmp = VVar.GxAmp; + ADC = VVar.ADC; + Ext = VVar.Ext; + KzTmp = VVar.Kz(1); + KyTmp = VVar.Ky(1); + KxTmp = VVar.Kx(1); + utsi = VVar.utsi; + rfi = VVar.rfi; + Gzi = VVar.Gzi; + Gyi = VVar.Gyi; + Gxi = VVar.Gxi; + ADCi = VVar.ADCi; + Exti = VVar.Exti; + TRCount = VVar.TRCount; + + end + Exti = Exti + 1; + VVar.Exti = Exti; + end + + if (flag(1)+flag(2)+flag(3)+flag(4)+flag(5)+flag(6) == 0)% /* reset VVar */ + rfAmp(1:TxCoilNum)=0; + rfPhase(1:TxCoilNum)=0; + rfFreq(1:TxCoilNum)=0; + %ippsZero_64f(rfAmp, *TxCoilNum); + %ippsZero_64f(rfPhase, *TxCoilNum); + %ippsZero_64f(rfFreq, *TxCoilNum); + GzAmp = 0; + GyAmp = 0; + GxAmp = 0; + ADC = 0; + Ext = 0; + VVar.GzAmp = GzAmp; + VVar.GyAmp = GyAmp; + VVar.GxAmp = GxAmp; + VVar.ADC = ADC; + VVar.Ext = Ext; + end + + %/* execute spin precessing */ + if (dt == 0)% /* end of time point */ + continue; + elseif (dt < 0)% /* uncontinuous time point process */ + TRCount = TRCount + 1; + VVar.TRCount = TRCount; + fprintf("TR Counts: %d of %d\n", TRCount, TRNum); + + switch (RunMode) + case 1 %/* spin rotation simulation & rf simulation */ + disp("RunMode 1: NOT IMPLEMENTED!"); + %MzsBase = MzBase; + %MysBase = MyBase; + %MxsBase = MxBase; + %ippsCopy_32f(MzBase, MzsBase+(*utsi)*(*TypeNum)*(*SpinNum)*SpinMxSliceNum*SpinMxNum, (*TypeNum)*(*SpinNum)*SpinMxSliceNum*SpinMxNum); + %ippsCopy_32f(MyBase, MysBase+(*utsi)*(*TypeNum)*(*SpinNum)*SpinMxSliceNum*SpinMxNum, (*TypeNum)*(*SpinNum)*SpinMxSliceNum*SpinMxNum); + %ippsCopy_32f(MxBase, MxsBase+(*utsi)*(*TypeNum)*(*SpinNum)*SpinMxSliceNum*SpinMxNum, (*TypeNum)*(*SpinNum)*SpinMxSliceNum*SpinMxNum); + Muts(utsi+1) = Muts(utsi); + break; + end + continue; + end + + for Typei = 0:(TypeNum-1) + for Spini = 0:(SpinNum-1) + %/* signal acquisition */ + if (ADC == 1) + % SpinMxNum = num pixels (x,y) + % SpinMxSliceNum = num slices (3) + % SpinNum = 1?, Typei = 0? + %Mx = MxBase + (Typei*(*SpinNum)*SpinMxSliceNum*SpinMxNum + Spini*SpinMxSliceNum*SpinMxNum); + %My = MyBase + (Typei*(*SpinNum)*SpinMxSliceNum*SpinMxNum + Spini*SpinMxSliceNum*SpinMxNum); + offset = (Typei*SpinNum*SpinMxSliceNum*SpinMxNum + Spini*SpinMxSliceNum*SpinMxNum+1); + Mx = MxBase(offset:(offset+SpinMxSliceNum*SpinMxNum-1)); + My = MyBase(offset:(offset+SpinMxSliceNum*SpinMxNum-1)); + for RxCoili = 0:(RxCoilNum-1)% /* signal acquisition per Rx coil */ + %/* RxCoil sensitivity */ + if (RxCoilDefault == 0) + disp("RxCoilDefault==0: NOT IMPLEMENTED"); + %ippsMul_32f(Mx, RxCoilx+RxCoili*SpinMxSliceNum*SpinMxNum, buffer1, SpinMxSliceNum*SpinMxNum); + %ippsMul_32f(My, RxCoily+RxCoili*SpinMxSliceNum*SpinMxNum, buffer2, SpinMxSliceNum*SpinMxNum); + %ippsAdd_32f(buffer1, buffer2, buffer3, SpinMxSliceNum*SpinMxNum); + + %ippsMul_32f(Mx, RxCoily+RxCoili*SpinMxSliceNum*SpinMxNum, buffer1, SpinMxSliceNum*SpinMxNum); + %ippsMulC_32f_I(-1, buffer1, SpinMxSliceNum*SpinMxNum); + %ippsMul_32f(My, RxCoilx+RxCoili*SpinMxSliceNum*SpinMxNum, buffer2, SpinMxSliceNum*SpinMxNum); + %ippsAdd_32f(buffer1, buffer2, buffer4, SpinMxSliceNum*SpinMxNum); + + %ippsCopy_32f(buffer3, buffer1, SpinMxSliceNum*SpinMxNum); + %ippsCopy_32f(buffer4, buffer2, SpinMxSliceNum*SpinMxNum); + else + %ippsCopy_32f(Mx, buffer1, SpinMxSliceNum*SpinMxNum); + %ippsCopy_32f(My, buffer2, SpinMxSliceNum*SpinMxNum); + %ippsCopy_32f(Mx, buffer3, SpinMxSliceNum*SpinMxNum); + %ippsCopy_32f(My, buffer4, SpinMxSliceNum*SpinMxNum); + buffer1 = Mx(1:SpinMxSliceNum*SpinMxNum); + buffer2 = My(1:SpinMxSliceNum*SpinMxNum); + buffer3 = Mx(1:SpinMxSliceNum*SpinMxNum); + buffer4 = My(1:SpinMxSliceNum*SpinMxNum); + end + + %/* rfRef for demodulating rf Phase */ + if (rfRef ~= 0) + %ippsMulC_32f_I((float)cos(-(*rfRef)), buffer1, SpinMxSliceNum*SpinMxNum); + buffer1 = buffer1 .* cos(-rfRef); + %ippsMulC_32f_I((float)-sin(-(*rfRef)), buffer2, SpinMxSliceNum*SpinMxNum); + buffer2 = buffer2 .* -sin(-rfRef); + %ippsAdd_32f_I(buffer2, buffer1, SpinMxSliceNum*SpinMxNum); + buffer1 = buffer1 + buffer2; + + %ippsMulC_32f_I((float)sin(-(*rfRef)), buffer3, SpinMxSliceNum*SpinMxNum); + buffer3 = buffer3 .* sin(-rfRef); + %ippsMulC_32f_I((float)cos(-(*rfRef)), buffer4, SpinMxSliceNum*SpinMxNum); + buffer4 = buffer4 .* cos(-rfRef); + %ippsAdd_32f_I(buffer4, buffer3, SpinMxSliceNum*SpinMxNum); + buffer3 = buffer3 + buffer4; + else + %ippsCopy_32f(buffer4, buffer3, SpinMxSliceNum*SpinMxNum); + buffer3 = buffer4; + end + + %/* signal summation */ + %ippsSum_32f(buffer1, SpinMxSliceNum*SpinMxNum, buffer, ippAlgHintFast); + buffer = sum(buffer1); + %Sx[Typei*(*RxCoilNum)*(*SignalNum)+RxCoili*(*SignalNum)+Signali] += (double)*buffer; + Sx(Typei*RxCoilNum*SignalNum+RxCoili*SignalNum+Signali+1) = ... + Sx(Typei*RxCoilNum*SignalNum+RxCoili*SignalNum+Signali+1) + buffer; + VSig.Sx(Typei*RxCoilNum*SignalNum+RxCoili*SignalNum+Signali+1) = ... + Sx(Typei*RxCoilNum*SignalNum+RxCoili*SignalNum+Signali+1); + %ippsSum_32f(buffer3, SpinMxSliceNum*SpinMxNum, buffer, ippAlgHintFast); + buffer = sum(buffer3); + %Sy[Typei*(*RxCoilNum)*(*SignalNum)+RxCoili*(*SignalNum)+Signali] += (double)*buffer; + Sy(Typei*RxCoilNum*SignalNum+RxCoili*SignalNum+Signali+1) = ... + Sy(Typei*RxCoilNum*SignalNum+RxCoili*SignalNum+Signali+1) + buffer; + VSig.Sy(Typei*RxCoilNum*SignalNum+RxCoili*SignalNum+Signali+1) = ... + Sy(Typei*RxCoilNum*SignalNum+RxCoili*SignalNum+Signali+1); + end + end + + %/* set openMP for core inner loop */ + + %#pragma omp parallel num_threads(ThreadNum){ + + %#pragma omp for private(Slicei, Mz, My, Mx, dWRnd, Rho, T1, T2, dB0, Gzgrid, Gygrid, Gxgrid, TxCoilmg, TxCoilpe) + %/* need private clause to keep variables isolated, VERY IMPORTANT!!! Otherwise cause cross influence by threads */ + + for Slicei=0:(SpinMxSliceNum-1) + + %/* don't put Matlab function here, may cause Matlab crash! */ + + %/* set pointer offset for input */ + %Mz = MzBase + (Typei*(*SpinNum)*SpinMxSliceNum*SpinMxNum + Spini*SpinMxSliceNum*SpinMxNum + Slicei*SpinMxNum); + %My = MyBase + (Typei*(*SpinNum)*SpinMxSliceNum*SpinMxNum + Spini*SpinMxSliceNum*SpinMxNum + Slicei*SpinMxNum); + %Mx = MxBase + (Typei*(*SpinNum)*SpinMxSliceNum*SpinMxNum + Spini*SpinMxSliceNum*SpinMxNum + Slicei*SpinMxNum); + %dWRnd = dWRndBase + (Typei*(*SpinNum)*SpinMxSliceNum*SpinMxNum + Spini*SpinMxSliceNum*SpinMxNum + Slicei*SpinMxNum); + offset = (Typei*SpinNum*SpinMxSliceNum*SpinMxNum + Spini*SpinMxSliceNum*SpinMxNum + Slicei*SpinMxNum+1); + Mz = MzBase(offset:(offset+SpinMxNum-1)); + My = MyBase(offset:(offset+SpinMxNum-1)); + Mx = MxBase(offset:(offset+SpinMxNum-1)); + dWRnd = dWRndBase(offset:(offset+SpinMxNum-1)); + + %Rho = RhoBase+(Typei*SpinMxSliceNum*SpinMxNum + Slicei*SpinMxNum); + %T1 = T1Base+ (Typei*SpinMxSliceNum*SpinMxNum + Slicei*SpinMxNum); + %T2 = T2Base+ (Typei*SpinMxSliceNum*SpinMxNum + Slicei*SpinMxNum); + offset = (Typei*SpinMxSliceNum*SpinMxNum + Slicei*SpinMxNum+1); + Rho = RhoBase(offset:(offset+SpinMxNum-1)); + T1 = T1Base(offset:(offset+SpinMxNum-1)); + T2 = T2Base(offset:(offset+SpinMxNum-1)); + + %dB0 = dB0Base + Slicei*SpinMxNum; + %Gzgrid = GzgridBase + Slicei*SpinMxNum; + %Gygrid = GygridBase + Slicei*SpinMxNum; + %Gxgrid = GxgridBase + Slicei*SpinMxNum; + %TxCoilmg = TxCoilmgBase + Slicei*SpinMxNum; + %TxCoilpe = TxCoilpeBase + Slicei*SpinMxNum; + offset = Slicei*SpinMxNum+1; + dB0 = dB0Base(offset:(offset+SpinMxNum-1)); + Gzgrid = GzgridBase(offset:(offset+SpinMxNum-1)); + Gygrid = GygridBase(offset:(offset+SpinMxNum-1)); + Gxgrid = GxgridBase(offset:(offset+SpinMxNum-1)); + TxCoilmg = TxCoilmgBase(offset:(offset+SpinMxNum-1)); + TxCoilpe = TxCoilpeBase(offset:(offset+SpinMxNum-1)); + + % /* call spin discrete precessing */ + [Mx, My, Mz] = BlochKernelNormalCPU(Gyro, CS, SpinNum, Rho, T1, T2, Mz, My, Mx, ... + dB0, dWRnd, Gzgrid, Gygrid, Gxgrid, TxCoilmg, TxCoilpe, ... + dt, rfAmp, rfPhase, rfFreq, GzAmp, GyAmp, GxAmp, ... + Typei, SpinMxNum, SpinMxSliceNum, TxCoilNum); + MxBase(offset:(offset+SpinMxNum-1)) = Mx; + MyBase(offset:(offset+SpinMxNum-1)) = My; + MzBase(offset:(offset+SpinMxNum-1)) = Mz; + VObj.Mx(offset:(offset+SpinMxNum-1)) = Mx; + VObj.My(offset:(offset+SpinMxNum-1)) = My; + VObj.Mz(offset:(offset+SpinMxNum-1)) = Mz; + if (false && Slicei==0) + figure(2), imshow(abs(reshape(Mx,[100,80])),[]); + drawnow; + end + end + end + end + + if (ADC == 1) + Kz(Signali+1) = Kz(Signali+1) + KzTmp; + Ky(Signali+1) = Ky(Signali+1) + KyTmp; + Kx(Signali+1) = Kx(Signali+1) + KxTmp; + VSig.Kz(Signali+1) = Kz(Signali+1); + VSig.Ky(Signali+1) = Ky(Signali+1); + VSig.Kx(Signali+1) = Kx(Signali+1); + Signali = Signali + 1; + end + + %/* update Kz, Ky & Kx */ + KzTmp = KzTmp + (GzAmp)*(dt)*(Gyro/(2*pi)); + KyTmp = KyTmp + (GyAmp)*(dt)*(Gyro/(2*pi)); + KxTmp = KxTmp + (GxAmp)*(dt)*(Gyro/(2*pi)); + VVar.KzTmp = KzTmp; + VVar.KyTmp = KyTmp; + VVar.KxTmp = KxTmp; + + switch (RunMode) + case 1 %/* spin rotation simulation & rf simulation */ + disp("RunMode 1: NOT IMPLEMENTED!"); + %ippsCopy_32f(MzBase, MzsBase+(*utsi)*(*TypeNum)*(*SpinNum)*SpinMxSliceNum*SpinMxNum, (*TypeNum)*(*SpinNum)*SpinMxSliceNum*SpinMxNum); + %ippsCopy_32f(MyBase, MysBase+(*utsi)*(*TypeNum)*(*SpinNum)*SpinMxSliceNum*SpinMxNum, (*TypeNum)*(*SpinNum)*SpinMxSliceNum*SpinMxNum); + %ippsCopy_32f(MxBase, MxsBase+(*utsi)*(*TypeNum)*(*SpinNum)*SpinMxSliceNum*SpinMxNum, (*TypeNum)*(*SpinNum)*SpinMxSliceNum*SpinMxNum); + Muts(utsi+1) = Muts(utsi) + dt; + break; + end + end +catch me + disp(me.message) + +end +end \ No newline at end of file diff --git a/Src/Main/SimuPanel.m b/Src/Main/SimuPanel.m index f73fe24..5c754bc 100644 --- a/Src/Main/SimuPanel.m +++ b/Src/Main/SimuPanel.m @@ -94,6 +94,12 @@ function SimuPanel_OpeningFcn(hObject, eventdata, handles, varargin) %System Hardware Check try handles.CPUInfo=DoCPUOSChk; + if ismac & ~isfield(handles.CPUInfo,'NumThreads') + handles.CPUInfo.NumThreads = 4*handles.CPUInfo.NumProcessors; + end + handles.CPU_MATLAB_uimenu=uimenu(handles.SelectPU_uimenu,'Label','MATLAB Simulation (SLOW!!)'); + set(handles.CPU_MATLAB_uimenu,'Callback',{@Uimenu_ChkFunction, handles.CPU_MATLAB_uimenu}); + set(handles.CPU_MATLAB_uimenu,'Checked','off'); if ~isempty(handles.CPUInfo) handles.CPU_uimenu=uimenu(handles.SelectPU_uimenu,'Label',handles.CPUInfo.Name); set(handles.CPU_uimenu,'Callback',{@Uimenu_ChkFunction, handles.CPU_uimenu}); @@ -894,14 +900,47 @@ function Scan_pushbutton_Callback(hObject, eventdata, handles) end set(hObject,'String','Scan...'); pause(0.01); - DoScanAtCPU; % global (VSeq,VObj,VCtl,VMag,VCoi,VVar,VSig) are needed - else - error('No processing unit is selected.'); + DoScanAtCPU; % global (VSeq,VObj,VCtl,VMag,VCoi,VVar,VSig) are needed + Exe = 1; + end + end + end + + if Exe == 0 + if isfield(handles,'CPU_MATLAB_uimenu') + if strcmp(get(handles.CPU_MATLAB_uimenu,'Checked'),'on') + if isfield(VCtl, 'MT_Flag') + if strcmp(VCtl.MT_Flag, 'on') + error('CPU engine currently doesn''t support Magnetization Transfer model.'); + end + elseif isfield(VCtl, 'ME_Flag') + if strcmp(VCtl.ME_Flag, 'on') + error('CPU engine currently doesn''t support Multiple Pool Exchange model.'); + end + elseif isfield(VCtl, 'CEST_Flag') + if strcmp(VCtl.CEST_Flag, 'on') + error('CPU engine currently doesn''t support Chemical Exchange Saturation Transfer model.'); + end + elseif isfield(VCtl, 'GM_Flag') + if strcmp(VCtl.GM_Flag, 'on') + error('CPU engine currently doesn''t support General Model.'); + end + end + + if handles.BatchFlag==1 + error('DoScanAtCPU_MATLAB'); + end + set(hObject,'String','Scan...'); + pause(0.01); + DoScanAtCPU_MATLAB; % global (VSeq,VObj,VCtl,VMag,VCoi,VVar,VSig) are needed + Exe = 1; end - else - error('No processing unit is selected.'); end end + + if Exe == 0 + error('No processing unit is selected.'); + end set(hObject,'String','Post...'); pause(0.01); DoPostScan(handles);