visualization / tensorview

Commit e4667cee90ab54837ae9450fcd788539c0a4aa0f

Authored by David Mayerich
0 parents

first commit

Showing 7 changed files with 902 additions and 0 deletions   Show diff stats
Inline Side-by-side
CMakeLists.txt 0 → 100644
  Show/Hide comments   View file @e4667ce
  1 +++ a/CMakeLists.txt
  1 +#Specify the version being used aswell as the language
  2 +cmake_minimum_required(VERSION 3.12)
  3 +
  4 +#Name your project here
  5 +project(tensorview)
  6 +
  7 +#set the module directory
  8 +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}")
  9 +
  10 +#default to release mode
  11 +if(NOT CMAKE_BUILD_TYPE)
  12 + set(CMAKE_BUILD_TYPE Release)
  13 +endif(NOT CMAKE_BUILD_TYPE)
  14 +
  15 +#build the executable in the binary directory on MS Visual Studio
  16 +if ( MSVC )
  17 + SET( CMAKE_RUNTIME_OUTPUT_DIRECTORY_DEBUG "${OUTPUT_DIRECTORY}")
  18 + SET( CMAKE_RUNTIME_OUTPUT_DIRECTORY_RELEASE "${OUTPUT_DIRECTORY}")
  19 + SET( LIBRARY_OUTPUT_DIRECTORY_DEBUG "${OUTPUT_DIRECTORY}")
  20 + SET( LIBRARY_OUTPUT_DIRECTORY_RELEASE "${OUTPUT_DIRECTORY}")
  21 + add_definitions(-D_CRT_SECURE_NO_WARNINGS)
  22 + add_definitions(-D_SCL_SECURE_NO_WARNINGS)
  23 +endif ( MSVC )
  24 +#MAYBE REMOVE-----------------
  25 +#set C++11 flags if using GCC
  26 +if( CMAKE_COMPILER_IS_GNUCC )
  27 + SET( CMAKE_CXX_FLAGS "-std=c++11")
  28 + SET( CUDA_NVCC_FLAGS "-std=c++11")
  29 +endif( CMAKE_COMPILER_IS_GNUCC )
  30 +#-----------------------------
  31 +
  32 +#find packages-----------------------------------
  33 +
  34 +#find the pthreads package
  35 +find_package(Threads)
  36 +
  37 +#find the X11 package
  38 +find_package(X11)
  39 +
  40 +find_package(GLFW REQUIRED)
  41 +if(WIN32)
  42 + find_package(GLEW REQUIRED)
  43 + include_directories(${GLEW_INCLUDE_DIR})
  44 +endif(WIN32)
  45 +
  46 +#find and set up OpenGL
  47 +find_package(OpenGL REQUIRED)
  48 +
  49 +find_package(STIM REQUIRED)
  50 +
  51 +add_executable(tensorview
  52 + main.cpp
  53 + tensorfield.h
  54 +)
  55 +
  56 +#include include directories
  57 +include_directories(${GLFW_INCLUDE_DIRS}
  58 + ${OpenGL_INCLUDE_DIRS}
  59 + ${STIM_INCLUDE_DIRS}
  60 +)
  61 +
  62 +target_link_libraries(tensorview
  63 + ${GLFW_LIBRARIES}
  64 + ${OPENGL_gl_LIBRARY}
  65 + ${OPENGL_glu_LIBRARY}
  66 +)
  67 +if(WIN32)
  68 + target_link_libraries(tensorview ${GLEW_GLEW_LIBRARY})
  69 +endif(WIN32)
0 70 \ No newline at end of file
... ...
FindGLEW.cmake 0 → 100644
  Show/Hide comments   View file @e4667ce
  1 +++ a/FindGLEW.cmake
  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 ()
... ...
FindGLFW.cmake 0 → 100644
  Show/Hide comments   View file @e4667ce
  1 +++ a/FindGLFW.cmake
  1 +# Find GLFW 3
  2 +#
  3 +# GLFW_LIBRARIES
  4 +# GLFW_INCLUDE_DIRS.
  5 +# GLFW_FOUND
  6 +set(GLFW_ROOT $ENV{GLFW_ROOT})
  7 +IF(NOT UNIX)
  8 + IF(NOT GLFW_ROOT)
  9 + MESSAGE("ERROR: GLFW_ROOT must be set!")
  10 + ENDIF(NOT GLFW_ROOT)
  11 +
  12 + FIND_PATH(GLFW_INCLUDE_DIRS DOC "Path to GLFW include directory."
  13 + NAMES GLFW/glfw3.h
  14 + PATHS ${GLFW_ROOT}/include)
  15 + IF(MSVC15)
  16 + FIND_LIBRARY(GLFW_LIBRARIES DOC "Absolute path to GLFW library."
  17 + NAMES glfw3.lib
  18 + PATHS ${GLFW_ROOT}/lib-vc2015)
  19 + ELSEIF(MSVC13)
  20 + FIND_LIBRARY(GLFW_LIBRARIES DOC "Absolute path to GLFW library."
  21 + NAMES glfw3.lib
  22 + PATHS ${GLFW_ROOT}/lib-vc2013)
  23 + ELSEIF(MSVC12)
  24 + FIND_LIBRARY(GLFW_LIBRARIES DOC "Absolute path to GLFW library."
  25 + NAMES glfw3.lib
  26 + PATHS ${GLFW_ROOT}/lib-vc2012)
  27 + ELSEIF(MSVC10)
  28 + FIND_LIBRARY(GLFW_LIBRARIES DOC "Absolute path to GLFW library."
  29 + NAMES glfw3.lib
  30 + PATHS ${GLFW_ROOT}/lib-vc2010)
  31 + ELSEIF(MINGW)
  32 + IF(CMAKE_CL_64)
  33 + FIND_LIBRARY(GLFW_LIBRARIES DOC "Absolute path to GLFW library."
  34 + NAMES glfw3.dll
  35 + PATHS ${GLFW_ROOT}/lib-mingw-w64)
  36 + ELSE(CMAKE_CL_64)
  37 + FIND_LIBRARY(GLFW_LIBRARIES DOC "Absolute path to GLFW library."
  38 + NAMES glfw3.dll
  39 + PATHS ${GLFW_ROOT}/lib-mingw)
  40 + ENDIF(CMAKE_CL_64)
  41 + ELSE(MINGW)
  42 + # Default to latest version of VC libs
  43 + FIND_LIBRARY(GLFW_LIBRARIES DOC "Absolute path to GLFW library."
  44 + NAMES glfw3.lib
  45 + PATHS ${GLFW_ROOT}/lib-vc2015)
  46 + ENDIF(MSVC15)
  47 +ELSE(NOT UNIX)
  48 + FIND_PATH(GLFW_INCLUDE_DIRS DOC "Path to GLFW include directory."
  49 + NAMES GLFW/glfw3.h
  50 + PATHS
  51 + /usr/include
  52 + /usr/local/include
  53 + /usr/target/include
  54 + /sw/include
  55 + /opt/local/include)
  56 +
  57 + FIND_LIBRARY(GLFW_LIBRARIES DOC "Absolute path to GLFW library."
  58 + NAMES glfw3.dll
  59 + PATHS
  60 + /usr/local/lib
  61 + /usr/lib
  62 + /lib)
  63 +ENDIF(NOT UNIX)
  64 +
  65 +include(FindPackageHandleStandardArgs)
  66 +find_package_handle_standard_args(GLFW DEFAULT_MSG GLFW_LIBRARIES GLFW_INCLUDE_DIRS)
  67 +
  68 +mark_as_advanced(GLFW_INCLUDE_DIRS GLFW_LIBRARIES)
... ...
FindSTIM.cmake 0 → 100644
  Show/Hide comments   View file @e4667ce
  1 +++ a/FindSTIM.cmake
  1 +# finds the STIM library (downloads it if it isn't present)
  2 +# set STIMLIB_PATH to the directory containing the stim subdirectory (the stim repository)
  3 +
  4 +include(FindPackageHandleStandardArgs)
  5 +
  6 +set(STIM_ROOT $ENV{STIM_ROOT})
  7 +
  8 +IF(NOT UNIX)
  9 + IF(NOT STIM_ROOT)
  10 + MESSAGE("ERROR: STIM_ROOT environment variable must be set!")
  11 + ENDIF(NOT STIM_ROOT)
  12 +
  13 + FIND_PATH(STIM_INCLUDE_DIRS DOC "Path to GLFW include directory."
  14 + NAMES stim/image/image.h
  15 + PATHS ${STIM_ROOT})
  16 +ENDIF(NOT UNIX)
  17 +
  18 +find_package_handle_standard_args(STIM DEFAULT_MSG STIM_INCLUDE_DIRS)
  19 +
  20 +if(STIM_FOUND)
  21 + set(STIM_INCLUDE_DIRS ${STIM_INCLUDE_DIRS})
  22 +elseif(STIM_FOUND)
  23 + message("STIM library not found. Set the STIM_ROOT environment variable to the STIM location.")
  24 + message("STIMLIB can be found here: https://git.stim.ee.uh.edu/codebase/stimlib")
  25 +endif(STIM_FOUND)
  26 +
  27 +find_package_handle_standard_args(STIM DEFAULT_MSG STIM_INCLUDE_DIRS)
... ...
main.cpp 0 → 100644
  Show/Hide comments   View file @e4667ce
  1 +++ a/main.cpp
  1 +#include <GL/glew.h>
  2 +#include <GLFW/glfw3.h>
  3 +
  4 +#include <iostream>
  5 +#include <vector>
  6 +const double PI = 3.14159265358979323846;
  7 +
  8 +#include "tensorfield.h"
  9 +#include <stim/math/quaternion.h>
  10 +//#include <stim/>
  11 +
  12 +// Global variable to store the GLFW window
  13 +GLFWwindow* window;
  14 +tensorfield<float> T;
  15 +size_t z_slice = 0;
  16 +bool cout_frame = false;
  17 +
  18 +float diffuse_intensity = 0.7f;
  19 +float ambient_intensity = 0.3f;
  20 +
  21 +int glyph_type = 1;
  22 +float gamma = 3;
  23 +int glyph_resolution = 10;
  24 +float glyph_scale = 0.5;
  25 +
  26 +//glyph vertices
  27 +bool glyph_calculated = false;
  28 +
  29 +struct float3 {
  30 + float x;
  31 + float y;
  32 + float z;
  33 +};
  34 +
  35 +struct float2 {
  36 + float u;
  37 + float v;
  38 +};
  39 +
  40 +std::vector<float3> vertices;
  41 +std::vector<float3> normals;
  42 +std::vector<float2> texcoords;
  43 +std::vector<float> sin_theta;
  44 +std::vector<float> cos_theta;
  45 +std::vector<float> sin_phi;
  46 +std::vector<float> cos_phi;
  47 +
  48 +struct eigendecomposition {
  49 +
  50 + stim::vec3<float> v0;
  51 + stim::vec3<float> v2;
  52 + stim::vec3<float> lambda;
  53 +};
  54 +
  55 +inline float sgn(float x) {
  56 + if (x < 0)
  57 + return -1;
  58 + else if (x > 0)
  59 + return 1;
  60 + return 0;
  61 +}
  62 +
  63 +stim::vec3<float> triangle_norm(stim::vec3<float> p0, stim::vec3<float> p1, stim::vec3<float> p2) {
  64 + stim::vec3<float> a = p1 - p0;
  65 + stim::vec3<float> b = p2 - p0;
  66 + stim::vec3<float> n = a.cross(b);
  67 + return n;
  68 +}
  69 +
  70 +
  71 +// returns the fractional anisotropy and stores the spherical, linear, and planar components in cs, cl, and cp
  72 +float get_anisotropy(stim::vec3<float> lambda, float& cs, float& cl, float& cp) {
  73 +
  74 + //calculate the denominator for each specific anisotropy
  75 + float denom = lambda[0] + lambda[1] + lambda[2];
  76 +
  77 + cl = (lambda[2] - lambda[1]) / denom;
  78 + cp = 2 * (lambda[1] - lambda[0]) / denom;
  79 + cs = 3 * lambda[0] / denom;
  80 +
  81 + return 0;
  82 +}
  83 +
  84 +//calculate the eigenvectors and eigenvalues for the input pixel (x, y, z)
  85 +// eigenvalues are embedded in the length of the eigenvector
  86 +void get_pixel_eigen(eigendecomposition* e, size_t x, size_t y, size_t z) {
  87 +
  88 + e->v0[0] = T(x, y, z, 0, 0);
  89 + e->v0[1] = T(x, y, z, 1, 0);
  90 + e->v0[2] = T(x, y, z, 2, 0);
  91 +
  92 + //temporarily store v1 in order to get the eigenvalue
  93 + // (the eigenvector is redundant because it's the cross product of v0 and v1)
  94 + stim::vec3<float> v1;
  95 + v1[0] = T(x, y, z, 0, 1);
  96 + v1[1] = T(x, y, z, 1, 1);
  97 + v1[2] = T(x, y, z, 2, 1);
  98 +
  99 + e->v2[0] = T(x, y, z, 0, 2);
  100 + e->v2[1] = T(x, y, z, 1, 2);
  101 + e->v2[2] = T(x, y, z, 2, 2);
  102 +
  103 + e->lambda[0] = e->v0.len();
  104 + e->lambda[1] = v1.len();
  105 + e->lambda[2] = e->v2.len();
  106 + e->v0 = e->v0 / e->lambda[0];
  107 + e->v2 = e->v2 / e->lambda[2];
  108 +}
  109 +
  110 +/// Create a rotation matrix to orient a glyph along the tensor direction. Orientations are based on the input vectors
  111 +/// v0, v1, and v2
  112 +void get_glyph_rotation_matrix(float* R, stim::vec3<float> v0, stim::vec3<float> v2) {
  113 +
  114 + stim::matrix_sq<float, 4> M;
  115 + M(0, 0) = v0[0];
  116 + M(1, 0) = v0[1];
  117 + M(2, 0) = v0[2];
  118 + M(3, 0) = 0;
  119 +
  120 + M(0, 2) = v2[0];
  121 + M(1, 2) = v2[1];
  122 + M(2, 2) = v2[2];
  123 + M(3, 2) = 0;
  124 +
  125 + stim::vec3<float> a(v0[0], v0[1], v0[2]);
  126 + stim::vec3<float> b(v2[0], v2[1], v2[2]);
  127 + stim::vec3<float> c = a.cross(b);
  128 +
  129 + M(0, 1) = c[0];
  130 + M(1, 1) = c[1];
  131 + M(2, 1) = c[2];
  132 + M(3, 1) = 0;
  133 +
  134 + M(0, 3) = 0;
  135 + M(1, 3) = 0;
  136 + M(2, 3) = 0;
  137 + M(3, 3) = 1;
  138 +
  139 + memcpy(R, M.M, 16 * sizeof(float));
  140 + return;
  141 +}
  142 +
  143 +
  144 +// Perform the necessary updates when the user reshapes the window
  145 +void window_reshape() {
  146 + int width, height;
  147 +
  148 + //get the context size
  149 + glfwGetFramebufferSize(window, &width, &height);
  150 +
  151 +
  152 + //create an OpenGL viewport
  153 + glViewport(0, 0, width, height);
  154 +
  155 + //set the default orthographic view (assuming that the window and image aspect ratio are identical)
  156 + float left = 0.0f;
  157 + float right = (float)T.shape[2];
  158 + float bottom = 0.0f;
  159 + float top = (float)T.shape[1];
  160 +
  161 + //create an orthographic projection
  162 + glMatrixMode(GL_PROJECTION);
  163 + glLoadIdentity();
  164 + gluOrtho2D(left, right, bottom, top);
  165 +
  166 + glMatrixMode(GL_MODELVIEW);
  167 + glLoadIdentity();
  168 +}
  169 +
  170 +void init() {
  171 + glEnable(GL_DEPTH_TEST);
  172 + glEnable(GL_CULL_FACE);
  173 + glFrontFace(GL_CCW);
  174 +}
  175 +
  176 +void generate_glyph_points(int resolution = 5) {
  177 + int sectorCount = 2 * resolution;
  178 + int stackCount = resolution;
  179 +
  180 + int num_vertices = (sectorCount + 1) * (stackCount + 1);
  181 + sin_theta.resize(num_vertices);
  182 + sin_phi.resize(num_vertices);
  183 + cos_theta.resize(num_vertices);
  184 + cos_phi.resize(num_vertices);
  185 +
  186 + float theta, phi; // vertex texCoord
  187 +
  188 + float sectorStep = 2 * (float)PI / sectorCount;
  189 + float stackStep = (float)PI / stackCount;
  190 +
  191 + size_t i = 0;
  192 + //calculate the vertex positions
  193 + for (int phi_i = 0; phi_i <= stackCount; ++phi_i)
  194 + {
  195 +
  196 + // add (sectorCount+1) vertices per stack
  197 + // the first and last vertices have same position and normal, but different tex coords
  198 + for (int theta_i = 0; theta_i <= sectorCount; ++theta_i)
  199 + {
  200 +
  201 + // calculate the spherical coordinates of the vertex
  202 + theta = (float)theta_i / sectorCount * 2 * (float)PI;
  203 + phi = (float)phi_i / stackCount * (float)PI;
  204 +
  205 + //pre-compute the sine and cosine values used to create the superquadric
  206 + sin_theta[i] = sinf(theta);
  207 + sin_phi[i] = sinf(phi);
  208 + cos_theta[i] = cosf(theta);
  209 + cos_phi[i] = cosf(phi);
  210 + i++;
  211 + }
  212 + }
  213 +}
  214 +
  215 +inline float weird_exp(float x, float e) {
  216 + return sgn(x) * powf(fabs(x), e);
  217 +}
  218 +
  219 +stim::vec3<float> qx(size_t i, float alpha = 1, float beta = 1) {
  220 +
  221 + float sin_phi_beta = weird_exp(sin_phi[i], beta);
  222 + float sin_theta_alpha = weird_exp(sin_theta[i], alpha);
  223 + float cos_theta_alpha = weird_exp(cos_theta[i], alpha);
  224 + float cos_phi_beta = weird_exp(cos_phi[i], beta);
  225 +
  226 + stim::vec3<float> p;
  227 + p[0] = cos_phi_beta;
  228 + p[1] = -sin_theta_alpha * sin_phi_beta;
  229 + p[2] = cos_theta_alpha * sin_phi_beta;
  230 + return p;
  231 +}
  232 +
  233 +void render_triangle(stim::vec3<float> p0, stim::vec3<float> p1, stim::vec3<float> p2) {
  234 + stim::vec3<float> n = triangle_norm(p0, p1, p2);
  235 +
  236 + glNormal3f(n[0], n[1], n[2]);
  237 + glVertex3f(p0[0], p0[1], p0[2]);
  238 + glVertex3f(p1[0], p1[1], p1[2]);
  239 + glVertex3f(p2[0], p2[1], p2[2]);
  240 +}
  241 +
  242 +//render a glyph (0 = ellipsoid, 1 = superquadric)
  243 +void render_glyph(eigendecomposition* e, int glyph_type = 0, int resolution = 5) {
  244 + glMatrixMode(GL_MODELVIEW);
  245 + glPushMatrix();
  246 +
  247 + stim::vec3<float> norm_lambda = e->lambda.norm();
  248 + glScalef(glyph_scale, glyph_scale, glyph_scale);
  249 + glScalef(norm_lambda[2], norm_lambda[1], norm_lambda[0]);
  250 +
  251 + if (glyph_type == 0) {
  252 + GLUquadric* q = gluNewQuadric();
  253 + gluSphere(q, 1, 2 * resolution, resolution);
  254 + }
  255 + else if (glyph_type == 1) {
  256 + float radius = 0.5;
  257 + int sectorCount = 2 * resolution;
  258 + int stackCount = resolution;
  259 +
  260 + if (!glyph_calculated) {
  261 + generate_glyph_points(resolution);
  262 + glyph_calculated = true;
  263 + }
  264 +
  265 + //get the anisotropy values
  266 + float fa, cs, cl, cp, alpha, beta;
  267 + fa = get_anisotropy(e->lambda, cs, cl, cp);
  268 +
  269 + if (cl >= cp) {
  270 + alpha = powf(1 - cp, gamma);
  271 + beta = powf(1 - cl, gamma);
  272 + }
  273 + else {
  274 + alpha = powf(1 - cl, gamma);
  275 + beta = powf(1 - cp, gamma);
  276 + }
  277 +
  278 + //draw the sphere
  279 + int k1, k2;
  280 + float3 p[3];
  281 + float2 s[3];
  282 + float3 n;
  283 +
  284 + //draw the glyph
  285 + glBegin(GL_TRIANGLES);
  286 + for (int i = 0; i < stackCount; ++i)
  287 + {
  288 + k1 = i * (sectorCount + 1); // beginning of current stack
  289 + k2 = k1 + sectorCount + 1; // beginning of next stack
  290 +
  291 + for (int j = 0; j < sectorCount; ++j, ++k1, ++k2)
  292 + {
  293 + // 2 triangles per sector excluding first and last stacks
  294 + // k1 => k2 => k1+1
  295 + if (i != 0)
  296 + {
  297 +
  298 + stim::vec3<float> p0 = qx(k1, alpha, beta);
  299 + stim::vec3<float> p1 = qx(k2, alpha, beta);
  300 + stim::vec3<float> p2 = qx(k1 + 1, alpha, beta);
  301 +
  302 + render_triangle(p0, p1, p2);
  303 + }
  304 +
  305 + // k1+1 => k2 => k2+1
  306 + if (i != (stackCount - 1))
  307 + {
  308 + stim::vec3<float> p0 = qx(k1 + 1, alpha, beta);
  309 + stim::vec3<float> p1 = qx(k2, alpha, beta);
  310 + stim::vec3<float> p2 = qx(k2 + 1, alpha, beta);
  311 +
  312 + render_triangle(p0, p1, p2);
  313 + }
  314 + }
  315 + }
  316 + glEnd();
  317 + }
  318 + glPopMatrix();
  319 +}
  320 +
  321 +void lighting() {
  322 +
  323 + GLfloat light_position[] = { 1, 1, -1.0, 0.0 };
  324 + GLfloat light_ambient[] = { ambient_intensity, ambient_intensity, ambient_intensity, 1.0 };
  325 + GLfloat light_diffuse[] = { diffuse_intensity, diffuse_intensity, diffuse_intensity, 1.0 };
  326 + glShadeModel(GL_SMOOTH);
  327 +
  328 + glLightfv(GL_LIGHT0, GL_POSITION, light_position);
  329 + glLightfv(GL_LIGHT0, GL_DIFFUSE, light_diffuse);
  330 + glLightfv(GL_LIGHT0, GL_AMBIENT, light_ambient);
  331 +
  332 + glEnable(GL_LIGHTING);
  333 + glEnable(GL_LIGHT0);
  334 +
  335 + glColorMaterial(GL_FRONT, GL_AMBIENT_AND_DIFFUSE);
  336 + glEnable(GL_COLOR_MATERIAL);
  337 + glEnable(GL_NORMALIZE);
  338 +
  339 +}
  340 +void render() {
  341 + size_t zi = z_slice;
  342 +
  343 + float dc = 1.0f / T.shape[1];
  344 + float c[3];
  345 + eigendecomposition e;
  346 + float R[16];
  347 +
  348 +
  349 + glMatrixMode(GL_MODELVIEW_MATRIX);
  350 +
  351 + lighting();
  352 +
  353 + for (size_t yi = 0; yi < T.shape[1]; yi++) {
  354 + for (size_t xi = 0; xi < T.shape[2]; xi++) {
  355 + get_pixel_eigen(&e, xi, yi, zi);
  356 +
  357 + c[0] = abs(e.v0[0]);
  358 + c[1] = abs(e.v0[1]);
  359 + c[2] = abs(e.v0[2]);
  360 +
  361 + glPushMatrix();
  362 +
  363 +
  364 + glTranslatef((float)xi + 0.5f, (float)yi + 0.5f, 0);
  365 +
  366 + stim::vec3<float> v0(e.v0[0], e.v0[1], e.v0[2]);
  367 + stim::vec3<float> v2(e.v2[0], e.v2[1], e.v2[2]);
  368 + get_glyph_rotation_matrix(R, v0,v2);
  369 + glMultMatrixf((GLfloat*)R);
  370 +
  371 + glColor3f(c[0], c[1], c[2]);
  372 + render_glyph(&e, glyph_type, glyph_resolution);
  373 +
  374 +
  375 + glPopMatrix();
  376 + }
  377 + }
  378 + cout_frame = false;
  379 +
  380 +}
  381 +
  382 +void key_callback(GLFWwindow* window, int key, int scancode, int action, int mods)
  383 +{
  384 + if (key == GLFW_KEY_RIGHT && action != GLFW_RELEASE)
  385 + z_slice++;
  386 + else if (key == GLFW_KEY_LEFT && action != GLFW_RELEASE)
  387 + z_slice--;
  388 + if (z_slice >= T.shape[0]) z_slice = 0;
  389 + if (z_slice < 0) z_slice = T.shape[0] - 1;
  390 +}
  391 +
  392 +void display_rotation_matrix(eigendecomposition* e) {
  393 + float R[16];
  394 + stim::vec3<float> v0(e->v0[0], e->v0[1], e->v0[2]);
  395 + stim::vec3<float> v2(e->v2[0], e->v2[1], e->v2[2]);
  396 + get_glyph_rotation_matrix(R, v0, v2);
  397 + for (int r = 0; r < 4; r++) {
  398 + for (int c = 0; c < 4; c++) {
  399 + std::cout << R[r * 4 + c] << " ";
  400 + }
  401 + std::cout << std::endl;
  402 + }
  403 +}
  404 +
  405 +
  406 +void display_eigen(eigendecomposition* e) {
  407 + std::cout << "v0 = (" << e->v0[0] << ", " << e->v0[1] << ", " << e->v0[2] << ")" << std::endl;
  408 +}
  409 +
  410 +void mouse_button_callback(GLFWwindow* window, int button, int action, int mods) {
  411 + //if the user clicks inside the window, display information about the tensor field
  412 + if (button == GLFW_MOUSE_BUTTON_LEFT && action == GLFW_PRESS) {
  413 + double xpos, ypos;
  414 + glfwGetCursorPos(window, &xpos, &ypos); //get the position of the mouse pointer
  415 + int width, height;
  416 + glfwGetFramebufferSize(window, &width, &height);
  417 + int pixel_x = (int)(xpos / width * T.shape[2]);
  418 + int pixel_y = T.shape[1] - (int)(ypos / height * (int)T.shape[1]) - 1;
  419 + std::cout << "----------------------------" << std::endl;
  420 + std::cout << "x: " << pixel_x << " y: " << pixel_y << std::endl;
  421 + std::cout << "----------------------------" << std::endl;
  422 + eigendecomposition e;
  423 + get_pixel_eigen(&e, pixel_x, pixel_y, z_slice);
  424 + display_eigen(&e);
  425 + std::cout << "===========" << std::endl;
  426 + display_rotation_matrix(&e);
  427 + }
  428 +}
  429 +
  430 +void reshape_callback(GLFWwindow* window, int width, int height) {
  431 + window_reshape();
  432 + render();
  433 +}
  434 +
  435 +
  436 +int main(int argc, char** argv) {
  437 + std::string filename(argv[1]);
  438 + if (T.load_tira(filename) != TIRA_SUCCESS) {
  439 + return -1;
  440 + }
  441 +
  442 + /* Initialize the library */
  443 + if (!glfwInit())
  444 + return -1;
  445 +
  446 +
  447 + /* Create a windowed mode window and its OpenGL context */
  448 + window = glfwCreateWindow(700, 700, "Hello World", NULL, NULL);
  449 + if (!window)
  450 + {
  451 + glfwTerminate();
  452 + return -1;
  453 + }
  454 +
  455 + /* Make the window's context current */
  456 + glfwMakeContextCurrent(window);
  457 +
  458 + glfwSetKeyCallback(window, key_callback);
  459 + glfwSetFramebufferSizeCallback(window, reshape_callback);
  460 + glfwSetMouseButtonCallback(window, mouse_button_callback);
  461 +
  462 + GLenum err = glewInit();
  463 + //deal with a GLEW initialization failure
  464 + if (GLEW_OK != err)
  465 + std::cout << "GLEW Error: " << glewGetErrorString(err) << std::endl;
  466 +
  467 + //set up the window viewport
  468 + window_reshape();
  469 +
  470 + //initialize the OpenGL render details
  471 + init();
  472 +
  473 + /* Loop until the user closes the window */
  474 + while (!glfwWindowShouldClose(window))
  475 + {
  476 + /* Render here */
  477 + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
  478 +
  479 + render();
  480 +
  481 + /* Swap front and back buffers */
  482 + glfwSwapBuffers(window);
  483 +
  484 + /* Poll for and process events */
  485 + glfwPollEvents();
  486 + }
  487 +
  488 + glfwTerminate();
  489 + return 0;
  490 +}
0 491 \ No newline at end of file
... ...
tensorfield.cpp 0 → 100644
  Show/Hide comments   View file @e4667ce
  1 +++ a/tensorfield.cpp
... ...
tensorfield.h 0 → 100644
  Show/Hide comments   View file @e4667ce
  1 +++ a/tensorfield.h
  1 +#ifndef TIRA_TENSORFIELD
  2 +#define TIRA_TENSORFIELD
  3 +
  4 +#include <valarray>
  5 +#include <string>
  6 +#include <fstream>
  7 +#include <initializer_list>
  8 +
  9 +enum tira_floattype {unused0, unused1, float16, unused3, float32, unused4, unused5, unused6, unused7, float64};
  10 +enum tira_error {TIRA_SUCCESS, TIRA_UNDEFINED, TIRA_FILE_OPEN_ERROR, TIRA_UNKNOWN_FILE_TYPE};
  11 +enum tira_filetype {unused, spectrum, tensor2};
  12 +struct tira_header {
  13 + std::valarray<size_t> shape; //stores the resolution of the tensor field along each dimension
  14 + tira_filetype filetype; //stores the file type
  15 + tira_floattype datatype; //stores the data type of the tensor field
  16 + float version; //file version
  17 +};
  18 +
  19 +tira_error tira_load_header(std::ifstream &infile, tira_header &header) {
  20 +
  21 + //see if the file represents a TIRA structure
  22 + unsigned char TIRA[5];
  23 + TIRA[4] = 0;
  24 + infile.read((char*)TIRA, 4);
  25 + if (strcmp((char*)TIRA, "TIRA") != 0)
  26 + return TIRA_UNKNOWN_FILE_TYPE;
  27 +
  28 + //get the file type
  29 + infile.read((char*)&(header.filetype), sizeof(tira_filetype));
  30 +
  31 + //get the version
  32 + infile.read((char*)&(header.version), sizeof(float));
  33 +
  34 + //get the number of dimensions in the field
  35 + size_t ndim;
  36 + infile.read((char*)&ndim, sizeof(size_t));
  37 + header.shape.resize(ndim);
  38 +
  39 + //read each dimension
  40 + for (int d = 0; d < ndim; d++)
  41 + infile.read((char*)&(header.shape[d]), sizeof(size_t));
  42 +
  43 + //read the data type
  44 + infile.read((char*)&(header.datatype), sizeof(tira_floattype));
  45 +
  46 + return TIRA_SUCCESS;
  47 +}
  48 +
  49 +template <typename T>
  50 +class tensorfield {
  51 +public:
  52 +
  53 + std::valarray<size_t> shape; //number of samples in each dimension
  54 + std::valarray<size_t> cum_offsets; //cumulative factors for fast array accesses
  55 + tira_filetype filetype; //stores the file type
  56 + tira_floattype datatype; //stores the data type of the tensor field
  57 + T* ptr; //pointer to the raw tensor data
  58 +
  59 + void init() {
  60 + size_t ndims = shape.size(); //get the number of dimensions
  61 + cum_offsets.resize(ndims); //allocate space for the cumulative offsets
  62 + size_t cum = 1;
  63 + for (int i = 0; i < ndims; i++) {
  64 + cum_offsets[ndims - 1 - i] = cum;
  65 + cum *= shape[ndims - 1 - i];
  66 + }
  67 + }
  68 +
  69 + size_t bytes() {
  70 + size_t product = 1;
  71 + for (size_t i = 0; i < shape.size(); i++)
  72 + product *= shape[i];
  73 + return product * (size_t)datatype;
  74 + }
  75 +
  76 + // Function turns a series of array indices into a 1D index to the corresponding element
  77 + inline size_t idx(std::initializer_list<size_t> coords) {
  78 + size_t index = 0;
  79 + size_t i = 0;
  80 + for (auto c : coords) {
  81 + index += c * cum_offsets[i];
  82 + i++;
  83 + }
  84 + return index;
  85 + }
  86 +
  87 + tira_error load_tira(std::string filename) {
  88 +
  89 + //open the file for binary reading
  90 + std::ifstream infile(filename, std::ios::binary);
  91 + if (!infile) return TIRA_FILE_OPEN_ERROR;
  92 +
  93 + tira_header header;
  94 + tira_load_header(infile, header);
  95 + shape = header.shape;
  96 + filetype = header.filetype;
  97 + datatype = header.datatype;
  98 +
  99 + //allocate space in memory for the tensor field
  100 + size_t size_bytes = bytes();
  101 + ptr = (T*)malloc(size_bytes);
  102 +
  103 + //copy data from the file to memory
  104 + infile.read((char*)ptr, bytes());
  105 +
  106 + //close the file
  107 + infile.close();
  108 +
  109 + //initialize all internal variables
  110 + init();
  111 +
  112 + return TIRA_SUCCESS;
  113 + }
  114 +
  115 + T& operator()(size_t x, size_t y, size_t z, size_t r, size_t c) {
  116 + size_t i = idx({ z, y, x, r, c });
  117 + return ptr[i];
  118 + }
  119 +
  120 +};
  121 +
  122 +#endif
0 123 \ No newline at end of file
... ...