Commit 39a92d0390bbd9240b2493d5a98afd9f17f0633b

Authored by Pavel Govyadinov
2 parents 8823488b 9b563709

Merge branch 'master' of git.stim.ee.uh.edu:codebase/stimlib into Graph

cmake/FindFANN.cmake 0 → 100644
  1 +#
  2 +# Windows users: define the GLEW_PATH environment variable to point
  3 +# to the directory containing:
  4 +# include/fann.h
  5 +# lib/*fann.lib
  6 +
  7 +
  8 +# FANN_FOUND - system has fann
  9 +# FANN_INCLUDE_DIRS - the fann include directory
  10 +# FANN_LIBRARIES - Link these to use fann
  11 +# FANN_DEFINITIONS - Compiler switches required for using fann
  12 +#
  13 +
  14 +if(FANN_LIBRARIES AND FANN_INCLUDE_DIRS)
  15 + set(FANN_FOUND TRUE)
  16 +else()
  17 + find_path(FANN_INCLUDE_DIR
  18 + NAMES
  19 + fann.h
  20 + PATHS
  21 + $ENV{FANN_PATH}/include
  22 + ${FANN_DIR}/include
  23 + /usr/include
  24 + /usr/local/include
  25 + /opt/local/include
  26 + /sw/include
  27 + )
  28 +
  29 + set( _libraries fann doublefann fixedfann floatfann )
  30 +
  31 + foreach( _lib ${_libraries} )
  32 + string( TOUPPER ${_lib} _name )
  33 +
  34 + find_library(${_name}_LIBRARY
  35 + NAMES
  36 + ${_lib}
  37 + PATHS
  38 + $ENV{FANN_PATH}/lib
  39 + ${FANN_DIR}/lib
  40 + /usr/lib
  41 + /usr/local/lib
  42 + /opt/local/lib
  43 + /sw/lib
  44 + )
  45 +
  46 + endforeach()
  47 +
  48 +
  49 + set(FANN_INCLUDE_DIRS
  50 + ${FANN_INCLUDE_DIR}
  51 + )
  52 +
  53 + set(FANN_LIBRARIES
  54 + ${FANN_LIBRARIES}
  55 + ${FANN_LIBRARY}
  56 + ${DOUBLEFANN_LIBRARY}
  57 + ${FIXEDFANN_LIBRARY}
  58 + ${FLOATFANN_LIBRARY}
  59 + )
  60 +
  61 + if( UNIX )
  62 + set( FANN_LIBRARIES ${FANN_LIBRARIES} m )
  63 + endif()
  64 +
  65 + if(FANN_INCLUDE_DIRS AND FANN_LIBRARIES)
  66 + set(FANN_FOUND TRUE)
  67 + endif()
  68 +
  69 + if(FANN_FOUND)
  70 + if(NOT FANN_FIND_QUIETLY)
  71 + message(STATUS "Found FANN:")
  72 + message(STATUS "FANN_INCLUDE_DIRS: ${FANN_INCLUDE_DIRS}")
  73 + message(STATUS "FANN_LIBRARIES: ${FANN_LIBRARIES}")
  74 + endif()
  75 + else()
  76 + if(FANN_FIND_REQUIRED)
  77 + message(FATAL_ERROR "Could not find FANN")
  78 + endif()
  79 + endif()
  80 +
  81 + mark_as_advanced(FANN_INCLUDE_DIRS FANN_LIBRARIES)
  82 +endif()
... ...
cmake/FindGLEW.cmake 0 → 100644
  1 +# Copyright (c) 2012-2016 DreamWorks Animation LLC
  2 +#
  3 +# All rights reserved. This software is distributed under the
  4 +# Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
  5 +#
  6 +# Redistributions of source code must retain the above copyright
  7 +# and license notice and the following restrictions and disclaimer.
  8 +#
  9 +# * Neither the name of DreamWorks Animation nor the names of
  10 +# its contributors may be used to endorse or promote products derived
  11 +# from this software without specific prior written permission.
  12 +#
  13 +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  14 +# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  15 +# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
  16 +# A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
  17 +# OWNER OR CONTRIBUTORS BE LIABLE FOR ANY INDIRECT, INCIDENTAL,
  18 +# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
  19 +# LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  20 +# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
  21 +# THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  22 +# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
  23 +# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  24 +# IN NO EVENT SHALL THE COPYRIGHT HOLDERS' AND CONTRIBUTORS' AGGREGATE
  25 +# LIABILITY FOR ALL CLAIMS REGARDLESS OF THEIR BASIS EXCEED US$250.00.
  26 +#
  27 +
  28 +#-*-cmake-*-
  29 +# - Find GLEW
  30 +#
  31 +# Author : Nicholas Yue yue.nicholas@gmail.com
  32 +#
  33 +# This auxiliary CMake file helps in find the GLEW headers and libraries
  34 +#
  35 +# GLEW_FOUND set if Glew is found.
  36 +# GLEW_INCLUDE_DIR GLEW's include directory
  37 +# GLEW_glew_LIBRARY GLEW libraries
  38 +# GLEW_glewmx_LIBRARY GLEWmx libraries (Mulitple Rendering Context)
  39 +
  40 +FIND_PACKAGE ( PackageHandleStandardArgs )
  41 +
  42 +FIND_PATH( GLEW_LOCATION include/GL/glew.h
  43 + "$ENV{GLEW_ROOT}"
  44 + NO_DEFAULT_PATH
  45 + NO_SYSTEM_ENVIRONMENT_PATH
  46 + )
  47 +
  48 +FIND_PACKAGE_HANDLE_STANDARD_ARGS ( GLEW
  49 + REQUIRED_VARS GLEW_LOCATION
  50 + )
  51 +
  52 +IF ( GLEW_LOCATION )
  53 +
  54 + SET( GLEW_INCLUDE_DIR "${GLEW_LOCATION}/include" CACHE STRING "GLEW include path")
  55 +
  56 + SET ( ORIGINAL_CMAKE_FIND_LIBRARY_SUFFIXES ${CMAKE_FIND_LIBRARY_SUFFIXES})
  57 + IF (GLEW_USE_STATIC_LIBS)
  58 + IF (APPLE)
  59 + SET(CMAKE_FIND_LIBRARY_SUFFIXES ".a")
  60 + FIND_LIBRARY ( GLEW_LIBRARY_PATH GLEW PATHS ${GLEW_LOCATION}/lib
  61 + NO_DEFAULT_PATH
  62 + NO_SYSTEM_ENVIRONMENT_PATH
  63 + )
  64 + FIND_LIBRARY ( GLEWmx_LIBRARY_PATH GLEWmx PATHS ${GLEW_LOCATION}/lib
  65 + NO_DEFAULT_PATH
  66 + NO_SYSTEM_ENVIRONMENT_PATH
  67 + )
  68 + # MESSAGE ( "APPLE STATIC" )
  69 + # MESSAGE ( "GLEW_LIBRARY_PATH = " ${GLEW_LIBRARY_PATH} )
  70 + ELSEIF (WIN32)
  71 + # Link library
  72 + SET(CMAKE_FIND_LIBRARY_SUFFIXES ".lib")
  73 + FIND_LIBRARY ( GLEW_LIBRARY_PATH GLEW32S PATHS ${GLEW_LOCATION}/lib )
  74 + FIND_LIBRARY ( GLEWmx_LIBRARY_PATH GLEW32MXS PATHS ${GLEW_LOCATION}/lib )
  75 + ELSE (APPLE)
  76 + SET(CMAKE_FIND_LIBRARY_SUFFIXES ".a")
  77 + FIND_LIBRARY ( GLEW_LIBRARY_PATH GLEW PATHS ${GLEW_LOCATION}/lib
  78 + NO_DEFAULT_PATH
  79 + NO_SYSTEM_ENVIRONMENT_PATH
  80 + )
  81 + FIND_LIBRARY ( GLEWmx_LIBRARY_PATH GLEWmx PATHS ${GLEW_LOCATION}/lib
  82 + NO_DEFAULT_PATH
  83 + NO_SYSTEM_ENVIRONMENT_PATH
  84 + )
  85 + # MESSAGE ( "LINUX STATIC" )
  86 + # MESSAGE ( "GLEW_LIBRARY_PATH = " ${GLEW_LIBRARY_PATH} )
  87 + ENDIF (APPLE)
  88 + ELSE ()
  89 + IF (APPLE)
  90 + SET(CMAKE_FIND_LIBRARY_SUFFIXES ".dylib")
  91 + FIND_LIBRARY ( GLEW_LIBRARY_PATH GLEW PATHS ${GLEW_LOCATION}/lib )
  92 + FIND_LIBRARY ( GLEWmx_LIBRARY_PATH GLEWmx PATHS ${GLEW_LOCATION}/lib )
  93 + ELSEIF (WIN32)
  94 + # Link library
  95 + SET(CMAKE_FIND_LIBRARY_SUFFIXES ".lib")
  96 + FIND_LIBRARY ( GLEW_LIBRARY_PATH GLEW32 PATHS ${GLEW_LOCATION}/lib )
  97 + FIND_LIBRARY ( GLEWmx_LIBRARY_PATH GLEW32mx PATHS ${GLEW_LOCATION}/lib )
  98 + # Load library
  99 + SET(CMAKE_FIND_LIBRARY_SUFFIXES ".dll")
  100 + FIND_LIBRARY ( GLEW_DLL_PATH GLEW32 PATHS ${GLEW_LOCATION}/bin
  101 + NO_DEFAULT_PATH
  102 + NO_SYSTEM_ENVIRONMENT_PATH
  103 + )
  104 + FIND_LIBRARY ( GLEWmx_DLL_PATH GLEW32mx PATHS ${GLEW_LOCATION}/bin
  105 + NO_DEFAULT_PATH
  106 + NO_SYSTEM_ENVIRONMENT_PATH
  107 + )
  108 + ELSE (APPLE)
  109 + # Unices
  110 + FIND_LIBRARY ( GLEW_LIBRARY_PATH GLEW PATHS ${GLEW_LOCATION}/lib
  111 + NO_DEFAULT_PATH
  112 + NO_SYSTEM_ENVIRONMENT_PATH
  113 + )
  114 + FIND_LIBRARY ( GLEWmx_LIBRARY_PATH GLEWmx PATHS ${GLEW_LOCATION}/lib
  115 + NO_DEFAULT_PATH
  116 + NO_SYSTEM_ENVIRONMENT_PATH
  117 + )
  118 + ENDIF (APPLE)
  119 + ENDIF ()
  120 + # MUST reset
  121 + SET(CMAKE_FIND_LIBRARY_SUFFIXES ${ORIGINAL_CMAKE_FIND_LIBRARY_SUFFIXES})
  122 +
  123 + SET( GLEW_GLEW_LIBRARY ${GLEW_LIBRARY_PATH} CACHE STRING "GLEW library")
  124 + SET( GLEW_GLEWmx_LIBRARY ${GLEWmx_LIBRARY_PATH} CACHE STRING "GLEWmx library")
  125 +
  126 +ENDIF ()
... ...
cmake/FindGLUT.cmake 0 → 100644
  1 +#.rst:
  2 +# FindGLUT
  3 +# --------
  4 +#
  5 +# try to find glut library and include files.
  6 +#
  7 +# IMPORTED Targets
  8 +# ^^^^^^^^^^^^^^^^
  9 +#
  10 +# This module defines the :prop_tgt:`IMPORTED` targets:
  11 +#
  12 +# ``GLUT::GLUT``
  13 +# Defined if the system has GLUT.
  14 +#
  15 +# Result Variables
  16 +# ^^^^^^^^^^^^^^^^
  17 +#
  18 +# This module sets the following variables:
  19 +#
  20 +# ::
  21 +#
  22 +# GLUT_INCLUDE_DIR, where to find GL/glut.h, etc.
  23 +# GLUT_LIBRARIES, the libraries to link against
  24 +# GLUT_FOUND, If false, do not try to use GLUT.
  25 +#
  26 +# Also defined, but not for general use are:
  27 +#
  28 +# ::
  29 +#
  30 +# GLUT_glut_LIBRARY = the full path to the glut library.
  31 +# GLUT_Xmu_LIBRARY = the full path to the Xmu library.
  32 +# GLUT_Xi_LIBRARY = the full path to the Xi Library.
  33 +
  34 +#=============================================================================
  35 +# Copyright 2001-2009 Kitware, Inc.
  36 +#
  37 +# Distributed under the OSI-approved BSD License (the "License");
  38 +# see accompanying file Copyright.txt for details.
  39 +#
  40 +# This software is distributed WITHOUT ANY WARRANTY; without even the
  41 +# implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  42 +# See the License for more information.
  43 +#=============================================================================
  44 +# (To distribute this file outside of CMake, substitute the full
  45 +# License text for the above reference.)
  46 +
  47 +if (WIN32)
  48 + find_path( GLUT_INCLUDE_DIR NAMES GL/glut.h
  49 + PATHS $ENV{GLUT_ROOT_PATH}/include )
  50 +
  51 + if( CMAKE_SIZEOF_VOID_P EQUAL 8 )
  52 + find_library( GLUT_glut_LIBRARY NAMES freeglut
  53 + PATHS
  54 + $ENV{GLUT_ROOT_PATH}/lib/x64
  55 +
  56 + NO_DEFAULT_PATH
  57 + )
  58 + else( CMAKE_SIZEOF_VOID_P EQUAL 8 )
  59 + find_library( GLUT_glut_LIBRARY NAMES glut glut32 freeglut
  60 + PATHS
  61 + ${OPENGL_LIBRARY_DIR}
  62 + $ENV{GLUT_ROOT_PATH}/lib
  63 + )
  64 + endif( CMAKE_SIZEOF_VOID_P EQUAL 8 )
  65 +
  66 +else ()
  67 +
  68 + if (APPLE)
  69 + find_path(GLUT_INCLUDE_DIR glut.h ${OPENGL_LIBRARY_DIR})
  70 + find_library(GLUT_glut_LIBRARY GLUT DOC "GLUT library for OSX")
  71 + find_library(GLUT_cocoa_LIBRARY Cocoa DOC "Cocoa framework for OSX")
  72 +
  73 + if(GLUT_cocoa_LIBRARY AND NOT TARGET GLUT::Cocoa)
  74 + add_library(GLUT::Cocoa UNKNOWN IMPORTED)
  75 + # Cocoa should always be a Framework, but we check to make sure.
  76 + if(GLUT_cocoa_LIBRARY MATCHES "/([^/]+)\\.framework$")
  77 + set_target_properties(GLUT::Cocoa PROPERTIES
  78 + IMPORTED_LOCATION "${GLUT_cocoa_LIBRARY}/${CMAKE_MATCH_1}")
  79 + else()
  80 + set_target_properties(GLUT::Cocoa PROPERTIES
  81 + IMPORTED_LOCATION "${GLUT_cocoa_LIBRARY}")
  82 + endif()
  83 + endif()
  84 + else ()
  85 +
  86 + if (BEOS)
  87 +
  88 + set(_GLUT_INC_DIR /boot/develop/headers/os/opengl)
  89 + set(_GLUT_glut_LIB_DIR /boot/develop/lib/x86)
  90 +
  91 + else()
  92 +
  93 + find_library( GLUT_Xi_LIBRARY Xi
  94 + /usr/openwin/lib
  95 + )
  96 +
  97 + find_library( GLUT_Xmu_LIBRARY Xmu
  98 + /usr/openwin/lib
  99 + )
  100 +
  101 + if(GLUT_Xi_LIBRARY AND NOT TARGET GLUT::Xi)
  102 + add_library(GLUT::Xi UNKNOWN IMPORTED)
  103 + set_target_properties(GLUT::Xi PROPERTIES
  104 + IMPORTED_LOCATION "${GLUT_Xi_LIBRARY}")
  105 + endif()
  106 +
  107 + if(GLUT_Xmu_LIBRARY AND NOT TARGET GLUT::Xmu)
  108 + add_library(GLUT::Xmu UNKNOWN IMPORTED)
  109 + set_target_properties(GLUT::Xmu PROPERTIES
  110 + IMPORTED_LOCATION "${GLUT_Xmu_LIBRARY}")
  111 + endif()
  112 +
  113 + endif ()
  114 +
  115 + find_path( GLUT_INCLUDE_DIR GL/glut.h
  116 + /usr/include/GL
  117 + /usr/openwin/share/include
  118 + /usr/openwin/include
  119 + /opt/graphics/OpenGL/include
  120 + /opt/graphics/OpenGL/contrib/libglut
  121 + ${_GLUT_INC_DIR}
  122 + )
  123 +
  124 + find_library( GLUT_glut_LIBRARY glut
  125 + /usr/openwin/lib
  126 + ${_GLUT_glut_LIB_DIR}
  127 + )
  128 +
  129 + unset(_GLUT_INC_DIR)
  130 + unset(_GLUT_glut_LIB_DIR)
  131 +
  132 + endif ()
  133 +
  134 +endif ()
  135 +
  136 +FIND_PACKAGE_HANDLE_STANDARD_ARGS(GLUT REQUIRED_VARS GLUT_glut_LIBRARY GLUT_INCLUDE_DIR)
  137 +
  138 +if (GLUT_FOUND)
  139 + # Is -lXi and -lXmu required on all platforms that have it?
  140 + # If not, we need some way to figure out what platform we are on.
  141 + set( GLUT_LIBRARIES
  142 + ${GLUT_glut_LIBRARY}
  143 + ${GLUT_Xmu_LIBRARY}
  144 + ${GLUT_Xi_LIBRARY}
  145 + ${GLUT_cocoa_LIBRARY}
  146 + )
  147 +
  148 + if(NOT TARGET GLUT::GLUT)
  149 + add_library(GLUT::GLUT UNKNOWN IMPORTED)
  150 + set_target_properties(GLUT::GLUT PROPERTIES
  151 + INTERFACE_INCLUDE_DIRECTORIES "${GLUT_INCLUDE_DIR}")
  152 + if(GLUT_glut_LIBRARY MATCHES "/([^/]+)\\.framework$")
  153 + set_target_properties(GLUT::GLUT PROPERTIES
  154 + IMPORTED_LOCATION "${GLUT_glut_LIBRARY}/${CMAKE_MATCH_1}")
  155 + else()
  156 + set_target_properties(GLUT::GLUT PROPERTIES
  157 + IMPORTED_LOCATION "${GLUT_glut_LIBRARY}")
  158 + endif()
  159 +
  160 + if(TARGET GLUT::Xmu)
  161 + set_property(TARGET GLUT::GLUT APPEND
  162 + PROPERTY INTERFACE_LINK_LIBRARIES GLUT::Xmu)
  163 + endif()
  164 +
  165 + if(TARGET GLUT::Xi)
  166 + set_property(TARGET GLUT::GLUT APPEND
  167 + PROPERTY INTERFACE_LINK_LIBRARIES GLUT::Xi)
  168 + endif()
  169 +
  170 + if(TARGET GLUT::Cocoa)
  171 + set_property(TARGET GLUT::GLUT APPEND
  172 + PROPERTY INTERFACE_LINK_LIBRARIES GLUT::Cocoa)
  173 + endif()
  174 + endif()
  175 +
  176 + #The following deprecated settings are for backwards compatibility with CMake1.4
  177 + set (GLUT_LIBRARY ${GLUT_LIBRARIES})
  178 + set (GLUT_INCLUDE_PATH ${GLUT_INCLUDE_DIR})
  179 +endif()
  180 +
  181 +mark_as_advanced(
  182 + GLUT_INCLUDE_DIR
  183 + GLUT_glut_LIBRARY
  184 + GLUT_Xmu_LIBRARY
  185 + GLUT_Xi_LIBRARY
  186 + )
... ...
cmake/FindSTIM.cmake
1   -include(FindPackageHandleStandardArgs)
2   -
3   -set(STIM_INCLUDE_DIR $ENV{STIMLIB_PATH})
4   -
5   -find_package_handle_standard_args(STIM DEFAULT_MSG STIM_INCLUDE_DIR)
6   -
7   -if(STIM_FOUND)
8   - set(STIM_INCLUDE_DIRS ${STIM_INCLUDE_DIR})
9   -endif()
10 1 \ No newline at end of file
  2 +# finds the STIM library (downloads it if it isn't present)
  3 +# set STIMLIB_PATH to the directory containing the stim subdirectory (the stim repository)
  4 +
  5 +include(FindPackageHandleStandardArgs)
  6 +
  7 +set(STIM_INCLUDE_DIR $ENV{STIMLIB_PATH})
  8 +
  9 +find_package_handle_standard_args(STIM DEFAULT_MSG STIM_INCLUDE_DIR)
  10 +
  11 +if(STIM_FOUND)
  12 + set(STIM_INCLUDE_DIRS ${STIM_INCLUDE_DIR})
  13 +elseif(STIM_FOUND)
  14 + #if the STIM library isn't found, download it
  15 + #file(REMOVE_RECURSE ${CMAKE_BINARY_DIR}/stimlib) #remove the stimlib directory if it exists
  16 + #set(STIM_GIT "https://git.stim.ee.uh.edu/codebase/stimlib.git")
  17 + #execute_process(COMMAND git clone --depth 1 ${STIM_GIT} WORKING_DIRECTORY ${CMAKE_BINARY_DIR})
  18 + #set(STIM_INCLUDE_DIRS "${CMAKE_BINARY_DIR}/stimlib" CACHE TYPE PATH)
  19 + message("STIM library not found. Set the STIMLIB_PATH environment variable to the STIMLIB location.")
  20 + message("STIMLIB can be found here: https://git.stim.ee.uh.edu/codebase/stimlib")
  21 +endif(STIM_FOUND)
  22 +
  23 +find_package_handle_standard_args(STIM DEFAULT_MSG STIM_INCLUDE_DIR)
... ...
matlab/bsq2tensorflow.m 0 → 100644
  1 +function T = bsq2tensorflow(I, n)
  2 +
  3 + sx = size(I, 1);
  4 + sy = size(I, 2) / n; %get the size of the tensor along Y
  5 + sb = size(I, 3);
  6 +
  7 + T = zeros(sx * sy * sb, n); %allocate space for the output matrix
  8 + for i = 0:n-1
  9 + ti = I(:, i * sy + 1 : i * sy + sy, :);
  10 + T(:, i+1) = ti(:);
  11 + end
  12 +end
  13 +
  14 +
0 15 \ No newline at end of file
... ...
matlab/enviLoadRaw.m
1 1 %loads an ENVI file without any manipulation (changing orientation)
  2 +% enviLoadRaw(filename, headername)
2 3 function M = enviLoadRaw(filename, headername)
3 4  
4 5 %if a header isn't provided, assume it's just the filename
... ...
matlab/enviSaveRaw.m
1 1 %saves an ENVI file without any manipulation, assumes (X, Y, S)
  2 +% enviSaveRaw(M, filename, headername)
2 3 function enviSaveRaw(M, filename, headername)
3 4  
4 5 %if a header isn't provided, assume it's just the filename
... ...
matlab/readspe.m 0 → 100644
  1 +% Read images of TIFF, SPE2.2(WinSpec) and SPE3.0(Lightfield)
  2 +% Version: JTL Jun-9-2016
  3 +% ----------------- READ THIS FIRST !!!!! --------------------------------
  4 +% Change the file name to "readspe" before use
  5 +% Example:
  6 +% Z = readspe(filename)
  7 +% Z = readspe(filename,'info')
  8 +% Z = readspe(filename,frame_index)
  9 +% Z = readspe(filename,frame_index,'info')
  10 +% Input:
  11 +% filename - filename string, e.g. 'image.spe'
  12 +% frame_index - frame index, start from 1
  13 +% If you have multiple frames, use a "for" loop
  14 +% 'info' - flag to show file info, i.e. dimension, number of frames, version
  15 +% Output:
  16 +% Z - UINT16 image (convert to double if you need)
  17 +% ------------------------------------------------------------------------
  18 +% Z = readspe (filename,frame_index,'info')
  19 +function Z = readspe (filename,varargin)
  20 +
  21 +if exist(filename) == 2
  22 +
  23 + Nfr = 1; % default read first frame
  24 + if nargin >1
  25 + if isa(varargin{1},'numeric')
  26 + Nfr = varargin{1};
  27 + end
  28 + end
  29 +
  30 + [~,name,ext] = fileparts(filename);
  31 + switch upper(ext)
  32 + case '.TIFF'
  33 + file_ver = 'TIFF';
  34 + Z = imread(filename);
  35 + [Y,X] = size(Z);
  36 + % datatype = class(Z)
  37 +
  38 + case '.SPE'
  39 + fid = fopen(filename);
  40 + I = fread(fid,Inf,'uint8');
  41 + X = double(typecast(uint8(I(43:44)),'uint16'));
  42 + Y = double(typecast(uint8(I(657:658)),'uint16'));
  43 + fr = typecast(uint8(I(1447:1450)),'int32');
  44 + spe_ver = typecast(uint8(I(1993:1996)),'single');
  45 + file_ver = ['SPE ' num2str(spe_ver)];
  46 + datatypeN = typecast(uint8(I(109:110)),'int16');
  47 + switch datatypeN
  48 + case 0 % 32-bit float
  49 + datatype = 'single'; datalength = 4;
  50 + case 1 % 32-bit signed integer
  51 + datatype = 'int32'; datalength = 4;
  52 + case 2 % 16-bit signed integer
  53 + datatype = 'int16'; datalength = 2;
  54 + case 3 % 16-bit unsigned integer
  55 + datatype = 'uint16'; datalength = 2;
  56 + case 8 % 32-bit unsigned integer
  57 + datatype = 'uint32'; datalength = 4;
  58 + end
  59 + % A = I(4101:4100+X*Y*2); % Default read first frame
  60 + A = I(4101+X*Y*datalength*(Nfr-1):4100+X*Y*datalength*Nfr);
  61 + B = typecast(uint8(A),datatype); % important
  62 + Z = reshape(B,X,Y);
  63 + Z = Z';
  64 + fclose(fid);
  65 + end
  66 +
  67 + if nargin >1
  68 + if varargin{end} == 'info'
  69 + display(['X = ' num2str(X)]);
  70 + display(['Y = ' num2str(Y)]);
  71 + if(exist('fr','var'));display(['Number of Frames: ' num2str(fr)]);end;
  72 + display(['File version: ' file_ver]);
  73 + end
  74 + end
  75 +
  76 +elseif exist(filename) == 0
  77 + display('File does not exist!');
  78 +end
0 79 \ No newline at end of file
... ...
matlab/spe2envi.m 0 → 100644
  1 +function spe2envi(filemask, outfile)
  2 +
  3 + filelist = dir(filemask);
  4 +
  5 + %get a list of date numbers
  6 + datenums = cell2mat({filelist.datenum});
  7 +
  8 + %sort the file order based on acquisition time
  9 + [~, id] = sort(datenums);
  10 +
  11 + %get the number of files
  12 + Y = length(id); %size of the image along Y
  13 +
  14 + %load the first file to determine the spectral and X-axis size
  15 + temp = readspe(filelist(1).name);
  16 + X = size(temp, 1); %size of the image along X
  17 + B = size(temp, 2); %number of bands in the image
  18 +
  19 + %create the cube
  20 + I = zeros(X, Y, B);
  21 +
  22 + %for each line
  23 + for y = 1:Y
  24 +
  25 + %read a SPE file
  26 + img = readspe(filelist(id(y)).name);
  27 +
  28 + I(:, y, :) = permute(img, [1 3 2]);
  29 + end
  30 +
  31 + enviSaveRaw(single(I), outfile, [outfile '.hdr']);
  32 +
  33 +
  34 +
... ...
matlab/brewermap.m renamed to matlab/stimBrewerMap.m
matlab/stimLoadAgilent.m 0 → 100644
  1 +%Loads a standard Agilent ResPro binary file
  2 +% stimLoadAgilent(filename)
  3 +function S = stimLoadAgilent(filename)
  4 +
  5 + fid = fopen(filename);
  6 + fseek(fid, 9, 'bof');
  7 + Z = fread(fid, 1, 'uint16');
  8 + fseek(fid, 13, 'cof');
  9 + X = fread(fid, 1, 'uint16');
  10 + Y = fread(fid, 1, 'uint16');
  11 +
  12 + fseek(fid, 1020, 'bof');
  13 +
  14 + S = reshape(fread(fid, [X, Y * Z], 'float32'), [X, Y, Z]);
  15 +
  16 +
0 17 \ No newline at end of file
... ...
matlab/stimROC.m 0 → 100644
  1 +function [TPR, FPR, AUC] = stimROC(C, T)
  2 +%build an ROC curve
  3 +% C - class labels as an array of binary values (1 = true positive)
  4 +% T - threshold used for classification
  5 +
  6 + %sort the thresholds in descending order and get the indices
  7 + [~, I] = sort(T, 'descend');
  8 +
  9 + %sort the class labels in the same order as the thresholds
  10 + Cs = C(I);
  11 +
  12 + %calculate the number of measurements
  13 + M = size(C, 2);
  14 +
  15 + %calculate the number of positives
  16 + P = nnz(C);
  17 +
  18 + %calculate the number of negatives
  19 + N = M - P;
  20 +
  21 + %if all examples are positives or negatives, return a perfect score?
  22 + if P == M
  23 + error('ERROR: no positive observations');
  24 + end
  25 + if P == 0
  26 + error('ERROR: no negative observations');
  27 + end
  28 +
  29 + %allocate space for the ROC curve
  30 + TPR = zeros(1, M);
  31 + FPR = zeros(1, M);
  32 +
  33 +
  34 +
  35 + %calculate the number of inflection points
  36 + ip = 0;
  37 + for i = 2:M
  38 + if Cs(i) ~= Cs(i-1)
  39 + ip = ip + 1;
  40 + end
  41 + end
  42 +
  43 + %initialize the true and false positive rates to zero
  44 + TP = 0;
  45 + FP = 0;
  46 + for i = 1:M
  47 + if Cs(i) == 1
  48 + TP = TP + 1;
  49 + else
  50 + FP = FP + 1;
  51 + end
  52 +
  53 + TPR(i) = TP / P;
  54 + FPR(i) = FP / N;
  55 + end
  56 +
  57 + %calculate the area under the ROC curve
  58 + AUC = 0;
  59 + for i = 2:M
  60 + w = FPR(i) - FPR(i-1);
  61 + h = TPR(i);
  62 + AUC = AUC + w * h;
  63 + end
  64 +
  65 +
  66 +
  67 +
  68 +
  69 +
0 70 \ No newline at end of file
... ...
python/enviProcess.py 0 → 100644
  1 +#!/usr/bin/python3
  2 +
  3 +#import system processes
  4 +import subprocess, sys
  5 +
  6 +if len(sys.argv) > 1:
  7 + infile = int(sys.argv[1])
  8 +
  9 +basefile = infile + "-base"
  10 +normfile = infile + "-norm"
  11 +
  12 +runcommand = "hsiproc " + infile + basefile + " --baseline baseline.txt"
  13 +subprocess.call(runcommand, shell=True)
0 14 \ No newline at end of file
... ...
stim/biomodels/cellset.h
... ... @@ -117,7 +117,7 @@ public:
117 117 }
118 118  
119 119 /// Return the maximum value of a field in this cell set
120   - double max(std::string field){
  120 + double maximum(std::string field){
121 121 size_t idx = fields[field]; //get the field index
122 122 size_t ncells = cells.size(); //get the total number of cells
123 123 double maxval, val; //stores the current and maximum values
... ... @@ -130,7 +130,7 @@ public:
130 130 }
131 131  
132 132 /// Return the maximum value of a field in this cell set
133   - double min(std::string field){
  133 + double minimum(std::string field){
134 134 size_t idx = fields[field]; //get the field index
135 135 size_t ncells = cells.size(); //get the total number of cells
136 136 double minval, val; //stores the current and maximum values
... ...
stim/biomodels/network.h
... ... @@ -11,8 +11,8 @@
11 11 #include <stim/math/vec3.h>
12 12 #include <stim/visualization/obj.h>
13 13 #include <stim/visualization/cylinder.h>
14   -#include <ANN/ANN.h>
15   -#include <boost/tuple/tuple.hpp>
  14 +#include <stim/structures/kdtree.cuh>
  15 +#include <stim/cuda/cudatools/timer.h>
16 16  
17 17  
18 18 namespace stim{
... ... @@ -35,7 +35,7 @@ class network{
35 35 // default constructor
36 36 edge() : cylinder<T>()
37 37 {
38   - v[1] = -1; v[0] = -1;
  38 + v[1] = (unsigned)(-1); v[0] = (unsigned)(-1);
39 39 }
40 40 /// Constructor - creates an edge from a list of points by calling the stim::fiber constructor
41 41  
... ... @@ -57,7 +57,7 @@ class network{
57 57 /// Output the edge information as a string
58 58 std::string str(){
59 59 std::stringstream ss;
60   - ss<<"("<<cylinder<T>::size()<<")\tl = "<<this.length()<<"\t"<<v[0]<<"----"<<v[1];
  60 + ss<<"("<<cylinder<T>::size()<<")\tl = "<<this->length()<<"\t"<<v[0]<<"----"<<v[1];
61 61 return ss.str();
62 62 }
63 63  
... ... @@ -125,7 +125,9 @@ public:
125 125 return V.size();
126 126 }
127 127  
128   - std::vector<vertex> operator*(T s){
  128 + //scale the network by some constant value
  129 + // I don't think these work??????
  130 + /*std::vector<vertex> operator*(T s){
129 131 for (unsigned i=0; i< vertices; i ++ ){
130 132 V[i] = V[i] * s;
131 133 }
... ... @@ -139,10 +141,9 @@ public:
139 141 }
140 142 }
141 143 return V;
142   - }
  144 + }*/
143 145  
144 146 // Returns an average of branching index in the network
145   -
146 147 double BranchingIndex(){
147 148 double B=0;
148 149 for(unsigned v=0; v < V.size(); v ++){
... ... @@ -154,7 +155,6 @@ public:
154 155 }
155 156  
156 157 // Returns number of branch points in thenetwork
157   -
158 158 unsigned int BranchP(){
159 159 unsigned int B=0;
160 160 unsigned int c;
... ... @@ -168,7 +168,6 @@ public:
168 168 }
169 169  
170 170 // Returns number of end points (tips) in thenetwork
171   -
172 171 unsigned int EndP(){
173 172 unsigned int B=0;
174 173 unsigned int c;
... ... @@ -202,10 +201,11 @@ public:
202 201 // return s;
203 202 //}
204 203  
205   -
  204 + //Calculate Metrics---------------------------------------------------
206 205 // Returns an average of fiber/edge lengths in the network
207 206 double Lengths(){
208   - stim::vec<T> L;double sumLength = 0;
  207 + stim::vec<T> L;
  208 + double sumLength = 0;
209 209 for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network
210 210 L.push_back(E[e].length()); //append the edge length
211 211 sumLength = sumLength + E[e].length();
... ... @@ -269,8 +269,10 @@ public:
269 269 double avg = sumFractDim / E.size();
270 270 return avg;
271 271 }
272   - stim::cylinder<T> get_cylinder(unsigned f){
273   - return E[f]; //return the specified edge (casting it to a fiber)
  272 +
  273 + //returns a cylinder represented a given fiber (based on edge index)
  274 + stim::cylinder<T> get_cylinder(unsigned e){
  275 + return E[e]; //return the specified edge (casting it to a fiber)
274 276 }
275 277  
276 278 //load a network from an OBJ file
... ... @@ -385,11 +387,27 @@ public:
385 387 return n;
386 388 }
387 389  
  390 + //Copy the point cloud representing the centerline for the network into an array
  391 + void centerline_cloud(T* dst) {
  392 + size_t p; //stores the current edge point
  393 + size_t P; //stores the number of points in an edge
  394 + size_t i = 0; //index into the output array of points
  395 + for (size_t e = 0; e < E.size(); e++) { //for each edge in the network
  396 + P = E[e].size(); //get the number of points in this edge
  397 + for (p = 0; p < P; p++) {
  398 + dst[i * 3 + 0] = E[e][p][0];
  399 + dst[i * 3 + 1] = E[e][p][1];
  400 + dst[i * 3 + 2] = E[e][p][2];
  401 + i++;
  402 + }
  403 + }
  404 + }
  405 +
388 406 // gaussian function
389 407 float gaussianFunction(float x, float std=25){ return exp(-x/(2*std*std));} // by default std = 25
390 408  
391   - // stim 3d vector to annpoint of 3 dimensions
392   - void stim2ann(ANNpoint &a, stim::vec3<T> b){
  409 + // convert vec3 to array
  410 + void stim2array(float *a, stim::vec3<T> b){
393 411 a[0] = b[0];
394 412 a[1] = b[1];
395 413 a[2] = b[2];
... ... @@ -413,57 +431,81 @@ public:
413 431  
414 432 /// @param A is the network to compare to - the field is generated for A
415 433 /// @param sigma is the user-defined tolerance value - smaller values provide a stricter comparison
416   - stim::network<T> compare(stim::network<T> A, float sigma){
  434 + stim::network<T> compare(stim::network<T> A, float sigma, int device){
417 435  
418   - stim::network<T> R; //generate a network storing the result of the comparison
419   - R = (*this); //initialize the result with the current network
  436 + stim::network<T> R; //generate a network storing the result of the comparison
  437 + R = (*this); //initialize the result with the current network
420 438  
421   - //generate a KD-tree for network A
422   - float metric = 0.0; // initialize metric to be returned after comparing the networks
423   - ANNkd_tree* kdt; // initialize a pointer to a kd tree
424   - double **c; // centerline (array of double pointers) - points on kdtree must be double
425   - unsigned int n_data = A.total_points(); // set the number of points
426   - c = (double**) malloc(sizeof(double*) * n_data); // allocate the array pointer
427   - for(unsigned int i = 0; i < n_data; i++) // allocate space for each point of 3 dimensions
428   - c[i] = (double*) malloc(sizeof(double) * 3);
  439 + T *c; // centerline (array of double pointers) - points on kdtree must be double
  440 + size_t n_data = A.total_points(); // set the number of points
  441 + c = (T*) malloc(sizeof(T) * n_data * 3); //allocate an array to store all points in the data set
429 442  
430 443 unsigned t = 0;
431   - for(unsigned e = 0; e < A.E.size(); e++){ //for each edge in the network
432   - for(unsigned p = 0; p < A.E[e].size(); p++){ //for each point in the edge
  444 + for(unsigned e = 0; e < A.E.size(); e++){ //for each edge in the network
  445 + for(unsigned p = 0; p < A.E[e].size(); p++){ //for each point in the edge
433 446 for(unsigned d = 0; d < 3; d++){ //for each coordinate
434 447  
435   - c[t][d] = A.E[e][p][d];
  448 + c[t * 3 + d] = A.E[e][p][d]; //copy the point into the array c
436 449 }
437 450 t++;
438 451 }
439 452 }
440 453  
  454 + //generate a KD-tree for network A
  455 + //float metric = 0.0; // initialize metric to be returned after comparing the network
  456 + size_t MaxTreeLevels = 3; // max tree level
  457 +
  458 +#ifdef __CUDACC__
  459 + cudaSetDevice(device);
  460 + stim::cuda_kdtree<T, 3> kdt; // initialize a pointer to a kd tree
  461 +
441 462 //compare each point in the current network to the field produced by A
442   - ANNpointArray pts = (ANNpointArray)c; // create an array of data points of type double
443   - kdt = new ANNkd_tree(pts, n_data, 3); // build a KD tree using the annpointarray
444   - double eps = 0; // error bound
445   - ANNdistArray dists = new ANNdist[1]; // near neighbor distances
446   - ANNidxArray nnIdx = new ANNidx[1]; // near neighbor indices // allocate near neigh indices
  463 + kdt.create(c, n_data, MaxTreeLevels); // build a KD tree
  464 + T *dists = new T[1]; // near neighbor distances
  465 + size_t *nnIdx = new size_t[1]; // near neighbor indices // allocate near neigh indices
447 466  
448 467 stim::vec3<T> p0, p1;
449   - float m1;
450   - float M = 0; //stores the total metric value
451   - float L = 0; //stores the total network length
452   - ANNpoint queryPt = annAllocPt(3);
  468 + T m1;
  469 + //float M = 0; //stores the total metric value
  470 + //float L = 0; //stores the total network length
  471 + T* queryPt = new T[3];
453 472 for(unsigned e = 0; e < R.E.size(); e++){ //for each edge in A
454 473 R.E[e].add_mag(0); //add a new magnitude for the metric
455 474  
456 475 for(unsigned p = 0; p < R.E[e].size(); p++){ //for each point in the edge
457 476  
458 477 p1 = R.E[e][p]; //get the next point in the edge
459   - stim2ann(queryPt, p1);
460   - kdt->annkSearch( queryPt, 1, nnIdx, dists, eps); //find the distance between A and the current network
461   - m1 = 1.0f - gaussianFunction((float)dists[0], sigma); //calculate the metric value based on the distance
  478 + stim2array(queryPt, p1);
  479 + kdt.search(queryPt, 1, nnIdx, dists); //find the distance between A and the current network
  480 +
  481 + m1 = 1.0f - gaussianFunction((T)dists[0], sigma); //calculate the metric value based on the distance
462 482 R.E[e].set_mag(m1, p, 1); //set the error for the second point in the segment
463 483  
464 484 }
465 485 }
  486 +#else
  487 + stim::cpu_kdtree<T, 3> kdt;
  488 + kdt.create(c, n_data, MaxTreeLevels);
  489 + T *dists = new T[1]; // near neighbor distances
  490 + size_t *nnIdx = new size_t[1]; // near neighbor indices // allocate near neigh indices
  491 +
  492 + stim::vec3<T> p0, p1;
  493 + T m1;
  494 + T* queryPt = new T[3];
  495 + for(unsigned e = 0; e < R.E.size(); e++){ //for each edge in A
  496 + R.E[e].add_mag(0); //add a new magnitude for the metric
  497 +
  498 + for(unsigned p = 0; p < R.E[e].size(); p++){ //for each point in the edge
466 499  
  500 + p1 = R.E[e][p]; //get the next point in the edge
  501 + stim2array(queryPt, p1);
  502 + kdt.cpu_search(queryPt, 1, nnIdx, dists); //find the distance between A and the current network
  503 +
  504 + m1 = 1.0f - gaussianFunction((T)dists[0], sigma); //calculate the metric value based on the distance
  505 + R.E[e].set_mag(m1, p, 1); //set the error for the second point in the segment
  506 + }
  507 + }
  508 +#endif
467 509 return R; //return the resulting network
468 510 }
469 511  
... ... @@ -487,7 +529,7 @@ public:
487 529 void load_txt(std::string filename)
488 530 {
489 531 std::vector <std::string> file_contents;
490   - std::ifstream file(filename);
  532 + std::ifstream file(filename.c_str());
491 533 std::string line;
492 534 std::vector<unsigned> id2vert; //this list stores the vertex ID associated with each network vertex
493 535 //for each line in the text file, store them as strings in file_contents
... ... @@ -538,7 +580,7 @@ public:
538 580 for(unsigned int d = 0; d < 3; d++){
539 581 ss<<p[i][d];
540 582 }
541   - ss < "\n";
  583 + ss << "\n";
542 584 }
543 585 return ss.str();
544 586 }
... ... @@ -552,8 +594,8 @@ public:
552 594 void
553 595 to_txt(std::string filename)
554 596 {
555   - std::ofstream ofs(filename, std::ofstream::out | std::ofstream::app);
556   - int num;
  597 + std::ofstream ofs(filename.c_str(), std::ofstream::out | std::ofstream::app);
  598 + //int num;
557 599 ofs << (E.size()).str() << "\n";
558 600 for(unsigned int i = 0; i < E.size(); i++)
559 601 {
... ... @@ -566,7 +608,8 @@ public:
566 608 {
567 609 std::string str;
568 610 str = V[i].str();
569   - removeCharsFromString(str, "[],");
  611 + char temp[4] = "[],";
  612 + removeCharsFromString(str, temp);
570 613 ofs << str << "\n";
571 614 }
572 615 ofs.close();
... ...
stim/biomodels/network_dep.h
... ... @@ -4,7 +4,7 @@
4 4 #include <stim/math/vector.h>
5 5 #include <stim/visualization/obj.h>
6 6 #include <list>
7   -#include <ANN/ANN.h>
  7 +//#include <ANN/ANN.h>
8 8  
9 9 namespace stim{
10 10  
... ...
stim/cuda/cudatools/error.h
  1 +#ifndef STIM_CUDA_ERROR_H
  2 +#define STIM_CUDA_ERROR_H
  3 +
1 4 #include <stdio.h>
2 5 #include <iostream>
3 6 using namespace std;
4 7 #include "cuda_runtime.h"
5 8 #include "device_launch_parameters.h"
6 9 #include "cufft.h"
7   -
8   -#ifndef CUDA_HANDLE_ERROR_H
9   -#define CUDA_HANDLE_ERROR_H
  10 +#include "cublas_v2.h"
10 11  
11 12 //handle error macro
12   -static void HandleError( cudaError_t err, const char *file, int line ) {
  13 +static void cuHandleError( cudaError_t err, const char *file, int line ) {
13 14 if (err != cudaSuccess) {
14   - //FILE* outfile = fopen("cudaErrorLog.txt", "w");
15   - //fprintf(outfile, "%s in %s at line %d\n", cudaGetErrorString( err ), file, line );
16   - //fclose(outfile);
17 15 printf("%s in %s at line %d\n", cudaGetErrorString( err ), file, line );
18   - //exit( EXIT_FAILURE );
19 16  
20 17 }
21 18 }
22   -#define HANDLE_ERROR( err ) (HandleError( err, __FILE__, __LINE__ ))
  19 +#define HANDLE_ERROR( err ) (cuHandleError( err, __FILE__, __LINE__ ))
23 20  
24   -static void CufftError( cufftResult err )
  21 +static void cufftHandleError( cufftResult err, const char*file, int line )
25 22 {
26 23 if (err != CUFFT_SUCCESS)
27 24 {
... ... @@ -42,7 +39,29 @@ static void CufftError( cufftResult err )
42 39  
43 40 }
44 41 }
  42 +#define CUFFT_HANDLE_ERROR( err ) (cufftHandleError( err, __FILE__, __LINE__ ))
45 43  
  44 +static void cublasHandleError( cublasStatus_t err, const char*file, int line ){
  45 + if(err != CUBLAS_STATUS_SUCCESS){
  46 + if(err == CUBLAS_STATUS_NOT_INITIALIZED)
  47 + std::cout<<"CUBLAS_STATUS_NOT_INITIALIZED" <<" in file "<<file<<" line "<<std::endl;
  48 + else if(err == CUBLAS_STATUS_ALLOC_FAILED)
  49 + std::cout<<"CUBLAS_STATUS_ALLOC_FAILED" <<" in file "<<file<<" line "<<std::endl;
  50 + else if(err == CUBLAS_STATUS_INVALID_VALUE)
  51 + std::cout<<"CUBLAS_STATUS_INVALID_VALUE" <<" in file "<<file<<" line "<<std::endl;
  52 + else if(err == CUBLAS_STATUS_ARCH_MISMATCH)
  53 + std::cout<<"CUBLAS_STATUS_ARCH_MISMATCH" <<" in file "<<file<<" line "<<std::endl;
  54 + else if(err == CUBLAS_STATUS_MAPPING_ERROR)
  55 + std::cout<<"CUBLAS_STATUS_MAPPING_ERROR" <<" in file "<<file<<" line "<<std::endl;
  56 + else if(err == CUBLAS_STATUS_EXECUTION_FAILED)
  57 + std::cout<<"CUBLAS_STATUS_EXECUTION_FAILED" <<" in file "<<file<<" line "<<std::endl;
  58 + else if(err == CUBLAS_STATUS_INTERNAL_ERROR)
  59 + std::cout<<"CUBLAS_STATUS_INTERNAL_ERROR" <<" in file "<<file<<" line "<<std::endl;
  60 + else
  61 + std::cout<<"Unknown error"<<" in file "<<file<<" line "<<std::endl;
  62 + }
  63 +}
  64 +#define CUBLAS_HANDLE_ERROR( err ) (cublasHandleError( err, __FILE__, __LINE__ ))
46 65  
47 66  
48 67 #endif
... ...
stim/envi/agilent_binary.h
... ... @@ -4,13 +4,15 @@
4 4  
5 5 #include <string>
6 6 #include <fstream>
  7 +#include <complex>
7 8  
8 9 //CUDA
9   -#ifdef CUDA_FOUND
10   - #include <cuda_runtime.h>
11   - #include "cufft.h"
12   - #include <stim/cuda/cudatools/error.h>
13   -#endif
  10 +//#ifdef CUDA_FOUND
  11 +#include <cuda_runtime.h>
  12 +#include "cufft.h"
  13 +#include <stim/cuda/cudatools/error.h>
  14 +#include <stim/envi/envi_header.h>
  15 +//#endif
14 16  
15 17 namespace stim{
16 18  
... ... @@ -19,10 +21,10 @@ class agilent_binary{
19 21  
20 22 protected:
21 23 std::string fname;
22   - T* ptr;
23   - size_t R[3];
24   - static const size_t header = 1020;
25   - double Z[2];
  24 + T* ptr; //pointer to the image data
  25 + size_t R[3]; //size of the binary image in X, Y, and Z
  26 + static const size_t header = 1020; //header size
  27 + double Z[2]; //range of z values (position or wavelength)
26 28  
27 29 public:
28 30 size_t size(){
... ... @@ -42,6 +44,10 @@ public:
42 44 alloc();
43 45 }
44 46  
  47 + size_t dim(size_t i){
  48 + return R[i];
  49 + }
  50 +
45 51 /// Create a deep copy of an agileng_binary object
46 52 void deep_copy(agilent_binary<T>* dst, const agilent_binary<T>* src){
47 53 dst->alloc(src->R[0], src->R[1], src->R[2]); //allocate memory
... ... @@ -136,6 +142,42 @@ public:
136 142 return header;
137 143 }
138 144  
  145 + /// Subtract the mean from each pixel. Generally used for centering an interferogram.
  146 + void meancenter(){
  147 + size_t Z = R[2]; //store the number of bands
  148 + size_t XY = R[0] * R[1]; //store the number of pixels in the image
  149 + T sum = (T)0;
  150 + T mean;
  151 + for(size_t xy = 0; xy < XY; xy++){ //for each pixel
  152 + sum = 0;
  153 + for(size_t z = 0; z < Z; z++){ //for each band
  154 + sum += ptr[ z * XY + xy ]; //add the band value to a running sum
  155 + }
  156 + mean = sum / (T)Z; //calculate the pixel mean
  157 + for(size_t z = 0; z < Z; z++){
  158 + ptr[ z * XY + xy ] -= mean; //subtract the mean from each band
  159 + }
  160 + }
  161 + }
  162 +
  163 + /// adds n bands of zero padding to the end of the file
  164 + void zeropad(size_t n){
  165 + size_t newZ = R[2] + n;
  166 + T* temp = (T*) calloc(R[0] * R[1] * newZ, sizeof(T)); //allocate space for the new image
  167 + memcpy(temp, ptr, size() * sizeof(T)); //copy the old data to the new image
  168 +
  169 + free(ptr); //free the old data
  170 + ptr = temp; //swap in the new data
  171 + R[2] = newZ; //set the z-dimension to the new zero value
  172 + }
  173 +
  174 + //pads to the nearest power-of-two
  175 + void zeropad(){
  176 + size_t newZ = (size_t)pow(2, ceil(log(R[2])/log(2))); //find the nearest power-of-two
  177 + size_t n = newZ - R[2]; //calculate the number of bands to add
  178 + zeropad(n); //add the padding
  179 + }
  180 +
139 181 /// Calculate the absorbance spectrum from the transmission spectrum given a background
140 182 void absorbance(stim::agilent_binary<T>* background){
141 183 size_t N = size(); //calculate the number of values to be ratioed
... ... @@ -147,7 +189,7 @@ public:
147 189 ptr[i] = -log10(ptr[i] / background->ptr[i]);
148 190 }
149 191  
150   -#ifdef CUDA_FOUND
  192 +//#ifdef CUDA_FOUND
151 193 /// Perform an FFT and return a binary file with bands in the specified range
152 194 agilent_binary<T> fft(double band_min, double band_max, double ELWN = 15798, int UDR = 2){
153 195 auto total_start = std::chrono::high_resolution_clock::now();
... ... @@ -234,7 +276,22 @@ public:
234 276  
235 277 return result;
236 278 }
237   -#endif
  279 +
  280 + //saves the binary as an ENVI file with a BIP interleave format
  281 + int bip(T* bip_ptr){
  282 + //std::ofstream out(outfile.c_str(), std::ios::binary); //create a binary file stream for output
  283 + size_t XY = R[0] * R[1];
  284 + size_t B = R[2];
  285 + size_t b;
  286 +
  287 + for(size_t xy = 0; xy < XY; xy++){
  288 + for(b = 0; b < B; b++){
  289 + bip_ptr[xy * B + b] = ptr[b * XY + xy];
  290 + }
  291 + }
  292 + return 0;
  293 + }
  294 +//#endif
238 295  
239 296 };
240 297  
... ...
stim/envi/bil.h
... ... @@ -4,6 +4,7 @@
4 4 #include "../envi/envi_header.h"
5 5 #include "../envi/hsi.h"
6 6 #include "../math/fd_coefficients.h"
  7 +#include <stim/cuda/cudatools/error.h>
7 8 #include <cstring>
8 9 #include <utility>
9 10 #include <deque>
... ... @@ -118,7 +119,7 @@ public:
118 119 page++;
119 120 //if wavelength is larger than the last wavelength in header file
120 121 if (page == Z()) {
121   - band_index(p, Z()-1);
  122 + band_index(p, Z()-1, PROGRESS);
122 123 return true;
123 124 }
124 125 }
... ... @@ -224,10 +225,44 @@ public:
224 225 }
225 226  
226 227 //given a Y ,return a XZ slice
227   - bool read_plane_y(T * p, unsigned long long y){
  228 + bool read_plane_xz(T * p, size_t y){
228 229 return binary<T>::read_plane_2(p, y);
229 230 }
230 231  
  232 + //given a Y, return ZX slice (transposed such that the spectrum is the leading dimension)
  233 + int read_plane_zx(T* p, size_t y){
  234 + T* temp = (T*) malloc(X() * Z() * sizeof(T)); //allocate space to store the temporary xz plane
  235 + binary<T>::read_plane_2(temp, y); //load the plane from disk
  236 + size_t z, x;
  237 + for(z = 0; z < Z(); z++){
  238 + for(x = 0; x <= z; x++){
  239 + p[x * Z() + z] = temp[z * X() + x]; //copy to the destination frame
  240 + }
  241 + }
  242 + }
  243 +
  244 + //load a frame y into a pre-allocated double-precision array
  245 + int read_plane_xzd(double* f, size_t y){
  246 + size_t XB = X() * Z();
  247 + T* temp = (T*) malloc(XB * sizeof(T)); //create a temporary location to store the plane at current precision
  248 + if(!read_plane_y(temp, y)) return 1; //read the plane in its native format, if it fails return a 1
  249 + for(size_t i = 0; i < XB; i++) f[i] = temp[i]; //convert the plane to a double
  250 + return 0;
  251 + }
  252 +
  253 + //given a Y, return ZX slice (transposed such that the spectrum is the leading dimension)
  254 + int read_plane_zxd(double* p, size_t y){
  255 + T* temp = (T*) malloc(X() * Z() * sizeof(T)); //allocate space to store the temporary xz plane
  256 + binary<T>::read_plane_2(temp, y); //load the plane from disk
  257 + size_t z, x;
  258 + for(z = 0; z < Z(); z++){
  259 + for(x = 0; x < X(); x++){
  260 + p[x * Z() + z] = (double)temp[z * X() + x]; //copy to the destination frame
  261 + }
  262 + }
  263 + return 0;
  264 + }
  265 +
231 266  
232 267 /// Perform baseline correction given a list of baseline points and stores the result in a new BSQ file.
233 268  
... ... @@ -268,7 +303,7 @@ public:
268 303 for (unsigned long long k =0; k < Y(); k++)
269 304 {
270 305 //get the current y slice
271   - read_plane_y(c, k);
  306 + read_plane_xz(c, k);
272 307  
273 308 //initialize lownum, highnum, low, high
274 309 ai = w[0];
... ... @@ -369,7 +404,7 @@ public:
369 404  
370 405 for(unsigned long long j = 0; j < Y(); j++)
371 406 {
372   - read_plane_y(c, j);
  407 + read_plane_xz(c, j);
373 408 for(unsigned long long i = 0; i < B; i++)
374 409 {
375 410 for(unsigned long long m = 0; m < X(); m++)
... ... @@ -469,7 +504,7 @@ public:
469 504  
470 505 for ( unsigned long long i = 0; i < Y(); i++)
471 506 {
472   - read_plane_y(p, i);
  507 + read_plane_xz(p, i);
473 508 for ( unsigned long long k = 0; k < Z(); k++)
474 509 {
475 510 unsigned long long ks = k * X();
... ... @@ -863,7 +898,7 @@ public:
863 898  
864 899 for (unsigned long long i = 0; i < Y(); i++) //for each value in Y() (BIP should be X)
865 900 {
866   - read_plane_y(temp, i); //retrieve an ZX slice, stored in temp
  901 + read_plane_xz(temp, i); //retrieve an ZX slice, stored in temp
867 902 for ( unsigned long long j = 0; j < Z(); j++) //for each Z() (Y)
868 903 {
869 904 for (unsigned long long k = 0; k < X(); k++) //for each band
... ... @@ -933,7 +968,7 @@ public:
933 968 //for each slice along the y axis
934 969 for (unsigned long long y = 0; y < Y(); y++) //Select a page by choosing Y coordinate, Y()
935 970 {
936   - read_plane_y(slice, y); //retrieve an ZX page, store in "slice"
  971 + read_plane_xz(slice, y); //retrieve an ZX page, store in "slice"
937 972  
938 973 //for each sample along X
939 974 for (unsigned long long x = 0; x < X(); x++) //Select a pixel by choosing X coordinate in the page, X()
... ... @@ -992,43 +1027,136 @@ public:
992 1027  
993 1028 /// @param p is a pointer to pre-allocated memory of size [B * sizeof(T)] that stores the mean spectrum
994 1029 /// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
995   - bool avg_band(double* p, unsigned char* mask = NULL, bool PROGRESS = false){
  1030 + bool mean_spectrum(double* m, double* std, unsigned char* mask = NULL, bool PROGRESS = false){
996 1031 unsigned long long XZ = X() * Z();
997 1032 unsigned long long XY = X() * Y();
998 1033 T* temp = (T*)malloc(sizeof(T) * XZ);
999   - for (unsigned long long j = 0; j < Z(); j++){
1000   - p[j] = 0;
1001   - }
  1034 + memset(m, 0, Z() * sizeof(double)); //initialize the mean to zero
  1035 + double* e_x2 = (double*)malloc(Z() * sizeof(double)); //allocate space for E[x^2]
  1036 + memset(e_x2, 0, Z() * sizeof(double)); //initialize E[x^2] to zero
1002 1037 //calculate vaild number in a band
1003   - unsigned long long count = 0;
1004   - for (unsigned long long j = 0; j < XY; j++){
1005   - if (mask == NULL || mask[j] != 0){
1006   - count++;
1007   - }
1008   - }
  1038 + size_t count = nnz(mask); //count the number of pixels in the mask
  1039 +
  1040 + double x; //create a register to store the pixel value
1009 1041 for (unsigned long long k = 0; k < Y(); k++){
1010   - read_plane_y(temp, k);
  1042 + read_plane_xz(temp, k);
1011 1043 unsigned long long kx = k * X();
1012 1044 for (unsigned long long i = 0; i < X(); i++){
1013 1045 if (mask == NULL || mask[kx + i] != 0){
1014 1046 for (unsigned long long j = 0; j < Z(); j++){
1015   - p[j] += temp[j * X() + i] / (double)count;
  1047 + x = temp[j * X() + i];
  1048 + m[j] += x / (double)count;
  1049 + e_x2[j] += x*x / (double)count;
1016 1050 }
1017 1051 }
1018 1052 }
1019 1053 if(PROGRESS) progress = (double)(k+1) / Y() * 100;
1020 1054 }
  1055 +
  1056 + for(size_t i = 0; i < Z(); i++) //calculate the standard deviation
  1057 + std[i] = sqrt(e_x2[i] - m[i] * m[i]);
  1058 +
1021 1059 free(temp);
1022 1060 return true;
1023 1061 }
1024 1062  
  1063 + int co_matrix_cublas(double* co, double* avg, unsigned char *mask, bool PROGRESS = false){
  1064 + cublasStatus_t stat;
  1065 + cublasHandle_t handle;
  1066 +
  1067 + progress = 0; //initialize the progress to zero (0)
  1068 + size_t XY = X() * Y(); //calculate the number of elements in a band image
  1069 + size_t XB = X() * Z();
  1070 + size_t B = Z(); //calculate the number of spectral elements
  1071 +
  1072 + double* F = (double*)malloc(sizeof(double) * B * X()); //allocate space for the frame that will be pulled from the file
  1073 + double* F_dev;
  1074 + HANDLE_ERROR(cudaMalloc(&F_dev, X() * B * sizeof(double))); //allocate space for the frame on the GPU
  1075 + double* s_dev; //declare a device pointer that will store the spectrum on the GPU
  1076 + double* A_dev; //declare a device pointer that will store the covariance matrix on the GPU
  1077 + double* avg_dev; //declare a device pointer that will store the average spectrum
  1078 + HANDLE_ERROR(cudaMalloc(&s_dev, B * sizeof(double))); //allocate space on the CUDA device for a spectrum
  1079 + HANDLE_ERROR(cudaMalloc(&A_dev, B * B * sizeof(double))); //allocate space on the CUDA device for the covariance matrix
  1080 + HANDLE_ERROR(cudaMemset(A_dev, 0, B * B * sizeof(double))); //initialize the covariance matrix to zero (0)
  1081 + HANDLE_ERROR(cudaMalloc(&avg_dev, XB * sizeof(double))); //allocate space on the CUDA device for the average spectrum
  1082 + for(size_t x = 0; x < X(); x++) //make multiple copies of the average spectrum in order to build a matrix
  1083 + HANDLE_ERROR(cudaMemcpy(&avg_dev[x * B], avg, B * sizeof(double), cudaMemcpyHostToDevice));
  1084 + //stat = cublasSetVector((int)B, sizeof(double), avg, 1, avg_dev, 1); //copy the average spectrum to the CUDA device
  1085 +
  1086 + double ger_alpha = 1.0/(double)XY; //scale the outer product by the inverse of the number of samples (mean outer product)
  1087 + double axpy_alpha = -1; //multiplication factor for the average spectrum (in order to perform a subtraction)
  1088 +
  1089 + CUBLAS_HANDLE_ERROR(stat = cublasCreate(&handle)); //create a cuBLAS instance
  1090 + if (stat != CUBLAS_STATUS_SUCCESS) return 1; //test the cuBLAS instance to make sure it is valid
  1091 +
  1092 + else std::cout<<"Using cuBLAS to calculate the mean covariance matrix..."<<std::endl;
  1093 + double beta = 1.0;
  1094 + size_t x, y;
  1095 + for(y = 0; y < Y(); y++){ //for each line
  1096 + read_plane_zxd(F, y); //read a frame from the file
  1097 + HANDLE_ERROR(cudaMemcpy(F_dev, F, XB * sizeof(double), cudaMemcpyHostToDevice)); //copy the frame to the GPU
  1098 + CUBLAS_HANDLE_ERROR(cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, (int)B, (int)X(), &axpy_alpha, avg_dev, (int)B, &beta, F_dev, (int)B, F_dev, (int)B));//subtract the mean spectrum
  1099 +
  1100 + for(x = 0; x < X(); x++)
  1101 + CUBLAS_HANDLE_ERROR(cublasDsyr(handle, CUBLAS_FILL_MODE_UPPER, (int)B, &ger_alpha, &F_dev[x*B], 1, A_dev, (int)B)); //perform an outer product
  1102 + if(PROGRESS) progress = (double)(y + 1) / Y() * 100;
  1103 + }
  1104 +
  1105 + cublasGetMatrix((int)B, (int)B, sizeof(double), A_dev, (int)B, co, (int)B); //copy the result from the GPU to the CPU
  1106 +
  1107 + cudaFree(A_dev); //clean up allocated device memory
  1108 + cudaFree(s_dev);
  1109 + cudaFree(avg_dev);
  1110 +
  1111 + for(unsigned long long i = 0; i < B; i++){ //copy the upper triangular portion to the lower triangular portion
  1112 + for(unsigned long long j = i+1; j < B; j++){
  1113 + co[B * i + j] = co[B * j + i];
  1114 + }
  1115 + }
  1116 +
  1117 + return 0;
  1118 +
  1119 +
  1120 +
  1121 + }
  1122 +
  1123 +
1025 1124 /// Calculate the covariance matrix for all masked pixels in the image.
1026 1125  
1027 1126 /// @param co is a pointer to pre-allocated memory of size [B * B] that stores the resulting covariance matrix
1028 1127 /// @param avg is a pointer to memory of size B that stores the average spectrum
1029 1128 /// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
1030   - bool co_matrix(double* co, double* avg, unsigned char *mask, bool PROGRESS = false){
  1129 + bool co_matrix(double* co, double* avg, unsigned char *mask, bool use_gpu = true, bool PROGRESS = false){
1031 1130 progress = 0;
  1131 +
  1132 + if(use_gpu){
  1133 + int dev_count;
  1134 + HANDLE_ERROR(cudaGetDeviceCount(&dev_count)); //get the number of CUDA devices
  1135 + std::cout<<"Number of CUDA devices: "<<dev_count<<std::endl; //output the number of CUDA devices
  1136 + cudaDeviceProp prop;
  1137 + int best_device_id = 0; //stores the best CUDA device
  1138 + float best_device_cc = 0.0f; //stores the compute capability of the best device
  1139 + std::cout<<"CUDA devices:"<<std::endl;
  1140 + for(int d = 0; d < dev_count; d++){ //for each CUDA device
  1141 + cudaGetDeviceProperties(&prop, d); //get the property of the first device
  1142 + float cc = prop.major + prop.minor / 10.0f; //calculate the compute capability
  1143 + std::cout<<"("<<prop.major<<"."<<prop.minor<<") "<<prop.name<<std::endl; //display the device information
  1144 + if(cc > best_device_cc){
  1145 + best_device_cc = cc; //if this is better than the previous device, use it
  1146 + best_device_id = d;
  1147 + }
  1148 + }
  1149 +
  1150 + if(dev_count > 0 && prop.major != 9999){ //if the first device is not an emulator
  1151 + std::cout<<"Using device "<<best_device_id<<std::endl;
  1152 + HANDLE_ERROR(cudaSetDevice(best_device_id));
  1153 + int status = co_matrix_cublas(co, avg, mask, PROGRESS); //use cuBLAS to calculate the covariance matrix
  1154 + if(status == 0) return true; //if the cuBLAS function returned correctly, we're done
  1155 + } //otherwise continue using the CPU
  1156 +
  1157 + std::cout<<"No supported CUDA devices found or cuBLAS failed. Using CPU"<<std::endl;
  1158 + }
  1159 +
1032 1160 //memory allocation
1033 1161 unsigned long long xy = X() * Y();
1034 1162 unsigned long long B = Z();
... ... @@ -1325,7 +1453,7 @@ public:
1325 1453 c = (T*)malloc( L ); //allocate space for the slice
1326 1454  
1327 1455 for(unsigned long long j = 0; j < Y(); j++){ //for each line
1328   - read_plane_y(c, j); //load the line into memory
  1456 + read_plane_xz(c, j); //load the line into memory
1329 1457 for(unsigned long long i = 0; i < B; i++){ //for each band
1330 1458 for(unsigned long long m = 0; m < X(); m++){ //for each sample
1331 1459 if( mask == NULL && mask[m + j * X()] ) //if the pixel is masked
... ... @@ -1355,7 +1483,7 @@ public:
1355 1483 c = (T*)malloc( L ); //allocate space for the slice
1356 1484  
1357 1485 for(unsigned long long j = 0; j < Y(); j++){ //for each line
1358   - read_plane_y(c, j); //load the line into memory
  1486 + read_plane_xz(c, j); //load the line into memory
1359 1487 for(unsigned long long i = 0; i < B; i++){ //for each band
1360 1488 for(unsigned long long m = 0; m < X(); m++){ //for each sample
1361 1489 if( mask == NULL && mask[m + j * X()] ) //if the pixel is masked
... ...
stim/envi/bip.h
... ... @@ -5,13 +5,16 @@
5 5 #include "../envi/bil.h"
6 6 #include "../envi/hsi.h"
7 7 #include <cstring>
  8 +#include <complex>
8 9 #include <utility>
9 10  
10 11 //CUDA
11   -#ifdef CUDA_FOUND
12   - #include <cuda_runtime.h>
13   - #include "cublas_v2.h"
14   -#endif
  12 +//#ifdef CUDA_FOUND
  13 +#include <stim/cuda/cudatools/error.h>
  14 +#include <cuda_runtime.h>
  15 +#include "cublas_v2.h"
  16 +#include "cufft.h"
  17 +//#endif
15 18  
16 19 namespace stim{
17 20  
... ... @@ -257,7 +260,7 @@ public:
257 260 }
258 261  
259 262 //given a Y ,return a ZX slice
260   - bool read_plane_y(T * p, unsigned long long y){
  263 + bool read_plane_y(T * p, size_t y){
261 264 return binary<T>::read_plane_2(p, y);
262 265 }
263 266  
... ... @@ -954,39 +957,43 @@ public:
954 957  
955 958 /// @param p is a pointer to pre-allocated memory of size [B * sizeof(T)] that stores the mean spectrum
956 959 /// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
957   - bool avg_band(double* p, unsigned char* mask = NULL, bool PROGRESS = false){
  960 + bool mean_spectrum(double* m, double* std, unsigned char* mask = NULL, bool PROGRESS = false){
958 961 unsigned long long XY = X() * Y(); //calculate the total number of pixels in the HSI
959 962 T* temp = (T*)malloc(sizeof(T) * Z()); //allocate space for the current spectrum to be read
960   - memset(p, 0, sizeof(double) * Z()); //initialize the average spectrum to zero (0)
961   - //for (unsigned j = 0; j < Z(); j++){
962   - // p[j] = 0;
963   - //}
  963 + memset(m, 0, Z() * sizeof(double)); //set the mean spectrum to zero
  964 + double* e_x2 = (double*)malloc(Z() * sizeof(double)); //allocate space for E[x^2]
  965 + memset(e_x2, 0, Z() * sizeof(double)); //set all values for E[x^2] to zero
964 966  
965 967 unsigned long long count = nnz(mask); //calculate the number of masked pixels
966   -
  968 + double x;
967 969 for (unsigned long long i = 0; i < XY; i++){ //for each pixel in the HSI
968 970 if (mask == NULL || mask[i] != 0){ //if the pixel is masked
969 971 pixel(temp, i); //get the spectrum
970 972 for (unsigned long long j = 0; j < Z(); j++){ //for each spectral component
971   - p[j] += (double)temp[j] / (double)count; //add the weighted value to the average
  973 + x = temp[j];
  974 + m[j] += x / (double)count; //add the weighted value to the average
  975 + e_x2[j] += x*x / (double)count;
972 976 }
973 977 }
974 978 if(PROGRESS) progress = (double)(i+1) / XY * 100; //increment the progress
975 979 }
976 980  
  981 + //calculate the standard deviation
  982 + for(size_t i = 0; i < Z(); i++)
  983 + std[i] = sqrt(e_x2[i] - m[i] * m[i]);
  984 +
977 985 free(temp);
978 986 return true;
979 987 }
980   -#ifdef CUDA_FOUND
  988 +//#ifdef CUDA_FOUND
981 989 /// Calculate the covariance matrix for masked pixels using cuBLAS
982 990 /// Note that cuBLAS only supports integer-sized arrays, so there may be issues with large spectra
983   - bool co_matrix_cublas(double* co, double* avg, unsigned char *mask, bool PROGRESS = false){
  991 + int co_matrix_cublas(double* co, double* avg, unsigned char *mask, bool PROGRESS = false){
984 992  
985 993 cudaError_t cudaStat;
986 994 cublasStatus_t stat;
987 995 cublasHandle_t handle;
988 996  
989   - progress = 0; //initialize the progress to zero (0)
990 997 unsigned long long XY = X() * Y(); //calculate the number of elements in a band image
991 998 unsigned long long B = Z(); //calculate the number of spectral elements
992 999  
... ... @@ -1004,10 +1011,9 @@ public:
1004 1011 double axpy_alpha = -1; //multiplication factor for the average spectrum (in order to perform a subtraction)
1005 1012  
1006 1013 stat = cublasCreate(&handle); //create a cuBLAS instance
1007   - if (stat != CUBLAS_STATUS_SUCCESS) { //test the cuBLAS instance to make sure it is valid
1008   - printf ("CUBLAS initialization failed\n");
1009   - return EXIT_FAILURE;
1010   - }
  1014 + if (stat != CUBLAS_STATUS_SUCCESS) return 1; //test the cuBLAS instance to make sure it is valid
  1015 +
  1016 + else std::cout<<"Using cuBLAS to calculate the mean covariance matrix..."<<std::endl;
1011 1017 for (unsigned long long xy = 0; xy < XY; xy++){ //for each pixel
1012 1018 if (mask == NULL || mask[xy] != 0){
1013 1019 pixeld(s, xy); //retreive the spectrum at the current xy pixel location
... ... @@ -1031,26 +1037,45 @@ public:
1031 1037 }
1032 1038 }
1033 1039  
1034   - return true;
  1040 + return 0;
1035 1041 }
1036   -#endif
  1042 +//#endif
1037 1043  
1038 1044 /// Calculate the covariance matrix for all masked pixels in the image with 64-bit floating point precision.
1039 1045  
1040 1046 /// @param co is a pointer to pre-allocated memory of size [B * B] that stores the resulting covariance matrix
1041 1047 /// @param avg is a pointer to memory of size B that stores the average spectrum
1042 1048 /// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
1043   - bool co_matrix(double* co, double* avg, unsigned char *mask, bool PROGRESS = false){
1044   -
1045   -#ifdef CUDA_FOUND
1046   - int dev_count;
1047   - cudaGetDeviceCount(&dev_count); //get the number of CUDA devices
1048   - cudaDeviceProp prop;
1049   - cudaGetDeviceProperties(&prop, 0); //get the property of the first device
1050   - if(dev_count > 0 && prop.major != 9999) //if the first device is not an emulator
1051   - return co_matrix_cublas(co, avg, mask, PROGRESS); //use cuBLAS to calculate the covariance matrix
1052   -#endif
  1049 + bool co_matrix(double* co, double* avg, unsigned char *mask, bool use_gpu = true, bool PROGRESS = false){
1053 1050 progress = 0;
  1051 +
  1052 + if(use_gpu){
  1053 + int dev_count;
  1054 + HANDLE_ERROR(cudaGetDeviceCount(&dev_count)); //get the number of CUDA devices
  1055 + std::cout<<"Number of CUDA devices: "<<dev_count<<std::endl; //output the number of CUDA devices
  1056 + cudaDeviceProp prop;
  1057 + int best_device_id = 0; //stores the best CUDA device
  1058 + float best_device_cc = 0.0f; //stores the compute capability of the best device
  1059 + std::cout<<"CUDA devices----"<<std::endl;
  1060 + for(int d = 0; d < dev_count; d++){ //for each CUDA device
  1061 + cudaGetDeviceProperties(&prop, d); //get the property of the first device
  1062 + float cc = prop.major + prop.minor / 10.0f; //calculate the compute capability
  1063 + std::cout<<d<<": ["<<prop.major<<"."<<prop.minor<<"] "<<prop.name<<std::endl; //display the device information
  1064 + if(cc > best_device_cc){
  1065 + best_device_cc = cc; //if this is better than the previous device, use it
  1066 + best_device_id = d;
  1067 + }
  1068 + }
  1069 +
  1070 + if(dev_count > 0 && prop.major != 9999){ //if the first device is not an emulator
  1071 + std::cout<<"Using device "<<best_device_id<<std::endl;
  1072 + HANDLE_ERROR(cudaSetDevice(best_device_id));
  1073 + int status = co_matrix_cublas(co, avg, mask, PROGRESS); //use cuBLAS to calculate the covariance matrix
  1074 + if(status == 0) return true; //if the cuBLAS function returned correctly, we're done
  1075 + } //otherwise continue using the CPU
  1076 +
  1077 + std::cout<<"No supported CUDA devices found or cuBLAS failed. Using CPU"<<std::endl;
  1078 + }
1054 1079 //memory allocation
1055 1080 unsigned long long XY = X() * Y();
1056 1081 unsigned long long B = Z();
... ... @@ -1092,10 +1117,10 @@ public:
1092 1117 }
1093 1118  
1094 1119  
1095   -#ifdef CUDA_FOUND
  1120 +//#ifdef CUDA_FOUND
1096 1121 /// Calculate the covariance matrix of Noise for masked pixels using cuBLAS
1097 1122 /// Note that cuBLAS only supports integer-sized arrays, so there may be issues with large spectra
1098   - bool coNoise_matrix_cublas(double* coN, double* avg, unsigned char *mask, bool PROGRESS = false){
  1123 + int coNoise_matrix_cublas(double* coN, double* avg, unsigned char *mask, bool PROGRESS = false){
1099 1124  
1100 1125 cudaError_t cudaStat;
1101 1126 cublasStatus_t stat;
... ... @@ -1123,11 +1148,9 @@ public:
1123 1148 double ger_alpha = 1.0/(double)XY; //scale the outer product by the inverse of the number of samples (mean outer product)
1124 1149 double axpy_alpha = -1; //multiplication factor for the average spectrum (in order to perform a subtraction)
1125 1150  
1126   - stat = cublasCreate(&handle); //create a cuBLAS instance
1127   - if (stat != CUBLAS_STATUS_SUCCESS) { //test the cuBLAS instance to make sure it is valid
1128   - printf ("CUBLAS initialization failed\n");
1129   - return EXIT_FAILURE;
1130   - }
  1151 + CUBLAS_HANDLE_ERROR(cublasCreate(&handle)); //create a cuBLAS instance
  1152 + if (stat != CUBLAS_STATUS_SUCCESS) return 1; //test the cuBLAS instance to make sure it is valid
  1153 +
1131 1154 for (unsigned long long xy = 0; xy < XY; xy++){ //for each pixel
1132 1155 if (mask == NULL || mask[xy] != 0){
1133 1156 pixeld(s, xy); //retreive the spectrum at the current xy pixel location
... ... @@ -1158,27 +1181,44 @@ public:
1158 1181 }
1159 1182 }
1160 1183  
1161   - return true;
  1184 + return 0;
1162 1185 }
1163   -#endif
  1186 +//#endif
1164 1187  
1165 1188 /// Calculate the covariance of noise matrix for all masked pixels in the image with 64-bit floating point precision.
1166 1189  
1167 1190 /// @param coN is a pointer to pre-allocated memory of size [B * B] that stores the resulting covariance matrix
1168 1191 /// @param avg is a pointer to memory of size B that stores the average spectrum
1169 1192 /// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
1170   - bool coNoise_matrix(double* coN, double* avg, unsigned char *mask, bool PROGRESS = false){
1171   -
1172   -#ifdef CUDA_FOUND
1173   - int dev_count;
1174   - cudaGetDeviceCount(&dev_count); //get the number of CUDA devices
1175   - cudaDeviceProp prop;
1176   - cudaGetDeviceProperties(&prop, 0); //get the property of the first device
1177   - if(dev_count > 0 && prop.major != 9999) //if the first device is not an emulator
1178   - return coNoise_matrix_cublas(coN, avg, mask, PROGRESS); //use cuBLAS to calculate the covariance matrix
1179   -#endif
1180   -
1181   -
  1193 + bool coNoise_matrix(double* coN, double* avg, unsigned char *mask, bool use_gpu = true, bool PROGRESS = false){
  1194 +
  1195 + if(use_gpu){
  1196 + int dev_count;
  1197 + HANDLE_ERROR(cudaGetDeviceCount(&dev_count)); //get the number of CUDA devices
  1198 + std::cout<<"Number of CUDA devices: "<<dev_count<<std::endl; //output the number of CUDA devices
  1199 + cudaDeviceProp prop;
  1200 + int best_device_id = 0; //stores the best CUDA device
  1201 + float best_device_cc = 0.0f; //stores the compute capability of the best device
  1202 + std::cout<<"CUDA devices:"<<std::endl;
  1203 + for(int d = 0; d < dev_count; d++){ //for each CUDA device
  1204 + cudaGetDeviceProperties(&prop, d); //get the property of the first device
  1205 + float cc = prop.major + prop.minor / 10.0f; //calculate the compute capability
  1206 + std::cout<<d<<": ("<<prop.major<<"."<<prop.minor<<") "<<prop.name<<std::endl; //display the device information
  1207 + if(cc > best_device_cc){
  1208 + best_device_cc = cc; //if this is better than the previous device, use it
  1209 + best_device_id = d;
  1210 + }
  1211 + }
  1212 +
  1213 + if(dev_count > 0 && prop.major != 9999){ //if the first device is not an emulator
  1214 + std::cout<<"Using device "<<best_device_id<<std::endl;
  1215 + HANDLE_ERROR(cudaSetDevice(best_device_id));
  1216 + int status = coNoise_matrix_cublas(coN, avg, mask, PROGRESS); //use cuBLAS to calculate the covariance matrix
  1217 + if(status == 0) return true; //if the cuBLAS function returned correctly, we're done
  1218 + } //otherwise continue using the CPU
  1219 +
  1220 + std::cout<<"cuBLAS initialization failed - using CPU"<<std::endl;
  1221 + }
1182 1222  
1183 1223 progress = 0;
1184 1224 //memory allocation
... ... @@ -1443,7 +1483,7 @@ public:
1443 1483 unsigned long long jump_sample = ( (Z() - b1) + b0 ) * sizeof(T);
1444 1484  
1445 1485 //distance between sample spectra in adjacent lines
1446   - unsigned long long jump_line = (X() - x1) * Z() * sizeof(T);
  1486 + unsigned long long jump_line = ( X() - x1 + x0 ) * Z() * sizeof(T);
1447 1487  
1448 1488  
1449 1489 //unsigned long long sp = y0 * X() + x0; //start pixel
... ... @@ -1682,7 +1722,117 @@ public:
1682 1722 return true;
1683 1723 }
1684 1724  
  1725 + int fft(std::string outname, size_t bandmin, size_t bandmax, size_t samples = 0, T* ratio = NULL, size_t rx = 0, size_t ry = 0, bool PROGRESS = false, int device = 0){
  1726 + if(device == -1){
  1727 + std::cout<<"ERROR: GPU required for FFT (uses cuFFT)."<<std::endl;
  1728 + exit(1);
  1729 + }
  1730 + if(samples == 0) samples = Z(); //if samples are specified, use all of them
  1731 + if(samples > Z()){
  1732 + std::cout<<"ERROR: stim::envi doesn't support FFT padding just yet."<<std::endl;
  1733 + exit(1);
  1734 + }
  1735 + int nd; //stores the number of CUDA devices
  1736 + HANDLE_ERROR(cudaGetDeviceCount(&nd)); //get the number of CUDA devices
  1737 + if(device >= nd){ //test for the existence of the requested device
  1738 + std::cout<<"ERROR: requested CUDA device for stim::envi::fft() doesn't exist"<<std::endl;
  1739 + exit(1);
  1740 + }
  1741 + HANDLE_ERROR(cudaSetDevice(device)); //set the CUDA device
  1742 + cudaDeviceProp prop;
  1743 + HANDLE_ERROR(cudaGetDeviceProperties(&prop, device)); //get the CUDA device properties
  1744 +
  1745 + size_t B = Z();
  1746 + size_t S = samples;
  1747 + size_t fft_size = S * sizeof(T); //number of bytes for each FFT