From c601924425187098e8481fa63618a2e69a3d80c0 Mon Sep 17 00:00:00 2001
From: bobarna <barnabas.borcsok@gmail.com>
Date: Thu, 13 May 2021 18:47:12 +0200
Subject: [PATCH] reworked whole codebase

---
 CMakeLists.txt                                |   60 +-
 data/plane.mtl                                |   12 +
 src/Constants.h                               |   10 +
 src/application.cpp                           |  124 --
 src/application.h                             |   68 -
 src/constants.h                               |    9 -
 .../geometry.cpp => geometries/Geometry.cpp}  |    2 +-
 src/geometries/Geometry.h                     |   27 +
 src/geometries/ObjGeometry.cpp                |   62 +
 src/geometries/ObjGeometry.h                  |   29 +
 .../PBDSimulation.cpp}                        |   67 +-
 .../PBDSimulation.h}                          |   27 +-
 src/geometries/ParamSurface.cpp               |   59 +
 src/geometries/ParamSurface.h                 |   42 +
 .../SPHSimulation.cpp}                        |  162 +-
 .../SPHSimulation.h}                          |   48 +-
 src/geometries/VertexData.cpp                 |   13 +
 src/geometries/VertexData.h                   |   18 +
 src/main.cpp                                  |  226 ++-
 src/objects/HeadObject.cpp                    |   23 +
 src/objects/HeadObject.h                      |   17 +
 .../objects/object.cpp => objects/Object.cpp} |   11 +-
 .../objects/object.h => objects/Object.h}     |   20 +-
 src/objects/PBDSimulationObject.cpp           |   22 +
 src/objects/PBDSimulationObject.h             |   24 +
 src/objects/Particle.h                        |   42 +
 src/objects/SPHSimulationObject.cpp           |   17 +
 src/objects/SPHSimulationObject.h             |   22 +
 src/rendering/Camera.cpp                      |   46 +
 src/rendering/Camera.h                        |   52 +
 src/rendering/RenderState.h                   |   18 +
 src/rendering/Scene.cpp                       |  113 ++
 src/rendering/Scene.h                         |   53 +
 src/rendering/camera.cpp                      |  101 -
 src/rendering/camera.h                        |   45 -
 src/rendering/geometry.h                      |   25 -
 src/rendering/lights/Light.cpp                |    5 +
 src/rendering/lights/Light.h                  |   13 +
 src/rendering/lights/light.cpp                |    6 -
 src/rendering/lights/light.h                  |   18 -
 src/rendering/materials/Material.cpp          |    1 +
 src/rendering/materials/Material.h            |   12 +
 src/rendering/materials/material.h            |   12 -
 src/rendering/objects/simulation_object.cpp   |   20 -
 src/rendering/objects/simulation_object.h     |   19 -
 src/rendering/render_state.h                  |   20 -
 .../{basic_shader.cpp => BasicShader.cpp}     |    8 +-
 .../shaders/{basic_shader.h => BasicShader.h} |    8 +-
 src/rendering/shaders/PhongShader.cpp         |   20 +
 src/rendering/shaders/PhongShader.h           |   99 +
 .../shaders/{shader.cpp => Shader.cpp}        |    4 +-
 src/rendering/shaders/Shader.h                |   25 +
 src/rendering/shaders/gpu_program.cpp         |    1 -
 src/rendering/shaders/shader.h                |   19 -
 src/rendering/textures/Texture.cpp            |    1 +
 .../textures/{texture.h => Texture.h}         |   36 +-
 src/rendering/textures/texture.cpp            |    5 -
 .../textures/texture_uniform_color.h          |   17 -
 src/rendering/vertex_data.cpp                 |    8 -
 src/rendering/vertex_data.h                   |   17 -
 src/scene.cpp                                 |   35 -
 src/scene.h                                   |   32 -
 src/simulation/particle.cpp                   |   25 -
 src/simulation/particle.h                     |   35 -
 src/simulation/simulation.h                   |   19 -
 src/utility/custom_math.h                     |   69 -
 src/utility/gl.h                              |    7 -
 .../gpu_program.h => utils/GPUProgram.h}      |   18 +-
 .../InputHandler.cpp}                         |   12 +-
 .../input_handler.h => utils/InputHandler.h}  |   34 +-
 src/utils/OBJReader.cpp                       |   89 +
 src/utils/OBJReader.h                         |   46 +
 src/utils/math.cpp                            |   20 +
 src/utils/math.h                              |  378 ++++
 src/utils/other/stb_image_write.h             | 1692 +++++++++++++++++
 src/utils/save_image.h                        |   30 +
 src/utils/util.cpp                            |   25 +
 src/utils/util.h                              |   16 +
 78 files changed, 3705 insertions(+), 987 deletions(-)
 create mode 100644 data/plane.mtl
 create mode 100644 src/Constants.h
 delete mode 100644 src/application.cpp
 delete mode 100644 src/application.h
 delete mode 100644 src/constants.h
 rename src/{rendering/geometry.cpp => geometries/Geometry.cpp} (91%)
 create mode 100644 src/geometries/Geometry.h
 create mode 100644 src/geometries/ObjGeometry.cpp
 create mode 100644 src/geometries/ObjGeometry.h
 rename src/{simulation/pbd/pbd_simulation.cpp => geometries/PBDSimulation.cpp} (70%)
 rename src/{simulation/pbd/pbd_simulation.h => geometries/PBDSimulation.h} (73%)
 create mode 100644 src/geometries/ParamSurface.cpp
 create mode 100644 src/geometries/ParamSurface.h
 rename src/{simulation/sph/sph_simulation.cpp => geometries/SPHSimulation.cpp} (54%)
 rename src/{simulation/sph/sph_simulation.h => geometries/SPHSimulation.h} (55%)
 create mode 100644 src/geometries/VertexData.cpp
 create mode 100644 src/geometries/VertexData.h
 create mode 100644 src/objects/HeadObject.cpp
 create mode 100644 src/objects/HeadObject.h
 rename src/{rendering/objects/object.cpp => objects/Object.cpp} (89%)
 rename src/{rendering/objects/object.h => objects/Object.h} (59%)
 create mode 100644 src/objects/PBDSimulationObject.cpp
 create mode 100644 src/objects/PBDSimulationObject.h
 create mode 100644 src/objects/Particle.h
 create mode 100644 src/objects/SPHSimulationObject.cpp
 create mode 100644 src/objects/SPHSimulationObject.h
 create mode 100644 src/rendering/Camera.cpp
 create mode 100644 src/rendering/Camera.h
 create mode 100644 src/rendering/RenderState.h
 create mode 100644 src/rendering/Scene.cpp
 create mode 100644 src/rendering/Scene.h
 delete mode 100644 src/rendering/camera.cpp
 delete mode 100644 src/rendering/camera.h
 delete mode 100644 src/rendering/geometry.h
 create mode 100644 src/rendering/lights/Light.cpp
 create mode 100644 src/rendering/lights/Light.h
 delete mode 100644 src/rendering/lights/light.cpp
 delete mode 100644 src/rendering/lights/light.h
 create mode 100644 src/rendering/materials/Material.cpp
 create mode 100644 src/rendering/materials/Material.h
 delete mode 100644 src/rendering/materials/material.h
 delete mode 100644 src/rendering/objects/simulation_object.cpp
 delete mode 100644 src/rendering/objects/simulation_object.h
 delete mode 100644 src/rendering/render_state.h
 rename src/rendering/shaders/{basic_shader.cpp => BasicShader.cpp} (73%)
 rename src/rendering/shaders/{basic_shader.h => BasicShader.h} (90%)
 create mode 100644 src/rendering/shaders/PhongShader.cpp
 create mode 100644 src/rendering/shaders/PhongShader.h
 rename src/rendering/shaders/{shader.cpp => Shader.cpp} (95%)
 create mode 100644 src/rendering/shaders/Shader.h
 delete mode 100644 src/rendering/shaders/gpu_program.cpp
 delete mode 100644 src/rendering/shaders/shader.h
 create mode 100644 src/rendering/textures/Texture.cpp
 rename src/rendering/textures/{texture.h => Texture.h} (73%)
 delete mode 100644 src/rendering/textures/texture.cpp
 delete mode 100644 src/rendering/textures/texture_uniform_color.h
 delete mode 100644 src/rendering/vertex_data.cpp
 delete mode 100644 src/rendering/vertex_data.h
 delete mode 100644 src/scene.cpp
 delete mode 100644 src/scene.h
 delete mode 100644 src/simulation/particle.cpp
 delete mode 100644 src/simulation/particle.h
 delete mode 100644 src/simulation/simulation.h
 delete mode 100644 src/utility/custom_math.h
 delete mode 100644 src/utility/gl.h
 rename src/{rendering/shaders/gpu_program.h => utils/GPUProgram.h} (96%)
 rename src/{utility/input_handler.cpp => utils/InputHandler.cpp} (82%)
 rename src/{utility/input_handler.h => utils/InputHandler.h} (63%)
 create mode 100644 src/utils/OBJReader.cpp
 create mode 100644 src/utils/OBJReader.h
 create mode 100644 src/utils/math.cpp
 create mode 100644 src/utils/math.h
 create mode 100644 src/utils/other/stb_image_write.h
 create mode 100644 src/utils/save_image.h
 create mode 100644 src/utils/util.cpp
 create mode 100644 src/utils/util.h

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 26be728..eb6206e 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,47 +1,47 @@
 cmake_minimum_required(VERSION 3.1)
+
 project(madid)
 
 set(CMAKE_CXX_STANDARD 17)
 #set(CMAKE_CXX_FLAGS "-Wall -Wextra -Werror -pedantic")
-set(CMAKE_CXX_FLAGS "-Wall -Wextra -pedantic -Wunused-parameter")
-set(OpenGL_GL_PREFERENCE "LEGACY")
+set(CMAKE_CXX_FLAGS "-Wall -Wextra -pedantic")
 
 add_executable(${PROJECT_NAME}
         src/main.cpp
-        src/constants.h 
-        src/utility/gl.h
-        src/application.cpp
-        src/application.h
-        src/scene.cpp
-        src/scene.h
-        src/rendering/camera.cpp
-        src/rendering/camera.h
-        src/utility/input_handler.h
-        src/utility/input_handler.cpp
-
-        src/simulation/particle.cpp
-        src/simulation/particle.h
-        src/simulation/pbd/pbd_simulation.cpp
-        src/simulation/pbd/pbd_simulation.h
-        src/simulation/sph/sph_simulation.cpp
-        src/simulation/sph/sph_simulation.h
-
-        src/rendering/geometry.cpp src/rendering/geometry.h
-        src/rendering/vertex_data.cpp src/rendering/vertex_data.h
-        src/rendering/objects/object.cpp src/rendering/objects/object.h src/rendering/objects/simulation_object.cpp src/rendering/objects/simulation_object.h src/simulation/simulation.h src/utility/custom_math.h src/rendering/render_state.h src/rendering/textures/texture.cpp src/rendering/textures/texture.h src/rendering/shaders/shader.cpp src/rendering/shaders/shader.h src/rendering/shaders/gpu_program.cpp src/rendering/shaders/gpu_program.h src/rendering/textures/texture_uniform_color.h src/rendering/materials/material.h src/rendering/lights/light.cpp src/rendering/lights/light.h src/rendering/shaders/basic_shader.cpp src/rendering/shaders/basic_shader.h)
-
+        src/utils/math.h src/utils/math.cpp
+        src/utils/util.h src/utils/util.cpp
+        src/geometries/PBDSimulation.h src/geometries/PBDSimulation.cpp
+        src/rendering/Camera.h src/rendering/Camera.cpp
+        src/objects/Object.h src/objects/Object.cpp
+        src/rendering/Scene.cpp
+        src/rendering/shaders/Shader.cpp
+        src/utils/InputHandler.cpp
+        src/geometries/Geometry.cpp
+        src/objects/PBDSimulationObject.cpp
+        src/objects/HeadObject.cpp
+        src/rendering/shaders/BasicShader.cpp src/rendering/shaders/BasicShader.h
+        src/rendering/shaders/PhongShader.cpp src/rendering/shaders/PhongShader.h
+        src/geometries/ParamSurface.cpp src/geometries/ParamSurface.h
+        src/geometries/VertexData.h src/rendering/materials/Material.cpp src/rendering/materials/Material.h src/rendering/textures/Texture.cpp src/rendering/textures/Texture.h src/rendering/lights/Light.cpp src/rendering/lights/Light.h src/geometries/VertexData.cpp src/utils/OBJReader.cpp src/utils/OBJReader.h src/geometries/ObjGeometry.cpp src/geometries/ObjGeometry.h
+        src/geometries/SPHSimulation.cpp src/geometries/SPHSimulation.h
+        src/objects/SPHSimulationObject.cpp src/objects/SPHSimulationObject.h)
 
 # opengl
 find_package(OpenGL REQUIRED)
+include_directories(${OPENGL_INCLUDE_DIRS})
+
+# glfw
+find_package(glfw3 REQUIRED)
+include_directories(${GLFW_INCLUDE_DIRS})
+link_libraries(${GLFW_LIBRARY_DIRS})
 
-link_directories(${GTKMM_LIBRARY_DIRS})
-include_directories(${GTKMM_INCLUDE_DIRS})
+# glew
+find_package(GLEW REQUIRED)
+include_directories(${GLEW_INCLUDE_DIRS})
 
 target_link_libraries(
         ${PROJECT_NAME}
-        GL
-        GLU
         glfw
-        GLEW
+        ${OPENGL_LIBRARIES}
+        ${GLEW_LIBRARIES}
 )
-
diff --git a/data/plane.mtl b/data/plane.mtl
new file mode 100644
index 0000000..e8374a5
--- /dev/null
+++ b/data/plane.mtl
@@ -0,0 +1,12 @@
+# Blender MTL File: 'None'
+# Material Count: 1
+
+newmtl Material
+Ns 323.999994
+Ka 1.000000 1.000000 1.000000
+Kd 0.800000 0.800000 0.800000
+Ks 0.500000 0.500000 0.500000
+Ke 0.000000 0.000000 0.000000
+Ni 1.450000
+d 1.000000
+illum 2
diff --git a/src/Constants.h b/src/Constants.h
new file mode 100644
index 0000000..c9ebce7
--- /dev/null
+++ b/src/Constants.h
@@ -0,0 +1,10 @@
+#ifndef BRAVE2_CONSTANTS_H
+#define BRAVE2_CONSTANTS_H
+
+static int WIDTH = 600;
+static int HEIGHT = 400;
+
+static float GRAVITY_ABS_VALUE = 0.98f;
+
+
+#endif //BRAVE2_CONSTANTS_H
diff --git a/src/application.cpp b/src/application.cpp
deleted file mode 100644
index e1dad28..0000000
--- a/src/application.cpp
+++ /dev/null
@@ -1,124 +0,0 @@
-#include "application.h"
-
-InputHandler *Application::inputHandler = InputHandler::GetInstance();
-bool Application::dragging = false;
-
-Application::Application() {
-    window = glfwCreateWindow(WIDTH, HEIGHT, NAME, NULL, NULL);
-}
-
-Application::~Application() {
-    glfwDestroyWindow(window);
-    glfwTerminate();
-}
-
-GLFWwindow *Application::get_window() {
-    return window;
-}
-
-void Application::init() {
-    this->set_callbacks();
-
-    glViewport(0, 0, this->width, this->height);
-    glMatrixMode(GL_MODELVIEW);
-    glEnable(GL_DEPTH_TEST);
-    glDepthMask(GL_TRUE);
-    glLoadIdentity();
-    glClearColor(255.0f / 255.0f, 25.0f / 255.0f, 25.0f / 255.0f, 1.0f);
-}
-
-void Application::error_callback(int error, const char *description) {
-    std::cerr << "Error: " << description << "." << std::endl;
-}
-
-void Application::resize_callback(GLFWwindow *window, int w, int h) {
-    // TODO set application's width and height
-    if (h < 1)
-        h = 1;
-    glViewport(0, 0, w, h);
-    glMatrixMode(GL_PROJECTION);
-    glLoadIdentity();
-    gluPerspective(45.0f, (float) w / (float) h, 0.1f, 1000.0f);
-    gluLookAt(0.0f, 0.0f, 30, 0.0f, 0.0f, 0.0f, 0.0f, 1.0f, 0.0f);
-    glMatrixMode(GL_MODELVIEW);
-}
-
-void Application::key_callback(GLFWwindow *window,
-                               int key,
-                               int scancode,
-                               int action,
-                               int mods) {
-//    std::cout << key << " was pressed" << std::endl;
-    if (key == GLFW_KEY_ESCAPE)
-        glfwSetWindowShouldClose(window, GLFW_TRUE);
-    if (action == GLFW_RELEASE)
-        InputHandler::KeyRelease(key);
-    else
-        InputHandler::KeyPress(key);
-
-    InputHandler::SetModifiers(mods);
-}
-
-
-void Application::mouse_click_callback(GLFWwindow *window,
-                                       int button,
-                                       int action,
-                                       int mods) {
-    switch (button) {
-        case GLFW_MOUSE_BUTTON_1:
-            dragging = (action == GLFW_PRESS);
-            break;
-    }
-}
-
-
-void Application::mouse_motion_callback(GLFWwindow *window,
-                                        double x,
-                                        double y) {
-    printf("mouse");
-}
-
-void Application::set_callbacks() {
-    glfwSetErrorCallback(error_callback);
-    glfwSetWindowSizeCallback(window, resize_callback);
-    glfwSetKeyCallback(window, key_callback);
-    glfwSetMouseButtonCallback(window, mouse_click_callback);
-    glfwSetCursorPosCallback(window, mouse_motion_callback);
-}
-
-void Application::start() {
-    while (!glfwWindowShouldClose(this->window)) {
-        this->run();
-        glfwSwapBuffers(window);
-        glfwPollEvents();
-    }
-    glfwTerminate();
-}
-
-void Application::run() {
-    // ticking every 24 FPS
-    bool tick = false;
-
-    double delta_time = glfwGetTime() - time_since_last_frame;
-    time_since_last_frame = glfwGetTime();
-
-    if (time_since_last_frame >= 1.0 / 24.0) {
-        glfwSetTime(0.0f);
-        time_since_last_frame = 0.0f;
-        tick = true;
-    }
-
-    // TODO think about rendering only in 24 FPS
-    // TODO guarantee delta_time to be infinitesimal
-    scene->update(delta_time);
-
-
-    //clear color and depth buffer
-    glClearColor(0, 0, 0, 1);
-    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
-    glLoadIdentity(); //load identity matrix
-
-    scene->render();
-}
-
-
diff --git a/src/application.h b/src/application.h
deleted file mode 100644
index fcb60d0..0000000
--- a/src/application.h
+++ /dev/null
@@ -1,68 +0,0 @@
-#ifndef MADID_APPLICATION_H
-#define MADID_APPLICATION_H
-
-#include <GL/glew.h>
-#include <GLFW/glfw3.h>
-#include <glm/glm.hpp>
-
-#include <iostream>
-#include <memory>
-
-#include "constants.h"
-#include "scene.h"
-#include "rendering/camera.h"
-#include "utility/input_handler.h"
-
-// TODO make this class a singleton
-class Application {
-    int width = WIDTH;
-    int height = HEIGHT;
-
-    static bool dragging;
-
-    // Callback functions
-    static void error_callback(int error, const char *description);
-
-    static void resize_callback(GLFWwindow *window, int w, int h);
-
-    static void key_callback(GLFWwindow *window, int key, int scancode,
-                             int action, int mods);
-
-    static void mouse_click_callback(GLFWwindow *window, int button, int action,
-                                     int mods);
-
-    static void mouse_motion_callback(GLFWwindow *window, double x, double y);
-
-    std::shared_ptr<Scene> scene = std::make_shared<Scene>();
-
-    static InputHandler *inputHandler;
-
-    void run();
-
-    double time_since_last_frame;
-
-    // TODO
-    //boolean flag, indicates whether a screen capture event should be performed
-    //on the next display event.
-    /* bool screen_capture; */
-    /* bool zoom_mode; */
-    /* bool pan_mode; */
-    /* bool trackball_mode; */
-public:
-    Application();
-
-    ~Application();
-
-    GLFWwindow *window;
-
-    GLFWwindow *get_window();
-
-    void init();
-
-    void set_callbacks();
-
-    void start();
-};
-
-
-#endif //MADID_APPLICATION_H
diff --git a/src/constants.h b/src/constants.h
deleted file mode 100644
index 21a5693..0000000
--- a/src/constants.h
+++ /dev/null
@@ -1,9 +0,0 @@
-#ifndef MADID_CONSTANTS_H
-#define MADID_CONSTANTS_H
-
-const int WIDTH = 500;
-const int HEIGHT = 500;
-
-char *const NAME = "Madid";
-
-#endif
diff --git a/src/rendering/geometry.cpp b/src/geometries/Geometry.cpp
similarity index 91%
rename from src/rendering/geometry.cpp
rename to src/geometries/Geometry.cpp
index eaf6ee8..302b06d 100644
--- a/src/rendering/geometry.cpp
+++ b/src/geometries/Geometry.cpp
@@ -1,4 +1,4 @@
-#include "geometry.h"
+#include "Geometry.h"
 
 Geometry::Geometry() {
     glGenVertexArrays(1, &vao);
diff --git a/src/geometries/Geometry.h b/src/geometries/Geometry.h
new file mode 100644
index 0000000..96f5b3f
--- /dev/null
+++ b/src/geometries/Geometry.h
@@ -0,0 +1,27 @@
+#ifndef BRAVE2_GEOMETRY_H
+#define BRAVE2_GEOMETRY_H
+
+#include "../utils/math.h"
+#include <GL/glew.h>
+#include "VertexData.h"
+
+class Geometry {
+protected:
+    GLuint vao, vbo;
+
+public:
+    std::vector<VertexData> vtxData;    // vertices on the CPU
+
+    Geometry();
+
+    virtual void Draw() = 0;
+
+    virtual VertexData GetVertexDataByUV(float u, float v) {
+        //TODO fix this horrible quick stuff
+        return vtxData[std::rand()%vtxData.size()];
+    };
+
+    ~Geometry();
+};
+
+#endif //BRAVE2_GEOMETRY_H
\ No newline at end of file
diff --git a/src/geometries/ObjGeometry.cpp b/src/geometries/ObjGeometry.cpp
new file mode 100644
index 0000000..ba4976b
--- /dev/null
+++ b/src/geometries/ObjGeometry.cpp
@@ -0,0 +1,62 @@
+#include "ObjGeometry.h"
+
+ObjGeometry::ObjGeometry(const std::string &objPath) :
+        Geometry(),
+        objReader(objPath, vtxData) {
+    Initialize();
+}
+
+void ObjGeometry::Initialize() {
+    glBufferData(GL_ARRAY_BUFFER, sizeof(VertexData) * vtxData.size(), &vtxData[0], GL_STATIC_DRAW);
+    // Enable the vertex attribute arrays
+    glEnableVertexAttribArray(0);  // attribute array 0 = POSITION
+    glEnableVertexAttribArray(1);  // attribute array 1 = NORMAL
+    glEnableVertexAttribArray(2);  // attribute array 2 = TEXCOORD
+    // attribute array, components/attribute, component type, normalize?, stride, offset
+    glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, sizeof(VertexData), (void *) offsetof(VertexData, position));
+    glVertexAttribPointer(1, 3, GL_FLOAT, GL_FALSE, sizeof(VertexData), (void *) offsetof(VertexData, normal));
+    glVertexAttribPointer(2, 2, GL_FLOAT, GL_FALSE, sizeof(VertexData), (void *) offsetof(VertexData, texcoord));
+}
+
+void ObjGeometry::Draw() {
+    glBindVertexArray(vao);
+    glDrawArrays(GL_TRIANGLES, 0, vtxData.size());
+}
+
+VertexData ObjGeometry::GetVertexDataByUV(float u, float v) {
+
+    auto calc_a_b_c = [](vec2 p, vec2 t1, vec2 t2, vec2 t3) {
+        float denom = (t2.y - t3.y) * (t1.x - t3.x) + (t3.x - t2.x) * (t1.y - t3.y);
+        float a = ((t2.y - t3.y) * (p.x - t3.x) + (t3.x - t2.x) * (p.y - t3.y)) / denom;
+        float b = ((t3.y - t1.y) * (p.x - t3.y) + (t1.x - t3.x) * (p.y - t3.y)) / denom;
+
+        return vec3(a, b, 1 - a - b);
+    };
+
+    // iterates through all of the triangles of the mesh
+    for (size_t i = 0; i < vtxData.size(); i += 3) {
+        VertexData t1 = vtxData.at(i);
+        VertexData t2 = vtxData.at(i + 1);
+        VertexData t3 = vtxData.at(i + 2);
+
+        // check if point is inside the triangle
+        // this could be in a function, but I'm doing premature optimization - cause of all evil...
+        // using baricentric method:
+        // p = a*t1 + b*t2 + c*t3
+        // p is inside triangle (t1, t2, t3) if and only if ( 0 <= a,b,c <= 1 )
+        vec3 a_b_c{calc_a_b_c({u, v}, t1.texcoord, t2.texcoord, t3.texcoord)};
+
+        // if point is inside the triangle
+        if (0 <= a_b_c.x && a_b_c.x <= 1 &&
+            0 <= a_b_c.y && a_b_c.y <= 1 &&
+            0 <= a_b_c.z && a_b_c.z <= 1) {
+            return {
+                    a_b_c.x * t1.position + a_b_c.y * t2.position + a_b_c.z * t3.position,
+                    a_b_c.x * t1.normal + a_b_c.y * t2.normal + a_b_c.z * t3.normal,
+                    {u, v}
+            };
+        }
+    }
+
+    return {(VertexData) false};
+}
diff --git a/src/geometries/ObjGeometry.h b/src/geometries/ObjGeometry.h
new file mode 100644
index 0000000..a9bffd2
--- /dev/null
+++ b/src/geometries/ObjGeometry.h
@@ -0,0 +1,29 @@
+#ifndef BRAVE2_OBJGEOMETRY_H
+#define BRAVE2_OBJGEOMETRY_H
+
+
+#include "../utils/math.h"
+#include "../objects/Object.h"
+#include "../utils/OBJReader.h"
+
+
+// This class assumes that the mesh is made up entirely of triangles
+// this is important for the GetVertexDataByUV function
+class ObjGeometry : public Geometry {
+    std::string objPath;
+    OBJReader objReader;
+
+public:
+    //TODO read every data from OBJ file
+    ObjGeometry(const std::string &objPath);
+
+    void Initialize();
+
+    void Draw() override;
+
+    // Assumes that the mesh is made up of triangles
+    VertexData GetVertexDataByUV(float u, float v) override;
+};
+
+
+#endif //BRAVE2_OBJGEOMETRY_H
diff --git a/src/simulation/pbd/pbd_simulation.cpp b/src/geometries/PBDSimulation.cpp
similarity index 70%
rename from src/simulation/pbd/pbd_simulation.cpp
rename to src/geometries/PBDSimulation.cpp
index af51ec3..195b941 100644
--- a/src/simulation/pbd/pbd_simulation.cpp
+++ b/src/geometries/PBDSimulation.cpp
@@ -1,4 +1,5 @@
-#include "pbd_simulation.h"
+#include <random>
+#include "PBDSimulation.h"
 
 void PBDSimulation::addForce(vec3 force) {
     externalForces += force;
@@ -12,7 +13,7 @@ PBDSimulation::PBDSimulation(size_t _nr_sims, size_t _nr_segments, float _l_seg)
 
     // placing the fibers of the cloth
     for (size_t i = 0; i < nrStrands; i++) {
-        vec3 currPos(-.5f, 0.f, lSeg * (float)i);
+        vec3 currPos(-0.15f, -0.15f, lSeg*0.5f * (float)i - 0.15);
 
         vec3 color = vec3(0.9f, 0.3f, 0.3f);
         strands.emplace_back(CreateFiber(nrSegments, lSeg, currPos * vec3(1, -1, 1), color));
@@ -80,8 +81,8 @@ void PBDSimulation::solve_distance_constraint(Particle *p1, Particle *p2, float
                 (length(p1->tmp_pos - p2->tmp_pos) - dist) *
                 (p1->tmp_pos - p2->tmp_pos) / length(p1->tmp_pos - p2->tmp_pos);
 
-    p1->tmp_pos += 0.8f*d_p1;
-    p2->tmp_pos += 0.8f*d_p2;
+    p1->tmp_pos += d_p1;
+    p2->tmp_pos += d_p2;
 }
 
 void PBDSimulation::solve_bending_constraint(Particle *p1, Particle *p2, float dist) {
@@ -92,8 +93,8 @@ void PBDSimulation::solve_bending_constraint(Particle *p1, Particle *p2, float d
                 (length(p1->tmp_pos - p2->tmp_pos) - dist) *
                 (p1->tmp_pos - p2->tmp_pos) / length(p1->tmp_pos - p2->tmp_pos);
 
-    p1->tmp_pos += 0.6f*d_p1;
-    p2->tmp_pos += 0.6f*d_p2;
+    p1->tmp_pos += 0.6 * d_p1;
+    p2->tmp_pos += 0.6 * d_p2;
 }
 
 void PBDSimulation::solve_collision_constraint(Particle *p, vec3 &q1, vec3 &q2, vec3 &q3) {
@@ -106,9 +107,11 @@ void PBDSimulation::solve_collision_constraint(Particle *p, vec3 &q1, vec3 &q2,
     p->tmp_pos += d_p;
 }
 
-void PBDSimulation::render() {
+void PBDSimulation::Draw() {
     std::vector<float> particlePosAndColor;
-    for (auto &strand : strands)
+    for (size_t curr_strand = 0; curr_strand < strands.size()-1; curr_strand ++) {
+        auto strand = strands[curr_strand];
+        auto next_strand = strands[curr_strand+1];
         for (size_t i = 0; i < strand.size() - 1; ++i) {
             particlePosAndColor.push_back(strand[i]->pos.x);
             particlePosAndColor.push_back(strand[i]->pos.y);
@@ -123,7 +126,37 @@ void PBDSimulation::render() {
             particlePosAndColor.push_back(strand[i + 1]->color.x);
             particlePosAndColor.push_back(strand[i + 1]->color.y);
             particlePosAndColor.push_back(strand[i + 1]->color.z);
+
+            particlePosAndColor.push_back(next_strand[i]->pos.x);
+            particlePosAndColor.push_back(next_strand[i]->pos.y);
+            particlePosAndColor.push_back(next_strand[i]->pos.z);
+            particlePosAndColor.push_back(next_strand[i]->color.x);
+            particlePosAndColor.push_back(next_strand[i]->color.y);
+            particlePosAndColor.push_back(next_strand[i]->color.z);
+
+            particlePosAndColor.push_back(next_strand[i]->pos.x);
+            particlePosAndColor.push_back(next_strand[i]->pos.y);
+            particlePosAndColor.push_back(next_strand[i]->pos.z);
+            particlePosAndColor.push_back(next_strand[i]->color.x);
+            particlePosAndColor.push_back(next_strand[i]->color.y);
+            particlePosAndColor.push_back(next_strand[i]->color.z);
+
+            particlePosAndColor.push_back(strand[i + 1]->pos.x);
+            particlePosAndColor.push_back(strand[i + 1]->pos.y);
+            particlePosAndColor.push_back(strand[i + 1]->pos.z);
+            particlePosAndColor.push_back(strand[i + 1]->color.x);
+            particlePosAndColor.push_back(strand[i + 1]->color.y);
+            particlePosAndColor.push_back(strand[i + 1]->color.z);
+
+            particlePosAndColor.push_back(next_strand[i+1]->pos.x);
+            particlePosAndColor.push_back(next_strand[i+1]->pos.y);
+            particlePosAndColor.push_back(next_strand[i+1]->pos.z);
+            particlePosAndColor.push_back(next_strand[i+1]->color.x);
+            particlePosAndColor.push_back(next_strand[i+1]->color.y);
+            particlePosAndColor.push_back(next_strand[i+1]->color.z);
         }
+    }
+
 
     glBindVertexArray(vao);
     glBindBuffer(GL_ARRAY_BUFFER, vbo);
@@ -134,14 +167,18 @@ void PBDSimulation::render() {
 
     // position attribute
     glEnableVertexAttribArray(0);
-    glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, sizeof(float) * 6, (void *) 0);
-
-    // color attribute
+    // normal attribute
     glEnableVertexAttribArray(1);
+
+    glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, sizeof(float) * 6, (void *) 0);
     glVertexAttribPointer(1, 3, GL_FLOAT, GL_FALSE, sizeof(float) * 6, (void *) (sizeof(float) * 3));
 
     glLineWidth(1.0f);
-    glDrawArrays(GL_LINES, 0, particlePosAndColor.size() / 6);
+    glDrawArrays(GL_TRIANGLES, 0, particlePosAndColor.size() / 6);
+}
+
+vec3 PBDSimulation::getExternalForces() const {
+    return externalForces;
 }
 
 std::vector<Particle *> PBDSimulation::CreateFiber(size_t n, float l, vec3 startPos, vec3 color) {
@@ -149,8 +186,10 @@ std::vector<Particle *> PBDSimulation::CreateFiber(size_t n, float l, vec3 start
     std::vector<Particle *> currentStrand;
 
     for (size_t i = 0; i < n; i++) {
-        // slightly random mass for better visual results
-        float m = (float) rand() / RAND_MAX * .35f + .15f;
+        // random mass
+        // A single strand of hair can weigh between .2 – .5 milligrams,
+        // which is 2.0e-7 kg == 10^(-7) kg
+        float m = util::randomOffsetf(.35f, .15f);
 
         // first and last particle's position is infinite
         if (i == 0 || i == n-1) currentStrand.push_back(new Particle(currPos, 0, color));
diff --git a/src/simulation/pbd/pbd_simulation.h b/src/geometries/PBDSimulation.h
similarity index 73%
rename from src/simulation/pbd/pbd_simulation.h
rename to src/geometries/PBDSimulation.h
index 90766cc..e5906e4 100644
--- a/src/simulation/pbd/pbd_simulation.h
+++ b/src/geometries/PBDSimulation.h
@@ -1,16 +1,11 @@
-#ifndef MADID_PBD_SIMULATION_H
-#define MADID_PBD_SIMULATION_H
+#ifndef BRAVE2_PBDSIMULATION_H
+#define BRAVE2_PBDSIMULATION_H
 
-#include <cstddef>
-#include <vector>
-#include <glm/glm.hpp>
-
-#include "../../utility/gl.h"
-#include "../particle.h"
-#include "../../rendering/geometry.h"
-#include "../simulation.h"
-
-typedef glm::vec3 vec3;
+#include "../utils/math.h"
+#include "../utils/util.h"
+#include "../objects/Particle.h"
+#include "Geometry.h"
+#include "../objects/HeadObject.h"
 
 class PBDSimulation : public Geometry {
     void solve_distance_constraint(Particle *p1, Particle *p2, float dist);
@@ -42,9 +37,13 @@ public:
 
     PBDSimulation(size_t _nr_sims, size_t _nr_segments, float _l_seg);
 
+    void propagateHead();
+
     void update(float dt);
 
-    void render() override;
+    void Draw();
+
+    vec3 getExternalForces() const;
 
     void resetExternalForces();
 
@@ -52,4 +51,4 @@ public:
 };
 
 
-#endif //MADID_PBD_SIMULATION_H
+#endif //BRAVE2_PBDSIMULATION_H
diff --git a/src/geometries/ParamSurface.cpp b/src/geometries/ParamSurface.cpp
new file mode 100644
index 0000000..c1126af
--- /dev/null
+++ b/src/geometries/ParamSurface.cpp
@@ -0,0 +1,59 @@
+#include "ParamSurface.h"
+
+ParamSurface::ParamSurface() {
+    nVtxPerStrip = nStrips = 0;
+}
+
+VertexData ParamSurface::GenVertexData(float u, float v) {
+    VertexData vtxData;
+    vtxData.texcoord = vec2(u, v);
+    Dnum2 X, Y, Z;
+    Dnum2 U(u, vec2(1, 0)), V(v, vec2(0, 1));
+    eval(U, V, X, Y, Z);
+    vtxData.position = vec3(X.f, Y.f, Z.f);
+    vec3 drdU(X.d.x, Y.d.x, Z.d.x), drdV(X.d.y, Y.d.y, Z.d.y);
+    vtxData.normal = cross(drdU, drdV);
+    return vtxData;
+}
+
+void ParamSurface::create(int N, int M) {
+    nVtxPerStrip = (M + 1) * 2;
+    nStrips = N;
+    for (int i = 0; i < N; i++) {
+        for (int j = 0; j <= M; j++) {
+            vtxData.push_back(GenVertexData((float) j / (float) M, (float) i / (float) N));
+            vtxData.push_back(GenVertexData((float) j / (float) M, (float) (i + 1) / (float) N));
+        }
+    }
+    glBufferData(GL_ARRAY_BUFFER, nVtxPerStrip * nStrips * sizeof(VertexData), &vtxData[0], GL_STATIC_DRAW);
+    // Enable the vertex attribute arrays
+    glEnableVertexAttribArray(0);  // attribute array 0 = POSITION
+    glEnableVertexAttribArray(1);  // attribute array 1 = NORMAL
+    glEnableVertexAttribArray(2);  // attribute array 2 = TEXCOORD0
+    // attribute array, components/attribute, component type, normalize?, stride, offset
+    glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, sizeof(VertexData), (void *) offsetof(VertexData, position));
+    glVertexAttribPointer(1, 3, GL_FLOAT, GL_FALSE, sizeof(VertexData), (void *) offsetof(VertexData, normal));
+    glVertexAttribPointer(2, 2, GL_FLOAT, GL_FALSE, sizeof(VertexData), (void *) offsetof(VertexData, texcoord));
+}
+
+void ParamSurface::Draw() {
+    glBindVertexArray(vao);
+    for (unsigned int i = 0; i < nStrips; i++) glDrawArrays(GL_TRIANGLE_STRIP, i * nVtxPerStrip, nVtxPerStrip);
+}
+
+VertexData ParamSurface::GetVertexDataByUV(float u, float v) {
+    return GenVertexData(u, v);
+}
+
+Sphere::Sphere() {
+    create(tessellationLevel, tessellationLevel);
+}
+
+void Sphere::eval(Dnum2 &U, Dnum2 &V, Dnum2 &X, Dnum2 &Y, Dnum2 &Z) {
+    U = U * 2.0f * (float) M_PI;
+    V = V * (float) M_PI;
+
+    X = Cos(U) * Sin(V);
+    Y = Sin(U) * Sin(V);
+    Z = Cos(V);
+}
\ No newline at end of file
diff --git a/src/geometries/ParamSurface.h b/src/geometries/ParamSurface.h
new file mode 100644
index 0000000..3132668
--- /dev/null
+++ b/src/geometries/ParamSurface.h
@@ -0,0 +1,42 @@
+#ifndef BRAVE2_PARAMSURFACE_H
+#define BRAVE2_PARAMSURFACE_H
+
+#include "Geometry.h"
+
+//TODO put this in a better place
+const int tessellationLevel = 10;
+
+/**
+ * ParamSurface
+ */
+
+class ParamSurface : public Geometry {
+    unsigned int nVtxPerStrip, nStrips;
+public:
+
+    ParamSurface();
+
+    virtual void eval(Dnum2 &U, Dnum2 &V, Dnum2 &X, Dnum2 &Y, Dnum2 &Z) = 0;
+
+    virtual VertexData GenVertexData(float u, float v);
+
+    void create(int N = tessellationLevel, int M = tessellationLevel);
+
+    void Draw();
+
+    VertexData GetVertexDataByUV(float u, float v) override;
+};
+
+/**
+ * Sphere
+ */
+
+class Sphere : public ParamSurface {
+public:
+    Sphere();
+
+    void eval(Dnum2 &U, Dnum2 &V, Dnum2 &X, Dnum2 &Y, Dnum2 &Z);
+};
+
+
+#endif //BRAVE2_PARAMSURFACE_H
diff --git a/src/simulation/sph/sph_simulation.cpp b/src/geometries/SPHSimulation.cpp
similarity index 54%
rename from src/simulation/sph/sph_simulation.cpp
rename to src/geometries/SPHSimulation.cpp
index 0d6a8e1..86e3f0d 100644
--- a/src/simulation/sph/sph_simulation.cpp
+++ b/src/geometries/SPHSimulation.cpp
@@ -1,6 +1,8 @@
-#include "sph_simulation.h"
+#include "SPHSimulation.h"
 
-SPHSimulation::SPHSimulation() {
+#include <cmath>
+
+SPHSimulation::SPHSimulation() : externalForces(.0f, .0f, .0f) {
     std::cout << "Initializing sph particle system with " << sph::N <<
               " particles." << std::endl;
 
@@ -10,9 +12,11 @@ SPHSimulation::SPHSimulation() {
                 if (particles.size() < sph::N) {
                     particles.emplace_back(
                             new Particle(
-                                    glm::linearRand(1.f, 1.1f),
-                                    glm::linearRand(1.f, 1.1f),
-                                    glm::linearRand(1.f, 1.1f)
+                                    vec3(util::randomBetween(-0.15f, 0.15f),
+                                     util::randomBetween(0.3f, 0.35),
+                                     util::randomBetween(-0.15f, 0.15f)),
+                                    1.0,
+                                    vec3(0.5f, 0.5f, 0.8f)
                             )
                     );
                 }
@@ -22,57 +26,45 @@ SPHSimulation::SPHSimulation() {
 }
 
 SPHSimulation::~SPHSimulation() {
-    for (auto &particle : particles) { delete particle; }
+    for (auto &particle : particles)
+        delete particle;
 }
 
-void SPHSimulation::compute_density_and_pressure() {
-    for (auto &p_i : particles) {
-        // compute density
-        p_i->rho = 0.f;
-        for (auto &p_j : particles) {
-            glm::vec3 r_ij = p_j->pos - p_i->pos;
-            // euclidean distance between p_i and p_j
-            float r = glm::l2Norm(p_j->pos, p_i->pos);
-            float r2 = glm::pow(r, 2.f);
-
-            // 0 <= r <= d
-            if (r2 <= sph::H2)
-                p_i->rho += sph::M * sph::POLY6 * glm::pow(sph::H2 - r2, 3.f);
-        }
-
-        // compute pressure
-        p_i->p = sph::GAS_CONST * (p_i->rho - sph::REST_DENS);
-    }
+void SPHSimulation::update(float t) {
+    compute_density_and_pressure();
+    compute_forces();
 
+    integrate();
 }
 
-void SPHSimulation::compute_forces() {
-    for (auto &p_i : particles) {
-        glm::vec3 f_pressure(0.f, 0.f, 0.f);
-        glm::vec3 f_viscosity(0.f, 0.f, 0.f);
-
-        for (auto &p_j : particles) {
-            if (&p_i == &p_j) continue;
-
-            glm::vec3 r_ijv = p_i->pos - p_j->pos;
-            float r_ij = glm::length(r_ijv);
-
-            if (r_ij <= sph::H) {
-                // compute pressure force contribution
+void SPHSimulation::Draw() {
+    std::vector<float> particlePosAndColor;
+    for (auto &p : particles) {
+        particlePosAndColor.push_back(p->pos.x);
+        particlePosAndColor.push_back(p->pos.y);
+        particlePosAndColor.push_back(p->pos.z);
+        particlePosAndColor.push_back(p->color.x);
+        particlePosAndColor.push_back(p->color.y);
+        particlePosAndColor.push_back(p->color.z);
+    }
 
-                f_pressure += -sph::M * (p_i->p + p_j->p) / (2 * p_j->rho) *
-                              sph::POLY6_GRAD_PRESS * r_ijv / r_ij * glm::pow(sph::H - r_ij, 2.0f);
+    glBindVertexArray(vao);
+    glBindBuffer(GL_ARRAY_BUFFER, vbo);
+    // _STREAM_ indicates that the data will be modified, and then used only a few (one) time.
+    // The buffer usage token might be changed to GL_DYNAMIC_DRAW when moving on to simulating positions on the GPU
+    glBufferData(GL_ARRAY_BUFFER, particlePosAndColor.size() * sizeof(float), particlePosAndColor.data(),
+                 GL_STREAM_DRAW);
 
-                f_viscosity += sph::MU * sph::M * (p_j->v - p_i->v) / p_j->rho *
-                               sph::POLY6_GRAD_VISC * (sph::H - r_ij);
-            }
+    // position attribute
+    glEnableVertexAttribArray(0);
+    glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, sizeof(float) * 6, (void *) 0);
 
-            glm::vec3 f_gravity = sph::G * p_i->rho;
+    // color attribute
+    glEnableVertexAttribArray(1);
+    glVertexAttribPointer(1, 3, GL_FLOAT, GL_FALSE, sizeof(float) * 6, (void *) (sizeof(float) * 3));
 
-            p_i->f = f_pressure + f_viscosity + f_gravity;
-            p_i->f = f_gravity + f_pressure;// * 0.001f;
-        }
-    }
+    glPointSize(3.f);
+    glDrawArrays(GL_POINTS, 0, particlePosAndColor.size() / 6);
 }
 
 void SPHSimulation::integrate() {
@@ -81,8 +73,8 @@ void SPHSimulation::integrate() {
         p->v += sph::DT * p->f / p->rho;
         p->pos += sph::DT * p->v;
 
-        const float max_b = 1.2f;
-        const float min_b = 0.9f;
+        const float max_b = 0.15f;
+        const float min_b = -0.15f;
 
         //boundary conditions
         if (p->pos.x - sph::EPS <= min_b) {
@@ -119,29 +111,67 @@ void SPHSimulation::integrate() {
 
         p->v *= 0.98;
     }
+
 }
 
-void SPHSimulation::update(float t) {
-    compute_density_and_pressure();
-    compute_forces();
+void SPHSimulation::compute_forces() {
+    for (auto &p_i : particles) {
+        vec3 f_pressure(0.f, 0.f, 0.f);
+        vec3 f_viscosity(0.f, 0.f, 0.f);
 
-//        if(t > sph::DT)
-//            for(float dt = 0.0f; dt < t; dt+=sph::DT)
-//                integrate();
-//        else
-    integrate();
+        for (auto &p_j : particles) {
+            if (&p_i == &p_j) continue;
+
+            vec3 r_ijv = p_i->pos - p_j->pos;
+            float r_ij = length(r_ijv);
+
+            if (r_ij <= sph::H) {
+                // compute pressure force contribution
+
+                f_pressure += -sph::M * (p_i->p + p_j->p) / (2 * p_j->rho) *
+                              sph::POLY6_GRAD_PRESS * r_ijv / r_ij * std::pow(sph::H - r_ij, 2.0f);
+
+                f_viscosity += sph::MU * sph::M * (p_j->v - p_i->v) / p_j->rho *
+                               sph::POLY6_GRAD_VISC * (sph::H - r_ij);
+            }
+
+            vec3 f_gravity = sph::G * p_i->rho;
+
+            p_i->f = f_pressure + f_viscosity + f_gravity;
+            p_i->f = f_gravity + f_pressure;// * 0.001f;
+        }
+    }
 }
 
-void SPHSimulation::render() {
-    // TODO MODERN OPENGL!!!!
-    glEnable(GL_POINT_SMOOTH);
-    glPointSize(1.f);
+void SPHSimulation::compute_density_and_pressure() {
+    for (auto &p_i : particles) {
+        // compute density
+        p_i->rho = 0.f;
+        for (auto &p_j : particles) {
+            vec3 r_ij = p_j->pos - p_i->pos;
+            // euclidean distance between p_i and p_j
+            float r = length(p_j->pos - p_i->pos);
+            float r2 = std::pow(r, 2.f);
+
+            // 0 <= r <= d
+            if (r2 <= sph::H2)
+                p_i->rho += sph::M * sph::POLY6 * std::pow(sph::H2 - r2, 3.f);
+        }
+
+        // compute pressure
+        p_i->p = sph::GAS_CONST * (p_i->rho - sph::REST_DENS);
+    }
+
+}
 
-    glColor4f(0.5f, 0.5f, 0.8f, 1.f);
+vec3 SPHSimulation::getExternalForces() const {
+    return vec3();
+}
 
-    glBegin(GL_POINTS);
-    for (auto &p : particles)
-        glVertex3f(p->pos.x, p->pos.y, p->pos.z);
-    glEnd();
+void SPHSimulation::resetExternalForces() {
+    externalForces = vec3(0, 0, 0);
+}
 
+void SPHSimulation::addForce(vec3 force) {
+    externalForces += force;
 }
diff --git a/src/simulation/sph/sph_simulation.h b/src/geometries/SPHSimulation.h
similarity index 55%
rename from src/simulation/sph/sph_simulation.h
rename to src/geometries/SPHSimulation.h
index 4b3c30c..b23e61e 100644
--- a/src/simulation/sph/sph_simulation.h
+++ b/src/geometries/SPHSimulation.h
@@ -1,21 +1,14 @@
-#ifndef MADID_SPH_SYSTEM
-#define MADID_SPH_SYSTEM
+#ifndef BRAVE2_SPHSIMULATION_H
+#define BRAVE2_SPHSIMULATION_H
 
-#include <glm/glm.hpp>
-#include <glm/gtc/random.hpp>
-#include <glm/gtx/norm.hpp>
-#include <iostream>
-#include <vector>
-
-#include "../../utility/gl.h"
-
-#include "../particle.h"
-#include "../../constants.h"
-#include "../simulation.h"
+#include "../utils/math.h"
+#include "../utils/util.h"
+#include "../objects/Particle.h"
+#include "Geometry.h"
 
 namespace sph {
     // external (gravitational force)
-    const static glm::vec3 G(0.f, -9.8f, 0.f);
+    const static vec3 G(0.f, -9.8f, 0.f);
     // rest density
     const static float REST_DENS = 1000.0f;
     // const for equation of state
@@ -37,13 +30,11 @@ namespace sph {
     // damping when colliding with the boundary
     const static float BOUND_DAMPING = -0.5f;
 
-    [[maybe_unused]] const static float PI = glm::pi<float>();
+    [[maybe_unused]] const static float PI = 3.1415926535f;
     // smoothing kernel as described in the Müller paper
-    const static float POLY6 = 315.f / (64.f * glm::pi<float>() * glm::pow(H, 9));
-    const static float POLY6_GRAD_PRESS = -45.f / (glm::pi<float>()
-                                                   * glm::pow(H, 6));
-    const static float POLY6_GRAD_VISC = 45.f / (glm::pi<float>() * glm::pow(H,
-                                                                             6));
+    const static float POLY6 = 315.f / (64.f * PI * pow(H, 9));
+    const static float POLY6_GRAD_PRESS = -45.f / (PI * pow(H, 6));
+    const static float POLY6_GRAD_VISC = 45.f / (PI * pow(H,6));
 
     // PARTICLES
     // number of particles
@@ -51,9 +42,9 @@ namespace sph {
 }
 
 
-class SPHSimulation : Simulation {
-
+class SPHSimulation : public Geometry {
     std::vector<Particle *> particles;
+    vec3 externalForces;
 
     void compute_density_and_pressure();
 
@@ -67,9 +58,16 @@ public:
     ~SPHSimulation();
 
     // step the system by t time
-    void update(float t) override;
+    void update(float t);
+
+    void Draw() override;
 
-    void render() override;
+    vec3 getExternalForces() const;
+
+    void resetExternalForces();
+
+    void addForce(vec3 force);
 };
 
-#endif
+
+#endif //BRAVE2_SPHSIMULATION_H
diff --git a/src/geometries/VertexData.cpp b/src/geometries/VertexData.cpp
new file mode 100644
index 0000000..9f8afec
--- /dev/null
+++ b/src/geometries/VertexData.cpp
@@ -0,0 +1,13 @@
+#include "VertexData.h"
+
+std::ostream &operator<<(std::ostream &out, const VertexData &vD) {
+    return out << vD.texcoord << " -> P: " << vD.position << " N: " << vD.normal;
+}
+
+VertexData::VertexData(vec3 p, vec3 n, vec2 uv) :
+        position(p), normal(n), texcoord(uv) {
+
+}
+
+VertexData::VertexData() {
+}
diff --git a/src/geometries/VertexData.h b/src/geometries/VertexData.h
new file mode 100644
index 0000000..5566618
--- /dev/null
+++ b/src/geometries/VertexData.h
@@ -0,0 +1,18 @@
+#ifndef BRAVE2_VERTEXDATA_H
+#define BRAVE2_VERTEXDATA_H
+
+#include "../utils/math.h"
+
+struct VertexData {
+    // needed for GetVertexByUV function
+    bool valid = true;
+    vec3 position, normal;
+    vec2 texcoord;
+    VertexData();
+    VertexData(vec3 p, vec3 n, vec2 uv);
+    explicit VertexData(bool valid, vec3 p = {0,0,0}, vec3 n = {0,0,0}, vec2 uv = {0,0}){};
+};
+
+std::ostream &operator<<(std::ostream &out, const VertexData &vD);
+
+#endif //BRAVE2_VERTEXDATA_H
diff --git a/src/main.cpp b/src/main.cpp
index dc3a544..cdd3eab 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,11 +1,134 @@
-#define GLFW_INCLUDE_NONE
 #include <GL/glew.h>
 #include <GLFW/glfw3.h>
-
+#include <cstdlib>
 #include <iostream>
-#include <memory>
 
-#include "application.h"
+#include "utils/math.h"
+
+#include "utils/other/stb_image_write.h"
+#include "utils/save_image.h"
+#include "rendering/Scene.h"
+#include "Constants.h"
+#include "utils/InputHandler.h"
+
+GLFWwindow *window;
+
+bool dragging = false;
+int keyArr[350];
+
+InputHandler *InputHandler = InputHandler::GetInstance();
+
+vec3 forceGenerated(0.0f, 0.0f, 0.0f);
+bool resetExternalForces = false;
+bool capturing = false;
+
+
+Scene Scene(WIDTH, HEIGHT);
+
+
+static void Initialize() {
+    glViewport(0, 0, WIDTH, HEIGHT);
+    glMatrixMode(GL_MODELVIEW);
+    glEnable(GL_DEPTH_TEST);
+    glDepthMask(GL_TRUE);
+    glLoadIdentity();
+//    glClearColor(196.0f / 255.0f, 233.0f / 255.0f, 241.0f / 255.0f, 1.0f);
+    glClearColor(1,1,1, 1.0f);
+}
+
+static void Update(GLFWwindow *window, float delta) {
+    if (keyArr[GLFW_KEY_ESCAPE])
+        glfwSetWindowShouldClose(window, 1);
+}
+
+static void RenderScene(GLFWwindow *window, float delta) {
+
+}
+
+static void Resize(GLFWwindow *window, int w, int h) {
+    if (h < 1)
+        h = 1;
+    glViewport(0, 0, w, h);
+    glMatrixMode(GL_PROJECTION);
+    glLoadIdentity();
+    gluPerspective(45.0f, (float) w / (float) h, 0.1f, 1000.0f);
+    gluLookAt(0.0f, 0.0f, 30, 0.0f, 0.0f, 0.0f, 0.0f, 1.0f, 0.0f);
+    glMatrixMode(GL_MODELVIEW);
+}
+
+static void getForce(char dir, float power = 0.05f) {
+    vec3 force(0.0f, 0.0f, 0.0f);
+    switch (dir) {
+        case 'x':
+            force = power * vec3(-1.f, 0.f, .0f);
+            break;
+        case 'X':
+            force = power * vec3(1.f, 0.f, .0f);
+            break;
+        case 'y':
+            force = power * vec3(.0f, -1.f, .0f);
+            break;
+        case 'Y':
+            force = power * vec3(.0f, 1.f, .0f);
+            break;
+        case 'z':
+            force = power * vec3(.0f, .0f, -1.f);
+            break;
+        case 'Z':
+            force = power * vec3(.0f, .0f, 1.f);
+            break;
+        default:
+            break;
+    }
+    forceGenerated += force;
+}
+
+static void KeyCallback(GLFWwindow *window, int key, int scancode, int action, int mods) {
+    if (action == GLFW_RELEASE)
+        InputHandler->KeyRelease(key);
+    else
+        InputHandler->KeyPress(key);
+
+    InputHandler->SetModifiers(mods);
+
+    if (action == GLFW_RELEASE && key == GLFW_KEY_C) {
+        if (capturing) capturing = false;
+        else capturing = true;
+    }
+
+}
+
+static void MouseClickCallback(GLFWwindow *window, int button, int action, int mods) {
+    switch (button) {
+        case GLFW_MOUSE_BUTTON_1:
+            dragging = (action == GLFW_PRESS);
+            break;
+    }
+}
+
+struct Mouse {
+    float x;
+    float y;
+
+    Mouse(float x, float y) {
+        this->x = x / WIDTH;
+        this->y = 1 - (y / HEIGHT);
+    }
+};
+
+Mouse prevMousePos(0, 0);
+
+static void MouseMotionCallback(GLFWwindow *window, double x, double y) {
+    if (!dragging) {
+        prevMousePos = Mouse(x, y);
+        return;
+    }
+    Mouse currMouse((float) x, (float) y);
+    vec2 deltaMouse(currMouse.x - prevMousePos.x, currMouse.y - prevMousePos.y);
+    prevMousePos = currMouse;
+    vec3 force((float) deltaMouse.x / 1.0f, (float) deltaMouse.y / 1.0f, 0.0f);
+    forceGenerated += force;
+}
 
 
 int main(int argc, char **argv) {
@@ -15,32 +138,101 @@ int main(int argc, char **argv) {
         return 1;
     }
 
-    /* Creating a GLFW application */
-    /* Application creates its GLFW window too */
-    std::shared_ptr<Application> application = std::make_shared<Application>();
-
-    if (!application->window) {
+    glfwInit();
+    window = glfwCreateWindow(WIDTH, HEIGHT, "Brave-2", NULL, NULL);
+    if (!window) {
         std::cerr << "ERROR: could not open window with GLFW3" << std::endl;
         glfwTerminate();
         return 1;
     }
 
-    glfwMakeContextCurrent(application->get_window());
+    glfwMakeContextCurrent(window);
 
     // start GLEW extension handler
-    if(glewInit() != GLEW_NO_ERROR) {
-        std::cerr << "Failed to initialize GLEW" << std::endl;
-        return -1;
-    }
     glewExperimental = GL_TRUE;
+    glewInit();
 
     // get version info
     const GLubyte *renderer = glGetString(GL_RENDERER); // get renderer string
-    const GLubyte *version = glGetString(GL_VERSION); // version as a string/home/bobarna/codes/madid
+    const GLubyte *version = glGetString(GL_VERSION); // version as a string
     std::cerr << "Renderer: " << renderer << std::endl;
     std::cerr << "OpenGL version supported: " << version << std::endl;
 
-    application->init();
-    application->start();
+    glfwSetWindowSizeCallback(window, Resize);
+    glfwSetKeyCallback(window, KeyCallback);
+    glfwSetMouseButtonCallback(window, MouseClickCallback);
+    glfwSetCursorPosCallback(window, MouseMotionCallback);
+
+    srand(time(NULL));
+
+    // Initialize
+    Initialize();
+    Scene.Build();
+
+    int imageNr = 0;
+    double timeElapsed;
+
+    // MAIN LOOP
+    while (!glfwWindowShouldClose(window)) {
+        // ticking every 24 FPS
+        bool tick = false;
+
+        double deltaTime = glfwGetTime() - timeElapsed;
+        timeElapsed = glfwGetTime();
+
+        if (timeElapsed >= 1.0f / 24.0f) {
+            glfwSetTime(0.0f);
+            timeElapsed = 0.0f;
+            tick = true;
+        }
+
+        //clear color and depth buffer
+        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+        glLoadIdentity(); //load identity matrix
+
+        // reset all external forces if gravity was toggled
+        if (resetExternalForces || !dragging) {
+            Scene.ResetExternalForces();
+            resetExternalForces = false;
+        }
+
+        //if there was force generated...
+        // ... add it to the simulation
+        if (forceGenerated.x != 0.0f ||
+            forceGenerated.y != 0.0f ||
+            forceGenerated.z != 0.0f) {
+            Scene.addForce(forceGenerated);
+            // ... and reset the force generated
+            forceGenerated = vec3(0, 0, 0);
+        }
+
+
+        //TODO guarantee dt to be infinitesimal
+        Scene.Update(deltaTime);
+        Scene.Render();
+
+        // big red dot in lower left corner
+//        if (capturing) {
+//            glPointSize(20.0f);
+//            glBegin(GL_POINTS);
+//            glColor3f(200.0f, 0.0f, 0.0f);
+//            glVertex3f(-0.8f, -0.8f, 1.0f);
+//            glEnd();
+//        }
+
+        if (tick && capturing) {
+            char path[100];
+            sprintf(path, "../renders/render%04d.bmp", imageNr++);
+            saveImage(path, window);
+            std::cout << path << ".bmp printed" << std::endl;
+        }
+
+
+        glfwSwapBuffers(window);
+        glfwPollEvents();
+    }
 
+    // close GL context and any other GLFW resources
+    glfwTerminate();
+    return 0;
 }
diff --git a/src/objects/HeadObject.cpp b/src/objects/HeadObject.cpp
new file mode 100644
index 0000000..50f1798
--- /dev/null
+++ b/src/objects/HeadObject.cpp
@@ -0,0 +1,23 @@
+#include "HeadObject.h"
+
+
+HeadObject::HeadObject(Shader *_shader, ObjGeometry *_geometry, Material *_material, Texture *_texture) :
+        Object(_shader, _geometry, _material, _texture) {}
+
+VertexData HeadObject::GetVertexDataByUV(float u, float v) {
+    mat4 M, Minv;
+    SetModelingTransform(M, Minv);
+    VertexData vD{reinterpret_cast<ObjGeometry*>(geometry)->GetVertexDataByUV(u, v)};
+    vec4 WP = vec4(vD.position.x, vD.position.y, vD.position.z, 1) * M;
+    vec4 WN = vec4(vD.normal.x, vD.normal.y, vD.normal.z, 0) * M;
+    vD.position = vec3(WP.x, WP.y, WP.z);
+    vD.normal = vec3(WN.x, WN.y, WN.z);
+
+    return vD;
+}
+
+ObjGeometry *HeadObject::getGeometry() {
+    return reinterpret_cast<ObjGeometry*>(geometry);
+}
+
+
diff --git a/src/objects/HeadObject.h b/src/objects/HeadObject.h
new file mode 100644
index 0000000..5c2bb54
--- /dev/null
+++ b/src/objects/HeadObject.h
@@ -0,0 +1,17 @@
+#ifndef BRAVE2_HEADOBJECT_H
+#define BRAVE2_HEADOBJECT_H
+
+#include "Object.h"
+#include "../geometries/ObjGeometry.h"
+
+class HeadObject : public Object {
+public:
+    HeadObject(Shader *_shader, ObjGeometry *_geometry, Material *_material, Texture *_texture);
+
+    VertexData GetVertexDataByUV(float u, float v);
+
+    ObjGeometry* getGeometry();
+};
+
+
+#endif //BRAVE2_HEADOBJECT_H
diff --git a/src/rendering/objects/object.cpp b/src/objects/Object.cpp
similarity index 89%
rename from src/rendering/objects/object.cpp
rename to src/objects/Object.cpp
index 8fe05cd..bfbd120 100644
--- a/src/rendering/objects/object.cpp
+++ b/src/objects/Object.cpp
@@ -1,4 +1,4 @@
-#include "object.h"
+#include "Object.h"
 
 Object::Object(Shader *_shader, Geometry *_geometry, Material *_material, Texture *_texture) :
         scale(vec3(1, 1, 1)), translation(vec3(0, 0, 0)), rotationAxis(0, 0, 1), rotationAngle(0) {
@@ -14,7 +14,7 @@ void Object::SetModelingTransform(mat4 &M, mat4 &Minv) {
            scale_matrix(vec3(1 / scale.x, 1 / scale.y, 1 / scale.z));
 }
 
-void Object::render(RenderState state) {
+void Object::Draw(RenderState state) {
     mat4 M, Minv;
     SetModelingTransform(M, Minv);
     state.M = M;
@@ -23,10 +23,10 @@ void Object::render(RenderState state) {
     state.material = material;
     state.texture = texture;
     shader->Bind(state);
-    geometry->render();
+    geometry->Draw();
 }
 
-void Object::update(float delta_t) {
+void Object::Animate(float delta_t) {
 }
 
 
@@ -41,4 +41,5 @@ void Object::Translate(vec3 t) {
 void Object::Rotate(vec3 axis, float angle) {
     rotationAxis = axis;
     rotationAngle = angle;
-}
\ No newline at end of file
+}
+
diff --git a/src/rendering/objects/object.h b/src/objects/Object.h
similarity index 59%
rename from src/rendering/objects/object.h
rename to src/objects/Object.h
index 1eef817..ec799b5 100644
--- a/src/rendering/objects/object.h
+++ b/src/objects/Object.h
@@ -1,13 +1,9 @@
-#ifndef MADID_OBJECT_H
-#define MADID_OBJECT_H
+#ifndef BRAVE2_OBJECT_H
+#define BRAVE2_OBJECT_H
 
-#include <glm/glm.hpp>
-#include "../../utility/custom_math.h"
-#include "../geometry.h"
-#include "../shaders/shader.h"
-#include "../materials/material.h"
-
-using namespace glm;
+#include "../utils/math.h"
+#include "../geometries/Geometry.h"
+#include "../rendering/shaders/Shader.h"
 
 class Object {
 protected:
@@ -23,9 +19,9 @@ public:
 
     virtual void SetModelingTransform(mat4 &M, mat4 &Minv);
 
-    void render(RenderState state);
+    void Draw(RenderState state);
 
-    virtual void update(float delta_t);
+    virtual void Animate(float delta_t);
 
     void Scale(vec3 s);
 
@@ -35,4 +31,4 @@ public:
 };
 
 
-#endif //MADID_OBJECT_H
+#endif //BRAVE2_OBJECT_H
diff --git a/src/objects/PBDSimulationObject.cpp b/src/objects/PBDSimulationObject.cpp
new file mode 100644
index 0000000..1e62327
--- /dev/null
+++ b/src/objects/PBDSimulationObject.cpp
@@ -0,0 +1,22 @@
+#include "PBDSimulationObject.h"
+
+PBDSimulationObject::PBDSimulationObject(Shader *_shader, PBDSimulation *_sim, Material* _material, Texture* _texture) :
+        Object(_shader, _sim, _material, _texture) {
+}
+
+void PBDSimulationObject::Animate(float delta_t) {
+    reinterpret_cast<PBDSimulation *>(geometry)->update(delta_t);
+}
+
+void PBDSimulationObject::ResetExternalForces() {
+    reinterpret_cast<PBDSimulation *>(geometry)->resetExternalForces();
+}
+
+void PBDSimulationObject::AddForce(vec3 f) {
+    reinterpret_cast<PBDSimulation *>(geometry)->addForce(f);
+}
+
+PBDSimulationObject::PBDSimulationObject(Shader *_shader, PBDSimulation *_sim) :
+        Object(_shader, _sim){
+
+}
diff --git a/src/objects/PBDSimulationObject.h b/src/objects/PBDSimulationObject.h
new file mode 100644
index 0000000..2fba8b6
--- /dev/null
+++ b/src/objects/PBDSimulationObject.h
@@ -0,0 +1,24 @@
+#ifndef BRAVE2_PBDSIMULATIONOBJECT_H
+#define BRAVE2_PBDSIMULATIONOBJECT_H
+
+#include "../geometries/PBDSimulation.h"
+#include "../rendering/shaders/Shader.h"
+#include "Object.h"
+#include "HeadObject.h"
+
+class PBDSimulationObject : public Object {
+
+public:
+    PBDSimulationObject(Shader *_shader, PBDSimulation *_sim, Material* _material, Texture* _texture);
+
+    PBDSimulationObject(Shader *_shader, PBDSimulation *_sim);
+
+    void ResetExternalForces();
+
+    void AddForce(vec3 f);
+
+    void Animate(float delta_t) override;
+};
+
+
+#endif //BRAVE2_PBDSIMULATIONOBJECT_H
diff --git a/src/objects/Particle.h b/src/objects/Particle.h
new file mode 100644
index 0000000..de2b3a6
--- /dev/null
+++ b/src/objects/Particle.h
@@ -0,0 +1,42 @@
+#ifndef BRAVE2_PARTICLE_H
+#define BRAVE2_PARTICLE_H
+
+#include "../utils/math.h"
+
+class Particle {
+public:
+    Particle(vec3 pos, float w, vec3 c) : pos(pos), tmp_pos(pos), w(w), v(0.0f), color(c) {
+
+    };
+
+    Particle(vec3 pos): pos(pos) {
+
+    };
+
+    vec3 pos;
+    vec3 tmp_pos;
+
+    float w; // inverse mass of the particle
+
+    vec3 v; //velocity
+
+    vec3 color;
+
+    // For SPH simulation
+    float p;            // pressure (varies)
+    float rho;          // density (varies)
+    vec3 f;           // current sum of external forces acting on the particle
+
+
+    /** gets the particle's mass
+     *
+     * @return the particle's mass. If infinity, -1 is returned.
+     */
+    float get_mass() {
+        if (w == 0.0f) return -1;
+        return 1 / w;
+    }
+};
+
+
+#endif //BRAVE2_PARTICLE_H
diff --git a/src/objects/SPHSimulationObject.cpp b/src/objects/SPHSimulationObject.cpp
new file mode 100644
index 0000000..9ceb8d9
--- /dev/null
+++ b/src/objects/SPHSimulationObject.cpp
@@ -0,0 +1,17 @@
+#include "SPHSimulationObject.h"
+
+SPHSimulationObject::SPHSimulationObject(Shader *_shader, SPHSimulation *_sim) :
+        Object(_shader, _sim) {
+}
+
+void SPHSimulationObject::Animate(float delta_t) {
+    reinterpret_cast<SPHSimulation *>(geometry)->update(delta_t);
+}
+
+void SPHSimulationObject::ResetExternalForces() {
+    reinterpret_cast<SPHSimulation *>(geometry)->resetExternalForces();
+}
+
+void SPHSimulationObject::AddForce(vec3 f) {
+    reinterpret_cast<SPHSimulation *>(geometry)->addForce(f);
+}
\ No newline at end of file
diff --git a/src/objects/SPHSimulationObject.h b/src/objects/SPHSimulationObject.h
new file mode 100644
index 0000000..00ca9c7
--- /dev/null
+++ b/src/objects/SPHSimulationObject.h
@@ -0,0 +1,22 @@
+#ifndef BRAVE2_SPHSIMULATIONOBJECT_H
+#define BRAVE2_SPHSIMULATIONOBJECT_H
+
+#include "../rendering/shaders/Shader.h"
+#include "Object.h"
+#include "HeadObject.h"
+#include "../geometries/SPHSimulation.h"
+
+class SPHSimulationObject : public Object {
+
+public:
+    SPHSimulationObject(Shader *_shader, SPHSimulation *_sim);
+
+    void ResetExternalForces();
+
+    void AddForce(vec3 f);
+
+    void Animate(float delta_t) override;
+};
+
+
+#endif //BRAVE2_SPHSIMULATIONOBJECT_H
diff --git a/src/rendering/Camera.cpp b/src/rendering/Camera.cpp
new file mode 100644
index 0000000..f210513
--- /dev/null
+++ b/src/rendering/Camera.cpp
@@ -0,0 +1,46 @@
+#include "Camera.h"
+
+Camera::Camera(vec3 _wEye, vec3 _wLookAt, vec3 _wVup, int windowWidth, int windowHeight) {
+    asp = (float) windowWidth / (float) windowHeight;
+    fov = 75.0f * (float) M_PI / 180.0f;
+    fp = 0.01f; //TODO is this right??
+    bp = 10;
+
+    wEye = _wEye;
+    wLookAt = _wLookAt;
+    wVup = _wVup;
+}
+
+mat4 Camera::V() const {
+    vec3 w = normalize(wEye - wLookAt);
+    vec3 u = normalize(cross(wVup, w));
+    vec3 v = cross(w, u);
+    return translate_matrix(-wEye) * mat4(
+            u.x, v.x, w.x, 0,
+            u.y, v.y, w.y, 0,
+            u.z, v.z, w.z, 0,
+            0, 0, 0, 1
+    );
+}
+
+mat4 Camera::P() const {
+    return mat4(1.0f / (tanf(fov / 2) * asp), 0, 0, 0,
+                0, 1.0f / tanf(fov / 2), 0, 0,
+                0, 0, -(fp + bp) / (bp - fp), -1,
+                0, 0, -2 * fp * bp / (bp - fp), 0);
+}
+
+void Camera::Translate(vec3 dir) {
+    vec3 forward = dir * vec3(0, 0, 1.f);
+
+    wEye = normalize(wEye + dir) * length(wEye) + forward;
+}
+
+RenderState Camera::getState() {
+    RenderState state;
+    state.wEye = wEye;
+    state.V = V();
+    state.P = P();
+
+    return state;
+}
diff --git a/src/rendering/Camera.h b/src/rendering/Camera.h
new file mode 100644
index 0000000..9b24161
--- /dev/null
+++ b/src/rendering/Camera.h
@@ -0,0 +1,52 @@
+#ifndef BRAVE2_CAMERA_H
+#define BRAVE2_CAMERA_H
+
+
+#include "../utils/math.h"
+#include "RenderState.h"
+
+class Camera {
+public:
+    /// Position of the camera in world coordinates.
+    vec3 wEye;
+
+    /// Location of the center of the window.
+    vec3 wLookAt;
+
+    /** @brief
+     * To find the vertical position of the window.
+     * If wVup is not perpendicular to (wLookAt-wEye),
+     * then only it's perpendicular position is used.
+     */
+    vec3 wVup;
+
+    // Field of View - vertical size of the window
+    float fov;
+
+    // Aspect Ratio - vertical/horizontal
+    float asp;
+
+    // Front Clipping Plane
+    float fp;
+
+    // Back Clipping Plane
+    float bp;
+
+public:
+    Camera(vec3 _wEye, vec3 _wLookAt, vec3 _wVup, int windowWidth, int windowHeight);
+
+    // View Matrix: translates the center to the origin
+    mat4 V() const;
+
+    // Projection Matrix
+    mat4 P() const;
+
+    RenderState getState();
+
+    // Translates the camera by dir
+    void Translate(vec3 dir);
+
+};
+
+
+#endif //BRAVE2_CAMERA_H
diff --git a/src/rendering/RenderState.h b/src/rendering/RenderState.h
new file mode 100644
index 0000000..057af06
--- /dev/null
+++ b/src/rendering/RenderState.h
@@ -0,0 +1,18 @@
+#ifndef BRAVE2_RENDERSTATE_H
+#define BRAVE2_RENDERSTATE_H
+
+#include "../utils/math.h"
+#include "../utils/GPUProgram.h"
+#include "materials/Material.h"
+#include "lights/Light.h"
+
+struct RenderState {
+    mat4 MVP, M, Minv, V, P;
+    vec3 wEye;
+    Texture *texture;
+    Material *material;
+
+    std::vector<Light> lights;
+};
+
+#endif //BRAVE2_RENDERSTATE_H
diff --git a/src/rendering/Scene.cpp b/src/rendering/Scene.cpp
new file mode 100644
index 0000000..b9a90a2
--- /dev/null
+++ b/src/rendering/Scene.cpp
@@ -0,0 +1,113 @@
+#include "Scene.h"
+
+
+Scene::Scene(int w, int h) : camera(vec3(0, 0.f, .5), // Camera position (wEye)
+                                    vec3(0, 0, 0), // wLookat
+                                    vec3(0, 1, 0), // wVup
+                                    w, h)    {
+}
+
+void Scene::Build() {
+    size_t nrSims = 100;
+    size_t nrSegments = 50;
+    float lSeg = 0.006f;
+
+    // Lights
+    lights.resize(3);
+    lights[0].wLightPos = vec4(0.0f, 3.0f, 1.5f, 0);    // ideal point -> directional light source
+    lights[0].La = vec3(0.0f, 7.0f, 7.0f);
+    lights[0].Le = vec3(1, 1, 1);
+
+    lights[1].wLightPos = vec4(2.0f, 3.0f, 1.0f, 0);    // ideal point -> directional light source
+    lights[1].La = vec3(0.f, 3.5f, 3.5f);
+    lights[1].Le = vec3(1, 1, 1);
+
+    lights[2].wLightPos = vec4(0, 0, 1, 0);    // ideal point -> directional light source
+    lights[2].La = vec3(0.f, 3.5f, 3.5f);
+    lights[2].Le = vec3(1, 1, 1);
+
+    Shader *basicShader = new BasicShader();
+    basicShader->Bind(camera.getState());
+
+    Shader *phongShader = new PhongShader();
+
+    auto PBDSim = new PBDSimulation(nrSims, nrSegments, lSeg);
+    auto PBDsimulationObject = new PBDSimulationObject(basicShader, PBDSim);
+    PBDsims.push_back(PBDsimulationObject);
+
+    auto SPHSim = new SPHSimulation();
+    auto SPHsimulationObject = new SPHSimulationObject(basicShader, SPHSim);
+    SPHsims.push_back(SPHsimulationObject);
+
+}
+
+void Scene::Render() {
+    RenderState state = camera.getState();
+    state.lights = lights;
+
+    for (auto *so: PBDsims) so->Draw(state);
+    for (auto *so : SPHsims) so->Draw(state);
+    for (auto *o: objects) o->Draw(state);
+}
+
+//delta_t is infinitesimal
+void Scene::Update(float delta_t) {
+    // TODO if there was a keypress
+    // Handle keypresses
+    HandleKeyPress();
+
+    for (Object *o: objects) o->Animate(delta_t);
+    for (PBDSimulationObject *so : PBDsims) so->Animate(delta_t);
+    for (SPHSimulationObject *so : SPHsims) so->Animate(delta_t);
+}
+
+void Scene::TranslateCamera(vec3 t) {
+    camera.Translate(t);
+}
+
+void Scene::ResetExternalForces() {
+    for (auto sim : PBDsims)
+        sim->ResetExternalForces();
+    if (gravityOn)
+        for (auto sim : PBDsims)
+            sim->AddForce({0, -GRAVITY_ABS_VALUE, 0});
+}
+
+void Scene::addForce(vec3 f) {
+    for (auto sim : PBDsims)
+        sim->AddForce(f);
+}
+
+void Scene::HandleKeyPress() {
+    // G turns gravity on
+    if (inputHandler->IsPressed(GLFW_KEY_G) && !inputHandler->IsAltPressed())
+        gravityOn = true;
+    // Alt+G turns gravity off
+    if (inputHandler->IsPressed(GLFW_KEY_G) && inputHandler->IsAltPressed())
+        gravityOn = false;
+    if (inputHandler->IsPressed(GLFW_KEY_W) ||
+        inputHandler->IsPressed(GLFW_KEY_A) ||
+        inputHandler->IsPressed(GLFW_KEY_S) ||
+        inputHandler->IsPressed(GLFW_KEY_D) ||
+        inputHandler->IsPressed(GLFW_KEY_Q) ||
+        inputHandler->IsPressed(GLFW_KEY_E))
+        HandleCameraMove();
+}
+
+void Scene::HandleCameraMove() {
+    vec3 vec = vec3(0, 0, 0);
+    if (inputHandler->IsPressed(GLFW_KEY_W))
+        vec = vec3(0, 0, -1);
+    if (inputHandler->IsPressed(GLFW_KEY_A))
+        vec = vec3(-1, 0, 0);
+    if (inputHandler->IsPressed(GLFW_KEY_S))
+        vec = vec3(0, 0, 1);
+    if (inputHandler->IsPressed(GLFW_KEY_D))
+        vec = vec3(1, 0, 0);
+    if (inputHandler->IsPressed(GLFW_KEY_Q))
+        vec = vec3(0, 1, 0);
+    if (inputHandler->IsPressed(GLFW_KEY_E))
+        vec = vec3(0, -1, 0);
+
+    TranslateCamera(vec * 0.1);
+}
diff --git a/src/rendering/Scene.h b/src/rendering/Scene.h
new file mode 100644
index 0000000..db65def
--- /dev/null
+++ b/src/rendering/Scene.h
@@ -0,0 +1,53 @@
+#ifndef BRAVE2_SCENE_H
+#define BRAVE2_SCENE_H
+
+#include <vector>
+#include "../geometries/PBDSimulation.h"
+#include "../objects/PBDSimulationObject.h"
+#include "Camera.h"
+#include "../objects/Object.h"
+#include "../Constants.h"
+#include "../utils/InputHandler.h"
+#include "../geometries/SPHSimulation.h"
+#include "../objects/SPHSimulationObject.h"
+#include "shaders/BasicShader.h"
+#include "../geometries/ParamSurface.h"
+#include "shaders/PhongShader.h"
+#include "../geometries/ObjGeometry.h"
+
+
+#include <GLFW/glfw3.h>
+
+class Scene {
+    std::vector<Object *> objects;
+    std::vector<PBDSimulationObject *> PBDsims;
+    std::vector<SPHSimulationObject *> SPHsims;
+    Camera camera;
+    std::vector<Light> lights;
+
+    InputHandler *inputHandler = InputHandler::GetInstance();
+
+    bool gravityOn = true;
+
+public:
+    Scene(int w, int h);
+
+    void Build();
+
+    void Update(float delta_t);
+
+    void Render();
+
+    void ResetExternalForces();
+
+    void addForce(vec3 vec3);
+
+    void HandleKeyPress();
+
+    void HandleCameraMove();
+
+    void TranslateCamera(vec3 t);
+};
+
+
+#endif //BRAVE2_SCENE_H
diff --git a/src/rendering/camera.cpp b/src/rendering/camera.cpp
deleted file mode 100644
index 2ae75b6..0000000
--- a/src/rendering/camera.cpp
+++ /dev/null
@@ -1,101 +0,0 @@
-#include <iostream>
-#include "camera.h"
-#include "../utility/custom_math.h"
-
-Camera::Camera() :
-        aspect_ratio(1),
-        pos(1, 1, -3),
-        forward(0, 0, 1),
-        right(-1, 0, 0),
-        speed(10) {
-    inputHandler = InputHandler::GetInstance();
-}
-
-Camera::~Camera() {
-}
-
-
-void Camera::control(float delta_time) {
-    // process camera keys
-    if (InputHandler::IsPressed(GLFW_KEY_W))
-        pos += forward * speed * delta_time;
-
-    if (InputHandler::IsPressed(GLFW_KEY_S))
-        pos -= forward * speed * delta_time;
-
-    if (InputHandler::IsPressed(GLFW_KEY_A))
-        pos -= right * speed * delta_time;
-
-    if (InputHandler::IsPressed(GLFW_KEY_D))
-        pos += right * speed * delta_time;
-
-    if (InputHandler::IsPressed(GLFW_KEY_Q))
-        pos += glm::vec3(0, -1, 0) * speed * delta_time;
-
-    if (InputHandler::IsPressed(GLFW_KEY_E))
-        pos += glm::vec3(0, 1, 0) * speed * delta_time;
-
-    std::cout << pos.x << pos.y << pos.z << std::endl;
-}
-
-void Camera::startDrag(int x, int y) {
-    drag_start = glm::vec2(0, 0);
-}
-
-void Camera::drag(int x, int y) {
-    glm::vec2 d(drag_start.x - (float) x, drag_start.y - (float) y);
-
-    forward = glm::normalize(forward);
-    right = cross(forward, glm::vec3(0, 1, 0));
-    right = glm::normalize(right);
-
-    // TODO 
-}
-
-
-void Camera::glSetupCamera() const {
-    glMatrixMode(GL_PROJECTION);
-    glLoadIdentity();
-    gluPerspective(45.0f, aspect_ratio, 1.0f, 1000.0f);
-
-    glMatrixMode(GL_MODELVIEW);
-    glLoadIdentity();
-    gluLookAt(
-            pos.x, pos.y, pos.z,
-            pos.x + forward.x, pos.y + forward.y, pos.z + forward.z,
-            0.0f, 1.0f, 0.0f);
-}
-
-glm::mat4 Camera::V() const {
-    glm::vec3 w = normalize(forward);
-    glm::vec3 u = normalize(right);
-    glm::vec3 v = cross(w, u);
-    return translate_matrix(-pos) * glm::mat4(
-            u.x, v.x, w.x, 0,
-            u.y, v.y, w.y, 0,
-            u.z, v.z, w.z, 0,
-            0, 0, 0, 1
-    );
-}
-
-glm::mat4 Camera::P() const {
-    float fov = 1.38f;
-    float asp = aspect_ratio;
-    float fp = 0.2f;
-    float bp = 10;
-
-    return mat4(1.0f / (tanf(fov / 2) * asp), 0, 0, 0,
-                0, 1.0f / tanf(fov / 2), 0, 0,
-                0, 0, -(fp + bp) / (bp - fp), -1,
-                0, 0, -2 * fp * bp / (bp - fp), 0);
-}
-
-RenderState Camera::getState() {
-    RenderState state;
-    state.wEye = pos;
-    state.V = V();
-    state.P = P();
-
-    return state;
-}
-
diff --git a/src/rendering/camera.h b/src/rendering/camera.h
deleted file mode 100644
index 4db6e79..0000000
--- a/src/rendering/camera.h
+++ /dev/null
@@ -1,45 +0,0 @@
-#ifndef MADID_CAMERA_H
-#define MADID_CAMERA_H
-
-#include <glm/glm.hpp>
-#include "../utility/gl.h"
-#include "../utility/input_handler.h"
-#include "render_state.h"
-
-class Camera {
-    float aspect_ratio;
-    glm::vec3 pos;
-    glm::vec3 forward;
-    glm::vec3 right;
-    float speed;
-
-    InputHandler *inputHandler;
-
-    glm::vec2 drag_start = glm::vec2(0, 0);
-public:
-    void setAspectRatio(float ar) {
-        aspect_ratio = ar;
-    }
-
-    void control(float delta_time);
-
-    void startDrag(int x, int y);
-
-    void drag(int x, int y);
-
-    void glSetupCamera() const;
-
-    Camera();
-
-    ~Camera();
-
-    // View Matrix: translates the center to the origin
-    glm::mat4 V() const;
-
-    // Projection Matrix
-    glm::mat4 P() const;
-
-    RenderState getState();
-};
-
-#endif
diff --git a/src/rendering/geometry.h b/src/rendering/geometry.h
deleted file mode 100644
index 5915f15..0000000
--- a/src/rendering/geometry.h
+++ /dev/null
@@ -1,25 +0,0 @@
-#ifndef MADID_GEOMETRY_H
-#define MADID_GEOMETRY_H
-
-#include <glm/glm.hpp>
-
-#include <GL/glew.h>
-#include <vector>
-#include "vertex_data.h"
-
-class Geometry {
-protected:
-    GLuint vao, vbo;
-
-public:
-    std::vector<VertexData> vtxData;    // vertices on the CPU
-
-    Geometry();
-
-    virtual void render() = 0;
-
-    ~Geometry();
-};
-
-
-#endif //MADID_GEOMETRY_H
diff --git a/src/rendering/lights/Light.cpp b/src/rendering/lights/Light.cpp
new file mode 100644
index 0000000..253aff4
--- /dev/null
+++ b/src/rendering/lights/Light.cpp
@@ -0,0 +1,5 @@
+#include "Light.h"
+
+void Light::Animate(float t) {
+    //stay in place
+}
diff --git a/src/rendering/lights/Light.h b/src/rendering/lights/Light.h
new file mode 100644
index 0000000..e9e081b
--- /dev/null
+++ b/src/rendering/lights/Light.h
@@ -0,0 +1,13 @@
+#ifndef BRAVE2_LIGHT_H
+#define BRAVE2_LIGHT_H
+
+#include "../../utils/math.h"
+
+struct Light {
+    vec3 La, Le;
+    vec4 wLightPos; // homogenous coordinates, can be at ideal point
+
+    void Animate(float t);
+};
+
+#endif //BRAVE2_LIGHT_H
diff --git a/src/rendering/lights/light.cpp b/src/rendering/lights/light.cpp
deleted file mode 100644
index 1b6d9f6..0000000
--- a/src/rendering/lights/light.cpp
+++ /dev/null
@@ -1,6 +0,0 @@
-#include "light.h"
-
-
-void Light::Animate(float t) {
-    // stay in place
-}
diff --git a/src/rendering/lights/light.h b/src/rendering/lights/light.h
deleted file mode 100644
index 9c4fda6..0000000
--- a/src/rendering/lights/light.h
+++ /dev/null
@@ -1,18 +0,0 @@
-#ifndef MADID_LIGHT_H
-#define MADID_LIGHT_H
-
-#include <glm/glm.hpp>
-
-class Light {
-    // homogenous coordinates, can be at ideal point
-
-    void Animate(float t);
-
-public:
-    glm::vec3 La;
-    glm::vec3 Le;
-    glm::vec4 wLightPos;
-};
-
-
-#endif //MADID_LIGHT_H
diff --git a/src/rendering/materials/Material.cpp b/src/rendering/materials/Material.cpp
new file mode 100644
index 0000000..c2f8f9b
--- /dev/null
+++ b/src/rendering/materials/Material.cpp
@@ -0,0 +1 @@
+#include "Material.h"
diff --git a/src/rendering/materials/Material.h b/src/rendering/materials/Material.h
new file mode 100644
index 0000000..52d77ec
--- /dev/null
+++ b/src/rendering/materials/Material.h
@@ -0,0 +1,12 @@
+#ifndef BRAVE2_MATERIAL_H
+#define BRAVE2_MATERIAL_H
+
+
+#include "../../utils/math.h"
+
+struct Material {
+    vec3 kd, ks, ka;
+    float shininess;
+};
+
+#endif //BRAVE2_MATERIAL_H
diff --git a/src/rendering/materials/material.h b/src/rendering/materials/material.h
deleted file mode 100644
index ab82179..0000000
--- a/src/rendering/materials/material.h
+++ /dev/null
@@ -1,12 +0,0 @@
-#ifndef MADID_MATERIAL_H
-#define MADID_MATERIAL_H
-
-#include <glm/glm.hpp>
-
-struct Material {
-    glm::vec3 kd, ks, ka;
-    float shininess;
-};
-
-
-#endif //MADID_MATERIAL_H
diff --git a/src/rendering/objects/simulation_object.cpp b/src/rendering/objects/simulation_object.cpp
deleted file mode 100644
index 1cb3769..0000000
--- a/src/rendering/objects/simulation_object.cpp
+++ /dev/null
@@ -1,20 +0,0 @@
-#include "simulation_object.h"
-
-SimulationObject::SimulationObject(Shader *_shader, PBDSimulation *_sim):
-    Object(_shader, _sim) {
-
-}
-
-void SimulationObject::ResetExternalForces() {
-
-}
-
-void SimulationObject::AddForce(glm::vec3 f) {
-    reinterpret_cast<Simulation *>(geometry)->addForce(f);
-}
-
-void SimulationObject::update(float delta_t) {
-    reinterpret_cast<Simulation *>(geometry)->update(delta_t);
-}
-
-
diff --git a/src/rendering/objects/simulation_object.h b/src/rendering/objects/simulation_object.h
deleted file mode 100644
index fe1489c..0000000
--- a/src/rendering/objects/simulation_object.h
+++ /dev/null
@@ -1,19 +0,0 @@
-#ifndef MADID_SIMULATION_OBJECT_H
-#define MADID_SIMULATION_OBJECT_H
-
-#include "object.h"
-#include "../../simulation/pbd/pbd_simulation.h"
-
-class SimulationObject : public Object {
-public:
-    SimulationObject(Shader *_shader, PBDSimulation *_sim);
-
-    void ResetExternalForces();
-
-    void AddForce(glm::vec3 f);
-
-    void update(float delta_t) override;
-};
-
-
-#endif //MADID_SIMULATION_OBJECT_H
diff --git a/src/rendering/render_state.h b/src/rendering/render_state.h
deleted file mode 100644
index 4a4649a..0000000
--- a/src/rendering/render_state.h
+++ /dev/null
@@ -1,20 +0,0 @@
-#ifndef MADID_RENDER_STATE_H
-#define MADID_RENDER_STATE_H
-
-#include "materials/material.h"
-#include "lights/light.h"
-#include "textures/texture.h"
-
-using namespace glm;
-
-struct RenderState {
-    mat4 MVP, M, Minv, V, P;
-    vec3 wEye;
-    Texture *texture;
-    Material *material;
-
-    std::vector<Light> lights;
-};
-
-
-#endif //MADID_RENDER_STATE_H
diff --git a/src/rendering/shaders/basic_shader.cpp b/src/rendering/shaders/BasicShader.cpp
similarity index 73%
rename from src/rendering/shaders/basic_shader.cpp
rename to src/rendering/shaders/BasicShader.cpp
index 505178b..293b3ac 100644
--- a/src/rendering/shaders/basic_shader.cpp
+++ b/src/rendering/shaders/BasicShader.cpp
@@ -1,4 +1,8 @@
-#include "basic_shader.h"
+//
+// Created by bobarna on 2020. 11. 28..
+//
+
+#include "BasicShader.h"
 
 BasicShader::BasicShader() {
     create(vertexCode, fragmentCode, "fragmentColor");
@@ -7,4 +11,4 @@ BasicShader::BasicShader() {
 void BasicShader::Bind(RenderState state) {
     Use(); //make this program run
     setUniform(state.MVP, "MVP");
-}
\ No newline at end of file
+}
diff --git a/src/rendering/shaders/basic_shader.h b/src/rendering/shaders/BasicShader.h
similarity index 90%
rename from src/rendering/shaders/basic_shader.h
rename to src/rendering/shaders/BasicShader.h
index 959ab86..5eab910 100644
--- a/src/rendering/shaders/basic_shader.h
+++ b/src/rendering/shaders/BasicShader.h
@@ -1,7 +1,7 @@
-#ifndef MADID_BASIC_SHADER_H
-#define MADID_BASIC_SHADER_H
+#ifndef BRAVE2_BASICSHADER_H
+#define BRAVE2_BASICSHADER_H
 
-#include "shader.h"
+#include "Shader.h"
 
 class BasicShader : public Shader {
     const char *vertexCode = R"(#version 330 core
@@ -35,4 +35,4 @@ public:
 };
 
 
-#endif //MADID_BASIC_SHADER_H
+#endif //BRAVE2_BASICSHADER_H
diff --git a/src/rendering/shaders/PhongShader.cpp b/src/rendering/shaders/PhongShader.cpp
new file mode 100644
index 0000000..f5e549b
--- /dev/null
+++ b/src/rendering/shaders/PhongShader.cpp
@@ -0,0 +1,20 @@
+#include "PhongShader.h"
+
+void PhongShader::Bind(RenderState state) {
+    Use();        // make this program run
+    setUniform(state.MVP, "MVP");
+    setUniform(state.M, "M");
+    setUniform(state.Minv, "Minv");
+    setUniform(state.wEye, "wEye");
+    setUniform(*state.texture, std::string("diffuseTexture"));
+    setUniformMaterial(*state.material, "materials");
+
+    setUniform((int) state.lights.size(), "nLights");
+    for (unsigned int i = 0; i < state.lights.size(); i++) {
+        setUniformLight(state.lights[i], std::string("lights[") + std::to_string(i) + std::string("]"));
+    }
+}
+
+PhongShader::PhongShader() {
+    create(vertexSource, fragmentSource, "fragmentColor");
+}
diff --git a/src/rendering/shaders/PhongShader.h b/src/rendering/shaders/PhongShader.h
new file mode 100644
index 0000000..e776852
--- /dev/null
+++ b/src/rendering/shaders/PhongShader.h
@@ -0,0 +1,99 @@
+#ifndef BRAVE2_PHONGSHADER_H
+#define BRAVE2_PHONGSHADER_H
+
+
+#include "Shader.h"
+
+class PhongShader : public Shader {
+    const char *vertexSource = R"(
+		#version 330
+		precision highp float;
+
+		struct Light {
+			vec3 La, Le;
+			vec4 wLightPos;
+		};
+
+		uniform mat4  MVP, M, Minv; // MVP, Model, Model-inverse
+		uniform Light[8] lights;    // light sources
+		uniform int   nLights;
+		uniform vec3  wEye;         // pos of eye
+
+		layout(location = 0) in vec3  vtxPos;            // pos in modeling space
+		layout(location = 1) in vec3  vtxNorm;      	 // normal in modeling space
+		layout(location = 2) in vec2  vtxUV;
+
+		out vec3 wNormal;		    // normal in world space
+		out vec3 wView;             // view in world space
+		out vec3 wLight[8];		    // light dir in world space
+		out vec2 texcoord;
+
+		void main() {
+			gl_Position = vec4(vtxPos, 1) * MVP; // to NDC
+			// vectors for radiance computation
+			vec4 wPos = vec4(vtxPos, 1) * M;
+			for(int i = 0; i < nLights; i++) {
+				wLight[i] = lights[i].wLightPos.xyz * wPos.w - wPos.xyz * lights[i].wLightPos.w;
+			}
+		    wView  = wEye * wPos.w - wPos.xyz;
+		    wNormal = (Minv * vec4(vtxNorm, 0)).xyz;
+		    texcoord = vtxUV;
+		}
+	)";
+
+    // fragment shader in GLSL
+    const char *fragmentSource = R"(
+		#version 330
+		precision highp float;
+
+		struct Light {
+			vec3 La, Le;
+			vec4 wLightPos;
+		};
+
+		struct Material {
+			vec3 kd, ks, ka;
+			float shininess;
+		};
+
+		uniform Material materials;
+		uniform Light[8] lights;    // light sources
+		uniform int   nLights;
+		uniform sampler2D diffuseTexture;
+
+		in  vec3 wNormal;       // interpolated world sp normal
+		in  vec3 wView;         // interpolated world sp view
+		in  vec3 wLight[8];     // interpolated world sp illum dir
+		in  vec2 texcoord;
+
+        out vec4 fragmentColor; // output goes to frame buffer
+
+		void main() {
+			vec3 N = normalize(wNormal);
+			vec3 V = normalize(wView);
+//			if (dot(N, V) < 0) N = -N;	// prepare for one-sided surfaces like Mobius or Klein
+			vec3 texColor = texture(diffuseTexture, texcoord).rgb;
+			vec3 ka = materials.ka * texColor;
+			vec3 kd = materials.kd * texColor;
+
+			vec3 radiance = vec3(0, 0, 0);
+			for(int i = 0; i < nLights; i++) {
+				vec3 L = normalize(wLight[i]);
+				vec3 H = normalize(L + V);
+				float cost = max(dot(N,L), 0), cosd = max(dot(N,H), 0);
+				// kd and ka are modulated by the texture
+				radiance += ka * lights[i].La +
+                           (kd * texColor * cost + materials.ks * pow(cosd, materials.shininess)) * lights[i].Le;
+			}
+			fragmentColor = vec4(radiance, 1);
+//            fragmentColor = vec4(N, 1);
+		}
+	)";
+public:
+    PhongShader();
+
+    void Bind(RenderState state);
+};
+
+
+#endif //BRAVE2_PHONGSHADER_H
diff --git a/src/rendering/shaders/shader.cpp b/src/rendering/shaders/Shader.cpp
similarity index 95%
rename from src/rendering/shaders/shader.cpp
rename to src/rendering/shaders/Shader.cpp
index fe56cda..579bbfe 100644
--- a/src/rendering/shaders/shader.cpp
+++ b/src/rendering/shaders/Shader.cpp
@@ -1,4 +1,4 @@
-#include "shader.h"
+#include "Shader.h"
 
 void Shader::setUniformMaterial(Material &material, const std::string &name) {
     setUniform(material.kd, name + ".kd");
@@ -11,4 +11,4 @@ void Shader::setUniformLight(const Light &light, const std::string &name) {
     setUniform(light.La, name + ".La");
     setUniform(light.Le, name + ".Le");
     setUniform(light.wLightPos, name + ".wLightPos");
-}
+}
\ No newline at end of file
diff --git a/src/rendering/shaders/Shader.h b/src/rendering/shaders/Shader.h
new file mode 100644
index 0000000..b668b3a
--- /dev/null
+++ b/src/rendering/shaders/Shader.h
@@ -0,0 +1,25 @@
+#ifndef BRAVE2_SHADER_H
+#define BRAVE2_SHADER_H
+
+#include <string>
+#include <GL/glew.h>
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include "../../utils/math.h"
+#include "../RenderState.h"
+#include "../../utils/GPUProgram.h"
+#include "../materials/Material.h"
+#include "../lights/Light.h"
+
+struct Shader : public GPUProgram {
+public:
+    virtual void Bind(RenderState state) = 0;
+
+    void setUniformMaterial(Material &material, const std::string &name);
+
+    void setUniformLight(const Light &light, const std::string &name);
+};
+
+
+#endif //BRAVE2_SHADER_H
diff --git a/src/rendering/shaders/gpu_program.cpp b/src/rendering/shaders/gpu_program.cpp
deleted file mode 100644
index 9eaff0f..0000000
--- a/src/rendering/shaders/gpu_program.cpp
+++ /dev/null
@@ -1 +0,0 @@
-#include "gpu_program.h"
diff --git a/src/rendering/shaders/shader.h b/src/rendering/shaders/shader.h
deleted file mode 100644
index 2cc8c83..0000000
--- a/src/rendering/shaders/shader.h
+++ /dev/null
@@ -1,19 +0,0 @@
-#ifndef MADID_SHADER_H
-#define MADID_SHADER_H
-
-#include "gpu_program.h"
-#include "../render_state.h"
-#include "../materials/material.h"
-#include "../lights/light.h"
-
-class Shader : public GPUProgram {
-public:
-    virtual void Bind(RenderState state) = 0;
-
-    void setUniformMaterial(Material &material, const std::string &name);
-
-    void setUniformLight(const Light &light, const std::string &name);
-};
-
-
-#endif //MADID_SHADER_H
diff --git a/src/rendering/textures/Texture.cpp b/src/rendering/textures/Texture.cpp
new file mode 100644
index 0000000..c741ecd
--- /dev/null
+++ b/src/rendering/textures/Texture.cpp
@@ -0,0 +1 @@
+#include "Texture.h"
diff --git a/src/rendering/textures/texture.h b/src/rendering/textures/Texture.h
similarity index 73%
rename from src/rendering/textures/texture.h
rename to src/rendering/textures/Texture.h
index 465584f..3a9c29e 100644
--- a/src/rendering/textures/texture.h
+++ b/src/rendering/textures/Texture.h
@@ -1,9 +1,7 @@
-#ifndef MADID_TEXTURE_H
-#define MADID_TEXTURE_H
+#ifndef BRAVE2_TEXTURE_H
+#define BRAVE2_TEXTURE_H
 
-#include "../../utility/custom_math.h"
-#include <iostream>
-#include "../../utility/gl.h"
+#include "../../utils/math.h"
 
 class Texture {
     std::vector<vec4> load(std::string pathname, bool transparent, int &width, int &height) {
@@ -76,5 +74,31 @@ public:
     }
 };
 
+class UniformColorTexture : public Texture {
+public:
+    UniformColorTexture(float r = 1, float g = 1, float b = 0, float a = 1) : Texture() {
+        std::vector<vec4> color = {r, g, b, a};
+
+        create(1, 1, color, GL_NEAREST);
+    }
+};
+
+//---------------------------
+class CheckerBoardTexture : public Texture {
+//---------------------------
+public:
+    CheckerBoardTexture(const int width = 0, const int height = 0,
+                        float c1r = 1, float c1g = 1, float c1b = 0, float c1a = 1,
+                        float c2r = 0, float c2g = 0, float c2b = 1, float c2a = 1) : Texture() {
+        std::vector<vec4> image(width * height);
+        const vec4 c1(c1r, c1g, c1b, c1a), c2(c2r, c2g, c2b, c2a);
+
+        for (int x = 0; x < width; x++)
+            for (int y = 0; y < height; y++)
+                image[y * width + x] = (x & 1) ^ (y & 1) ? c1 : c2;
+
+        create(width, height, image, GL_NEAREST);
+    }
+};
 
-#endif //MADID_TEXTURE_H
+#endif //BRAVE2_TEXTURE_H
diff --git a/src/rendering/textures/texture.cpp b/src/rendering/textures/texture.cpp
deleted file mode 100644
index ca7d364..0000000
--- a/src/rendering/textures/texture.cpp
+++ /dev/null
@@ -1,5 +0,0 @@
-//
-// Created by bobarna on 2021. 05. 13..
-//
-
-#include "texture.h"
diff --git a/src/rendering/textures/texture_uniform_color.h b/src/rendering/textures/texture_uniform_color.h
deleted file mode 100644
index 9705c11..0000000
--- a/src/rendering/textures/texture_uniform_color.h
+++ /dev/null
@@ -1,17 +0,0 @@
-#ifndef MADID_TEXTURE_UNIFORM_COLOR_H
-#define MADID_TEXTURE_UNIFORM_COLOR_H
-
-#include "texture.h"
-
-class UniformColorTexture : public Texture {
-public:
-    UniformColorTexture(float r = 1, float g = 1, float b = 0, float a = 1) : Texture() {
-        std::vector<vec4> color;
-        color.emplace_back(r, g, b, a);
-
-        create(1, 1, color, GL_NEAREST);
-    }
-};
-
-
-#endif //MADID_TEXTURE_UNIFORM_COLOR_H
diff --git a/src/rendering/vertex_data.cpp b/src/rendering/vertex_data.cpp
deleted file mode 100644
index 3a43a49..0000000
--- a/src/rendering/vertex_data.cpp
+++ /dev/null
@@ -1,8 +0,0 @@
-#include "vertex_data.h"
-
-VertexData::VertexData(glm::vec3 p, glm::vec3 n, glm::vec2 uv) :
-        position(p), normal(n), texcoord(uv) {
-}
-
-VertexData::VertexData() {
-}
\ No newline at end of file
diff --git a/src/rendering/vertex_data.h b/src/rendering/vertex_data.h
deleted file mode 100644
index a90397b..0000000
--- a/src/rendering/vertex_data.h
+++ /dev/null
@@ -1,17 +0,0 @@
-
-#ifndef MADID_VERTEX_DATA_H
-#define MADID_VERTEX_DATA_H
-
-#include <glm/glm.hpp>
-
-struct VertexData {
-    // needed for GetVertexByUV function
-    bool valid = true;
-    glm::vec3 position, normal;
-    glm::vec2 texcoord;
-    VertexData();
-    VertexData(glm::vec3 p, glm::vec3 n, glm::vec2 uv);
-    explicit VertexData(bool valid, glm::vec3 p = {0,0,0}, glm::vec3 n = {0,0,0}, glm::vec2 uv = {0,0}){};
-};
-
-#endif //MADID_VERTEX_DATA_H
diff --git a/src/scene.cpp b/src/scene.cpp
deleted file mode 100644
index d59df65..0000000
--- a/src/scene.cpp
+++ /dev/null
@@ -1,35 +0,0 @@
-#include "scene.h"
-#include "rendering/shaders/basic_shader.h"
-
-Scene::Scene() {
-    camera->glSetupCamera();
-
-    size_t nrSims = 100;
-    size_t nrSegments = 50;
-    float lSeg = 0.003f;
-
-    auto pbd_sim = new PBDSimulation(nrSims, nrSegments, lSeg);
-
-    Shader *basicShader = new BasicShader();
-
-    auto pbd_sim_object = std::make_shared<SimulationObject>(basicShader, pbd_sim);
-
-    simulations.push_back(pbd_sim_object);
-
-    std::cout << "Scene built successfully." << std::endl;
-}
-
-Scene::~Scene() {
-    std::cout << "Scene destroyed." << std::endl;
-}
-
-void Scene::update(double dt) {
-    camera->control(dt);
-
-    for (auto &s: simulations) s->update(dt);
-}
-
-void Scene::render() {
-    camera->glSetupCamera();
-    for (auto &s: simulations) s->render(camera->getState());
-}
diff --git a/src/scene.h b/src/scene.h
deleted file mode 100644
index 6ed4e1a..0000000
--- a/src/scene.h
+++ /dev/null
@@ -1,32 +0,0 @@
-#ifndef MADID_SCENE_H
-#define MADID_SCENE_H
-
-#include <iostream>
-#include <memory>
-
-#include "rendering/objects/object.h"
-#include "rendering/objects/simulation_object.h"
-#include "rendering/camera.h"
-
-class Scene {
-    bool draw_normals;
-    bool draw_surface;
-    bool draw_velocity;
-    bool use_emitter;
-
-    std::shared_ptr<Camera> camera = std::make_shared<Camera>();
-
-    // TODO shared_ptr
-    std::vector<std::shared_ptr<SimulationObject>> simulations;
-
-public:
-    Scene();
-
-    ~Scene();
-
-    void update(double dt);
-
-    void render();
-};
-
-#endif
diff --git a/src/simulation/particle.cpp b/src/simulation/particle.cpp
deleted file mode 100644
index fc9acd9..0000000
--- a/src/simulation/particle.cpp
+++ /dev/null
@@ -1,25 +0,0 @@
-#include "particle.h"
-
-Particle::Particle(float x, float y, float z) :
-        pos(x, y, z),
-        old_pos(pos),
-        radius(1), //TODO
-        radius_sqr(1 * 1), //TODO
-        v(0, 0, 0),
-        a(0, 0, 0),
-        f(0, 0, 0),
-        n(0, 0, 0),
-        w(0), // inverse mass, 1/inf is 0
-        rho(1),
-        p(0),
-        fixed(false) {
-}
-
-Particle::~Particle() {}
-
-Particle::Particle(Particle::vector pos, float w, Particle::vector color) :
-        pos(pos),
-        w(w),
-        color(color) {
-
-}
diff --git a/src/simulation/particle.h b/src/simulation/particle.h
deleted file mode 100644
index 624e3e1..0000000
--- a/src/simulation/particle.h
+++ /dev/null
@@ -1,35 +0,0 @@
-#ifndef MADID_PARTICLE
-#define MADID_PARTICLE
-
-#include <glm/glm.hpp>
-
-class Particle {
-public:
-    typedef glm::vec3 vector;
-    typedef vector point_type;
-
-    Particle(float x, float y, float z);
-
-    Particle(vector pos, float w, vector color);
-
-    ~Particle();
-
-    //TODO proper scope, maybe getter functions
-public:
-    point_type pos;     // position of the particle
-    point_type old_pos; // old position of the particle
-    float radius;       // radius of the particle
-    float radius_sqr;   // squared radius of the particle
-    vector v;           // current velocity
-    vector a;           // current acceleration
-    vector f;           // current sum of external forces acting on the particle
-    vector n;           // surface normal if close to/on the surface, else 0
-    float w;            // inverse mass (constant)
-    vector tmp_pos;     // temporary position (for PBD simulation)
-    vector color;       // (r,g,b) color of the particle
-    float rho;          // density (varies)
-    float p;            // pressure (varies)
-    bool fixed;         // Is this particled fixed (locked)?
-};
-
-#endif
diff --git a/src/simulation/simulation.h b/src/simulation/simulation.h
deleted file mode 100644
index 9310268..0000000
--- a/src/simulation/simulation.h
+++ /dev/null
@@ -1,19 +0,0 @@
-#ifndef MADID_SIMULATION_H
-#define MADID_SIMULATION_H
-
-
-#include "../rendering/geometry.h"
-
-class Simulation : public Geometry {
-public:
-    virtual void addForce(glm::vec3 force) = 0;
-    Simulation(): Geometry() { };
-
-    virtual void update(float dt) = 0;
-    virtual void render() = 0;
-
-    virtual void resetExternalForces() = 0;
-};
-
-
-#endif //MADID_SIMULATION_H
diff --git a/src/utility/custom_math.h b/src/utility/custom_math.h
deleted file mode 100644
index 31dbf91..0000000
--- a/src/utility/custom_math.h
+++ /dev/null
@@ -1,69 +0,0 @@
-#ifndef MADID_CUSTOM_MATH_H
-#define MADID_CUSTOM_MATH_H
-
-#include <glm/glm.hpp>
-#include <cmath>
-#include <vector>
-
-using namespace glm;
-
-inline mat4 translate_matrix(vec3 t) {
-    return mat4(vec4(1, 0, 0, 0),
-                vec4(0, 1, 0, 0),
-                vec4(0, 0, 1, 0),
-                vec4(t.x, t.y, t.z, 1));
-}
-
-inline mat4 scale_matrix(vec3 s) {
-    return mat4(vec4(s.x, 0, 0, 0),
-                vec4(0, s.y, 0, 0),
-                vec4(0, 0, s.z, 0),
-                vec4(0, 0, 0, 1));
-}
-
-inline mat4 rotation_matrix(float angle, vec3 w) {
-    float c = cosf(angle), s = sinf(angle);
-    w = normalize(w);
-    return mat4(vec4(c * (1 - w.x * w.x) + w.x * w.x, w.x * w.y * (1 - c) + w.z * s, w.x * w.z * (1 - c) - w.y * s, 0),
-                vec4(w.x * w.y * (1 - c) - w.z * s, c * (1 - w.y * w.y) + w.y * w.y, w.y * w.z * (1 - c) + w.x * s, 0),
-                vec4(w.x * w.z * (1 - c) + w.y * s, w.y * w.z * (1 - c) - w.x * s, c * (1 - w.z * w.z) + w.z * w.z, 0),
-                vec4(0, 0, 0, 1));
-}
-
-inline mat4 x_rotation_matrix(float angle) {
-    float s = sinf(angle);
-    float c = cosf(angle);
-    return mat4(
-            1.f, 0.f, 0.f, 0.f,
-            0.f, c, s, 0.f,
-            0.f, -s, c, 0.f,
-            0.f, 0.f, 0.f, 1.f
-    );
-}
-
-inline mat4 y_rotation_matrix(float angle) {
-    float s = sinf(angle);
-    float c = cosf(angle);
-    return mat4(
-            c, 0.f, -s, 0.f,
-            0.f, 1.f, 0.f, 0.f,
-            s, 0.f, c, 0.f,
-            0.f, 0.f, 0.f, 1.f
-    );
-}
-
-inline mat4 z_rotation_matrix(float angle) {
-    float s = sinf(angle);
-    float c = cosf(angle);
-    return mat4(
-            c, s, 0.f, 0.f,
-            -s, c, 0.f, 0.f,
-            0.f, 0.f, 1.f, 0.f,
-            0.f, 0.f, 0.f, 1.f
-    );
-}
-
-// automatic mat4 printing
-//std::ostream &operator<<(std::ostream &out, const mat4 &v);
-
-#endif //MADID_CUSTOM_MATH_H
diff --git a/src/utility/gl.h b/src/utility/gl.h
deleted file mode 100644
index dab3244..0000000
--- a/src/utility/gl.h
+++ /dev/null
@@ -1,7 +0,0 @@
-#ifndef MADID_GL_H
-#define MADID_GL_H
-
-#include <GL/glew.h>
-#include <GLFW/glfw3.h>
-
-#endif
diff --git a/src/rendering/shaders/gpu_program.h b/src/utils/GPUProgram.h
similarity index 96%
rename from src/rendering/shaders/gpu_program.h
rename to src/utils/GPUProgram.h
index d341e91..75f4b5a 100644
--- a/src/rendering/shaders/gpu_program.h
+++ b/src/utils/GPUProgram.h
@@ -1,12 +1,13 @@
-#ifndef MADID_GPU_PROGRAM_H
-#define MADID_GPU_PROGRAM_H
+#ifndef BRAVE2_GPUPROGRAM_H
+#define BRAVE2_GPUPROGRAM_H
 
-#include <glm/gtc/type_ptr.hpp>
-#include "../../utility/gl.h"
-#include "../../utility/custom_math.h"
-#include "../textures/texture.h"
+#include <GL/glew.h>
+#include "math.h"
+#include "../rendering/textures/Texture.h"
 
+//---------------------------
 class GPUProgram {
+//--------------------------
     unsigned int shaderProgramId = 0;
     unsigned int vertexShader = 0, geometryShader = 0, fragmentShader = 0;
     bool waitError = true;
@@ -155,7 +156,7 @@ public:
 
     void setUniform(const mat4 &mat, const std::string &name) {
         int location = getLocation(name);
-        if (location >= 0) glUniformMatrix4fv(location, 1, GL_TRUE, glm::value_ptr(mat));
+        if (location >= 0) glUniformMatrix4fv(location, 1, GL_TRUE, mat);
     }
 
     void setUniform(const Texture &texture, const std::string &samplerName, unsigned int textureUnit = 0) {
@@ -171,5 +172,4 @@ public:
     ~GPUProgram() { if (shaderProgramId > 0) glDeleteProgram(shaderProgramId); }
 };
 
-
-#endif //MADID_GPU_PROGRAM_H
+#endif //BRAVE2_GPUPROGRAM_H
diff --git a/src/utility/input_handler.cpp b/src/utils/InputHandler.cpp
similarity index 82%
rename from src/utility/input_handler.cpp
rename to src/utils/InputHandler.cpp
index d16f357..64f5f62 100644
--- a/src/utility/input_handler.cpp
+++ b/src/utils/InputHandler.cpp
@@ -1,12 +1,11 @@
-#include "input_handler.h"
+#include "InputHandler.h"
 
-InputHandler *InputHandler::singleton_ = nullptr;
-bool InputHandler::keyPressed[348];
-int InputHandler::modifiers = 0;
+InputHandler *InputHandler::singleton_ = nullptr;;
 
 InputHandler *InputHandler::GetInstance() {
-    if (singleton_ == nullptr)
+    if (singleton_ == nullptr) {
         singleton_ = new InputHandler();
+    }
 
     return singleton_;
 }
@@ -49,5 +48,4 @@ bool InputHandler::IsCapsLock() {
 
 bool InputHandler::IsNumLock() {
     return modifiers & GLFW_MOD_NUM_LOCK;
-}
-
+}
\ No newline at end of file
diff --git a/src/utility/input_handler.h b/src/utils/InputHandler.h
similarity index 63%
rename from src/utility/input_handler.h
rename to src/utils/InputHandler.h
index b3c7493..22995f0 100644
--- a/src/utility/input_handler.h
+++ b/src/utils/InputHandler.h
@@ -1,7 +1,8 @@
-#ifndef MADID_INPUT_HANDLER_H
-#define MADID_INPUT_HANDLER_H
+#ifndef BRAVE2_INPUTHANDLER_H
+#define BRAVE2_INPUTHANDLER_H
 
-#include "gl.h"
+#include <cstddef>
+#include <GLFW/glfw3.h>
 
 /**
  * The inputHandler class defines the `GetInstance` method that serves as an
@@ -11,8 +12,8 @@
 class InputHandler {
     // Input handling
     bool aKeyWasPressed;
-    static bool keyPressed[348];
-    static int modifiers;
+    bool keyPressed[348] = {false};
+    int modifiers;
 
     /**
      * The inputHandler's constructor should always be private to prevent direct
@@ -23,7 +24,6 @@ class InputHandler {
         modifiers = 0;
         for (bool &k : keyPressed)
             k = false;
-
     }
 
     static InputHandler *singleton_;
@@ -42,26 +42,26 @@ public:
 
     static InputHandler *GetInstance();
 
-    static void KeyPress(int key);
+    void KeyPress(int key);
 
-    static void SetModifiers(int m);
+    void SetModifiers(int m);
 
-    static void KeyRelease(int key);
+    void KeyRelease(int key);
 
-    static bool IsPressed(int key);
+    bool IsPressed(int key);
 
-    static bool IsShiftPressed();
+    bool IsShiftPressed();
 
-    static bool IsControlPressed();
+    bool IsControlPressed();
 
-    static bool IsAltPressed();
+    bool IsAltPressed();
 
-    static bool IsSuperPressed();
+    bool IsSuperPressed();
 
-    static bool IsCapsLock();
+    bool IsCapsLock();
 
-    static bool IsNumLock();
+    bool IsNumLock();
 
 };
 
-#endif //MADID_INPUT_HANDLER_H
+#endif //BRAVE2_INPUTHANDLER_H
diff --git a/src/utils/OBJReader.cpp b/src/utils/OBJReader.cpp
new file mode 100644
index 0000000..d2b41ed
--- /dev/null
+++ b/src/utils/OBJReader.cpp
@@ -0,0 +1,89 @@
+#include "OBJReader.h"
+
+OBJReader::OBJReader(const std::string &filePath, std::vector<VertexData> &out_vtxData) : is(filePath) {
+
+    if (!is) {
+        std::cerr << "Can't open OBJ file: " << filePath << std::endl;
+        return;
+    }
+
+    readData();
+
+    for (auto curr: out_vtxData)
+        std::cout << curr << std::endl;
+    for (unsigned int i : vertexIndices)
+        out_vtxData.emplace_back(temp_vertices[i - 1], temp_normals[i - 1], temp_uvs[i - 1]);
+
+}
+
+void OBJReader::readData() {
+    std::string lineType;
+
+    // read the file word by word
+    // the >> operator will only extract what can be considered a word from the stream,
+    // using whitespaces as separators
+    while (is >> lineType) {
+        if (lineType == "#")
+            // Comment, skip the line
+            std::getline(is, lineType);
+        else if (lineType == "v")
+            readVertexLine();
+        else if (lineType == "vt")
+            readUVLine();
+        else if (lineType == "vn")
+            readNormalLine();
+        else if (lineType == "f")
+            readFaceLine();
+        else {
+            // ignore everything else
+            std::cerr << "ignored line: " << lineType;
+            // get rest of the line
+            std::getline(is, lineType);
+            // output rest of the line
+            std::cerr << " " << lineType << std::endl;
+        }
+    }
+}
+
+void OBJReader::readVertexLine() {
+    // Vertex
+    vec3 v;
+    is >> v.x >> v.y >> v.z;
+    temp_vertices.push_back(v);
+}
+
+void OBJReader::readUVLine() {
+    // UV coordinates
+    vec2 uv;
+    is >> uv.x >> uv.y;
+    temp_uvs.push_back(uv);
+}
+
+void OBJReader::readNormalLine() {
+    // Normal vector
+    vec3 n;
+    is >> n.x >> n.y >> n.z;
+    temp_normals.push_back(n);
+}
+
+void OBJReader::readFaceLine() {
+// Face
+    std::string v1, v2, v3;
+    unsigned int v[3], uv[3], n[3];
+    is >> v[0] >> chlit('/') >> uv[0] >> chlit('/') >> n[0];
+    is >> v[1] >> chlit('/') >> uv[1] >> chlit('/') >> n[1];
+    is >> v[2] >> chlit('/') >> uv[2] >> chlit('/') >> n[2];
+    if (is.fail()) {
+        std::cerr << "can't read face data from OBJ file" << std::endl;
+        throw std::exception();
+    }
+    vertexIndices.push_back(v[0]);
+    vertexIndices.push_back(v[1]);
+    vertexIndices.push_back(v[2]);
+    uvIndices.push_back(uv[0]);
+    uvIndices.push_back(uv[1]);
+    uvIndices.push_back(uv[2]);
+    normalIndices.push_back(n[0]);
+    normalIndices.push_back(n[1]);
+    normalIndices.push_back(n[2]);
+}
diff --git a/src/utils/OBJReader.h b/src/utils/OBJReader.h
new file mode 100644
index 0000000..82968f9
--- /dev/null
+++ b/src/utils/OBJReader.h
@@ -0,0 +1,46 @@
+#ifndef BRAVE2_OBJREADER_H
+#define BRAVE2_OBJREADER_H
+
+#include <string>
+#include <vector>
+#include <fstream>
+#include "math.h"
+#include "../geometries/VertexData.h"
+
+class OBJReader {
+    std::vector<unsigned int> vertexIndices, normalIndices, uvIndices;
+    std::vector<vec3> temp_vertices;
+    std::vector<vec3> temp_normals;
+    std::vector<vec2> temp_uvs;
+
+    std::ifstream is;
+
+public:
+    OBJReader(const std::string &filePath, std::vector<VertexData>& out_vtxData);
+
+    void readData();
+
+    void readVertexLine();
+
+    void readUVLine();
+
+    void readNormalLine();
+
+    void readFaceLine();
+};
+
+struct chlit {
+    char c_target;
+
+    chlit(char _c) : c_target(_c) {}
+
+    friend std::istream &operator>>(std::istream &is, chlit x) {
+        char c;
+        if (is >> c && c != x.c_target)
+            is.setstate(std::iostream::failbit);
+        return is;
+    }
+
+};
+
+#endif //BRAVE2_OBJREADER_H
diff --git a/src/utils/math.cpp b/src/utils/math.cpp
new file mode 100644
index 0000000..6bd4668
--- /dev/null
+++ b/src/utils/math.cpp
@@ -0,0 +1,20 @@
+#include "math.h"
+
+std::ostream &operator<<(std::ostream &out, const vec2 &v) {
+    return out << "(" << v.x << ", " << v.y << ")";
+}
+
+std::ostream &operator<<(std::ostream &out, const vec3 &v) {
+    return out << "(" << v.x << ", " << v.y << ", " << v.z << ")";
+}
+
+std::ostream &operator<<(std::ostream &out, const vec4 &v) {
+    return out << "(" << v.x << ", " << v.y << ", " << v.z << ", " << v.w << ")";
+}
+
+std::ostream &operator<<(std::ostream &out, const mat4 &m) {
+    return out << m.rows[0] << std::endl <<
+               m.rows[1] << std::endl <<
+               m.rows[2] << std::endl <<
+               m.rows[3] << std::endl;
+}
\ No newline at end of file
diff --git a/src/utils/math.h b/src/utils/math.h
new file mode 100644
index 0000000..3345f6a
--- /dev/null
+++ b/src/utils/math.h
@@ -0,0 +1,378 @@
+#ifndef BRAVE2_MATH_H
+#define BRAVE2_MATH_H
+
+#include <cmath>
+#include <vector>
+#include <string>
+#include <iostream>
+
+#include <GL/glew.h>
+
+//--------------------------
+struct vec2 {
+//--------------------------
+    float x, y;
+
+    vec2(float x0 = 0, float y0 = 0) {
+        x = x0;
+        y = y0;
+    }
+
+    vec2 operator*(float a) const { return vec2(x * a, y * a); }
+
+    vec2 &operator*=(const float &a) {
+        (*this) = (*this) * a;
+        return *this;
+    }
+
+    vec2 operator/(float a) const { return vec2(x / a, y / a); }
+
+    vec2 &operator/=(const float &a) {
+        (*this) = (*this) / a;
+        return *this;
+    }
+
+    vec2 operator+(const vec2 &v) const { return vec2(x + v.x, y + v.y); }
+
+    vec2 &operator+=(const vec2 &rhs) {
+        (*this) = (*this) + rhs;
+        return *this;
+    }
+
+    vec2 operator-(const vec2 &v) const { return vec2(x - v.x, y - v.y); }
+
+    vec2 &operator-=(const vec2 &rhs) {
+        (*this) = (*this) - rhs;
+        return *this;
+    }
+
+    vec2 operator*(const vec2 &v) const { return vec2(x * v.x, y * v.y); }
+
+    vec2 &operator*=(const vec2 &rhs) {
+        (*this) = (*this) * rhs;
+        return *this;
+    }
+};
+
+inline float dot(const vec2 &v1, const vec2 &v2) {
+    return (v1.x * v2.x + v1.y * v2.y);
+}
+
+inline float length(const vec2 &v) { return sqrtf(dot(v, v)); }
+
+inline vec2 normalize(const vec2 &v) { return v * (1 / length(v)); }
+
+inline vec2 operator*(float a, const vec2 &v) { return vec2(v.x * a, v.y * a); }
+
+// automatic vec2 printing
+std::ostream &operator<<(std::ostream &out, const vec2 &v);
+
+//--------------------------
+struct vec3 {
+//--------------------------
+    float x, y, z;
+
+    vec3(float x0 = 0, float y0 = 0, float z0 = 0) {
+        x = x0;
+        y = y0;
+        z = z0;
+    }
+
+    vec3(vec2 v) {
+        x = v.x;
+        y = v.y;
+        z = 0;
+    }
+
+    vec3 operator*(float a) const { return vec3(x * a, y * a, z * a); }
+
+    vec3 operator*=(float const &a) {
+        (*this) = (*this) * a;
+        return *this;
+    }
+
+    vec3 operator/(float a) const { return vec3(x / a, y / a, z / a); }
+
+    vec3 operator/=(float const &a) {
+        (*this) = (*this) / a;
+        return *this;
+    }
+
+    vec3 operator+(const vec3 &v) const { return vec3(x + v.x, y + v.y, z + v.z); }
+
+    vec3 operator+=(const vec3 &v) {
+        (*this) = (*this) + v;
+        return *this;
+    }
+
+    vec3 operator-(const vec3 &v) const { return vec3(x - v.x, y - v.y, z - v.z); }
+
+    vec3 operator-=(const vec3 &v) {
+        (*this) = (*this) - v;
+        return *this;
+    }
+
+    vec3 operator*(const vec3 &v) const { return vec3(x * v.x, y * v.y, z * v.z); }
+
+    vec3 operator*=(const vec3 &v) {
+        (*this) = (*this) * v;
+        return *this;
+    }
+
+    vec3 operator-() const { return vec3(-x, -y, -z); }
+
+};
+
+inline float dot(const vec3 &v1, const vec3 &v2) { return (v1.x * v2.x + v1.y * v2.y + v1.z * v2.z); }
+
+inline float length(const vec3 &v) { return sqrtf(dot(v, v)); }
+
+inline vec3 normalize(const vec3 &v) { return v * (1 / length(v)); }
+
+inline vec3 cross(const vec3 &v1, const vec3 &v2) {
+    return vec3(v1.y * v2.z - v1.z * v2.y, v1.z * v2.x - v1.x * v2.z, v1.x * v2.y - v1.y * v2.x);
+}
+
+inline vec3 operator*(float a, const vec3 &v) { return vec3(v.x * a, v.y * a, v.z * a); }
+
+// automatic vec3 printing
+std::ostream &operator<<(std::ostream &out, const vec3 &v);
+
+//--------------------------
+struct vec4 {
+//--------------------------
+    float x, y, z, w;
+
+    vec4(float x0 = 0, float y0 = 0, float z0 = 0, float w0 = 0) {
+        x = x0;
+        y = y0;
+        z = z0;
+        w = w0;
+    }
+
+    float &operator[](int j) { return *(&x + j); }
+
+    float operator[](int j) const { return *(&x + j); }
+
+    vec4 operator*(float a) const { return vec4(x * a, y * a, z * a, w * a); }
+
+    vec4 operator*=(const float &a) {
+        (*this) = (*this) * a;
+        return *this;
+    }
+
+    vec4 operator/(float d) const { return vec4(x / d, y / d, z / d, w / d); }
+
+    vec4 operator/=(const float &d) {
+        (*this) = (*this) * d;
+        return *this;
+    }
+
+    vec4 operator+(const vec4 &v) const { return vec4(x + v.x, y + v.y, z + v.z, w + v.w); }
+
+    vec4 operator+=(const vec4 &v) {
+        (*this) = (*this) + v;
+        return *this;
+    }
+
+    vec4 operator-(const vec4 &v) const { return vec4(x - v.x, y - v.y, z - v.z, w - v.w); }
+
+    vec4 operator-=(const vec4 &v) {
+        (*this) = (*this) - v;
+        return *this;
+    }
+
+    vec4 operator*(const vec4 &v) const { return vec4(x * v.x, y * v.y, z * v.z, w * v.w); }
+
+    vec4 operator*=(const vec4 &v) {
+        (*this) = (*this) * v;
+        return *this;
+    }
+};
+
+inline float dot(const vec4 &v1, const vec4 &v2) {
+    return (v1.x * v2.x + v1.y * v2.y + v1.z * v2.z + v1.w * v2.w);
+}
+
+inline vec4 operator*(float a, const vec4 &v) {
+    return vec4(v.x * a, v.y * a, v.z * a, v.w * a);
+}
+
+// automatic vec4 printing
+std::ostream &operator<<(std::ostream &out, const vec4 &v);
+
+//---------------------------
+struct mat4 { // row-major matrix 4x4
+//---------------------------
+    vec4 rows[4];
+public:
+    mat4() {}
+
+    mat4(float m00, float m01, float m02, float m03,
+         float m10, float m11, float m12, float m13,
+         float m20, float m21, float m22, float m23,
+         float m30, float m31, float m32, float m33) {
+        rows[0][0] = m00;
+        rows[0][1] = m01;
+        rows[0][2] = m02;
+        rows[0][3] = m03;
+        rows[1][0] = m10;
+        rows[1][1] = m11;
+        rows[1][2] = m12;
+        rows[1][3] = m13;
+        rows[2][0] = m20;
+        rows[2][1] = m21;
+        rows[2][2] = m22;
+        rows[2][3] = m23;
+        rows[3][0] = m30;
+        rows[3][1] = m31;
+        rows[3][2] = m32;
+        rows[3][3] = m33;
+    }
+
+    mat4(vec4 it, vec4 jt, vec4 kt, vec4 ot) {
+        rows[0] = it;
+        rows[1] = jt;
+        rows[2] = kt;
+        rows[3] = ot;
+    }
+
+    vec4 &operator[](int i) { return rows[i]; }
+
+    vec4 operator[](int i) const { return rows[i]; }
+
+    operator float *() const { return (float *) this; }
+};
+
+inline vec4 operator*(const vec4 &v, const mat4 &mat) {
+    return v[0] * mat[0] + v[1] * mat[1] + v[2] * mat[2] + v[3] * mat[3];
+}
+
+inline mat4 operator*(const mat4 &left, const mat4 &right) {
+    mat4 result;
+    for (int i = 0; i < 4; i++) result.rows[i] = left.rows[i] * right;
+    return result;
+}
+
+inline mat4 translate_matrix(vec3 t) {
+    return mat4(vec4(1, 0, 0, 0),
+                vec4(0, 1, 0, 0),
+                vec4(0, 0, 1, 0),
+                vec4(t.x, t.y, t.z, 1));
+}
+
+inline mat4 scale_matrix(vec3 s) {
+    return mat4(vec4(s.x, 0, 0, 0),
+                vec4(0, s.y, 0, 0),
+                vec4(0, 0, s.z, 0),
+                vec4(0, 0, 0, 1));
+}
+
+inline mat4 rotation_matrix(float angle, vec3 w) {
+    float c = cosf(angle), s = sinf(angle);
+    w = normalize(w);
+    return mat4(vec4(c * (1 - w.x * w.x) + w.x * w.x, w.x * w.y * (1 - c) + w.z * s, w.x * w.z * (1 - c) - w.y * s, 0),
+                vec4(w.x * w.y * (1 - c) - w.z * s, c * (1 - w.y * w.y) + w.y * w.y, w.y * w.z * (1 - c) + w.x * s, 0),
+                vec4(w.x * w.z * (1 - c) + w.y * s, w.y * w.z * (1 - c) - w.x * s, c * (1 - w.z * w.z) + w.z * w.z, 0),
+                vec4(0, 0, 0, 1));
+}
+
+inline mat4 x_rotation_matrix(float angle) {
+    float s = sinf(angle);
+    float c = cosf(angle);
+    return mat4(
+            1.f, 0.f, 0.f, 0.f,
+            0.f, c, s, 0.f,
+            0.f, -s, c, 0.f,
+            0.f, 0.f, 0.f, 1.f
+    );
+}
+
+inline mat4 y_rotation_matrix(float angle) {
+    float s = sinf(angle);
+    float c = cosf(angle);
+    return mat4(
+            c, 0.f, -s, 0.f,
+            0.f, 1.f, 0.f, 0.f,
+            s, 0.f, c, 0.f,
+            0.f, 0.f, 0.f, 1.f
+    );
+}
+
+inline mat4 z_rotation_matrix(float angle) {
+    float s = sinf(angle);
+    float c = cosf(angle);
+    return mat4(
+            c, s, 0.f, 0.f,
+            -s, c, 0.f, 0.f,
+            0.f, 0.f, 1.f, 0.f,
+            0.f, 0.f, 0.f, 1.f
+    );
+}
+
+// automatic mat4 printing
+std::ostream &operator<<(std::ostream &out, const mat4 &v);
+
+/**
+ *  Dnum class
+ * @tparam T
+ * @param f0
+ * @param d0
+ */
+//---------------------------
+template<class T>
+struct Dnum { // Dual numbers for automatic derivation
+//---------------------------
+    float f; //function value
+    T d; //derivatives
+    Dnum(float f0 = 0, T d0 = T(0)) { f = f0, d = d0; }
+
+    Dnum operator+(Dnum r) { return Dnum(f + r.f, d + r.d); }
+
+    Dnum operator-(Dnum r) { return Dnum(f - r.f, d - r.d); }
+
+    Dnum operator*(Dnum r) {
+        return Dnum(f * r.f, f * r.d + d * r.f);
+    }
+
+    Dnum operator/(Dnum r) {
+        return Dnum(f / r.f, (r.f * d - r.d * f) / r.f / r.f);
+    }
+};
+
+template<class T>
+Dnum<T> C(float t) { return Dnum<T>(t, 1); }
+
+template<class T>
+Dnum<T> Exp(Dnum<T> g) { return Dnum<T>(expf(g.f), expf(g.f) * g.d); }
+
+template<class T>
+Dnum<T> Sin(Dnum<T> g) { return Dnum<T>(sinf(g.f), cosf(g.f) * g.d); }
+
+template<class T>
+Dnum<T> Cos(Dnum<T> g) { return Dnum<T>(cosf(g.f), -sinf(g.f) * g.d); }
+
+template<class T>
+Dnum<T> Tan(Dnum<T> g) { return Sin(g) / Cos(g); }
+
+template<class T>
+Dnum<T> Sinh(Dnum<T> g) { return Dnum<T>(sinh(g.f), cosh(g.f) * g.d); }
+
+template<class T>
+Dnum<T> Cosh(Dnum<T> g) { return Dnum<T>(cosh(g.f), sinh(g.f) * g.d); }
+
+template<class T>
+Dnum<T> Tanh(Dnum<T> g) { return Sinh(g) / Cosh(g); }
+
+template<class T>
+Dnum<T> Log(Dnum<T> g) { return Dnum<T>(logf(g.f), g.d / g.f); }
+
+template<class T>
+Dnum<T> Pow(Dnum<T> g, float n) {
+    return Dnum<T>(powf(g.f, n), n * powf(g.f, n - 1) * g.d);
+}
+
+typedef Dnum<vec2> Dnum2;
+
+
+#endif //BRAVE2_MATH_H
diff --git a/src/utils/other/stb_image_write.h b/src/utils/other/stb_image_write.h
new file mode 100644
index 0000000..7cba2b7
--- /dev/null
+++ b/src/utils/other/stb_image_write.h
@@ -0,0 +1,1692 @@
+/* stb_image_write - v1.15 - public domain - http://nothings.org/stb
+   writes out PNG/BMP/TGA/JPEG/HDR images to C stdio - Sean Barrett 2010-2015
+                                     no warranty implied; use at your own risk
+
+   Before #including,
+
+       #define STB_IMAGE_WRITE_IMPLEMENTATION
+
+   in the file that you want to have the implementation.
+
+   Will probably not work correctly with strict-aliasing optimizations.
+
+ABOUT:
+
+   This header file is a library for writing images to C stdio or a callback.
+
+   The PNG output is not optimal; it is 20-50% larger than the file
+   written by a decent optimizing implementation; though providing a custom
+   zlib compress function (see STBIW_ZLIB_COMPRESS) can mitigate that.
+   This library is designed for source code compactness and simplicity,
+   not optimal image file size or run-time performance.
+
+BUILDING:
+
+   You can #define STBIW_ASSERT(x) before the #include to avoid using assert.h.
+   You can #define STBIW_MALLOC(), STBIW_REALLOC(), and STBIW_FREE() to replace
+   malloc,realloc,free.
+   You can #define STBIW_MEMMOVE() to replace memmove()
+   You can #define STBIW_ZLIB_COMPRESS to use a custom zlib-style compress function
+   for PNG compression (instead of the builtin one), it must have the following signature:
+   unsigned char * my_compress(unsigned char *data, int data_len, int *out_len, int quality);
+   The returned data will be freed with STBIW_FREE() (free() by default),
+   so it must be heap allocated with STBIW_MALLOC() (malloc() by default),
+
+UNICODE:
+
+   If compiling for Windows and you wish to use Unicode filenames, compile
+   with
+       #define STBIW_WINDOWS_UTF8
+   and pass utf8-encoded filenames. Call stbiw_convert_wchar_to_utf8 to convert
+   Windows wchar_t filenames to utf8.
+
+USAGE:
+
+   There are five functions, one for each image file format:
+
+     int stbi_write_png(char const *filename, int w, int h, int comp, const void *data, int stride_in_bytes);
+     int stbi_write_bmp(char const *filename, int w, int h, int comp, const void *data);
+     int stbi_write_tga(char const *filename, int w, int h, int comp, const void *data);
+     int stbi_write_jpg(char const *filename, int w, int h, int comp, const void *data, int quality);
+     int stbi_write_hdr(char const *filename, int w, int h, int comp, const float *data);
+
+     void stbi_flip_vertically_on_write(int flag); // flag is non-zero to flip data vertically
+
+   There are also five equivalent functions that use an arbitrary write function. You are
+   expected to open/close your file-equivalent before and after calling these:
+
+     int stbi_write_png_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const void  *data, int stride_in_bytes);
+     int stbi_write_bmp_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const void  *data);
+     int stbi_write_tga_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const void  *data);
+     int stbi_write_hdr_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const float *data);
+     int stbi_write_jpg_to_func(stbi_write_func *func, void *context, int x, int y, int comp, const void *data, int quality);
+
+   where the callback is:
+      void stbi_write_func(void *context, void *data, int size);
+
+   You can configure it with these global variables:
+      int stbi_write_tga_with_rle;             // defaults to true; set to 0 to disable RLE
+      int stbi_write_png_compression_level;    // defaults to 8; set to higher for more compression
+      int stbi_write_force_png_filter;         // defaults to -1; set to 0..5 to force a filter mode
+
+
+   You can define STBI_WRITE_NO_STDIO to disable the file variant of these
+   functions, so the library will not use stdio.h at all. However, this will
+   also disable HDR writing, because it requires stdio for formatted output.
+
+   Each function returns 0 on failure and non-0 on success.
+
+   The functions create an image file defined by the parameters. The image
+   is a rectangle of pixels stored from left-to-right, top-to-bottom.
+   Each pixel contains 'comp' channels of data stored interleaved with 8-bits
+   per channel, in the following order: 1=Y, 2=YA, 3=RGB, 4=RGBA. (Y is
+   monochrome color.) The rectangle is 'w' pixels wide and 'h' pixels tall.
+   The *data pointer points to the first byte of the top-left-most pixel.
+   For PNG, "stride_in_bytes" is the distance in bytes from the first byte of
+   a row of pixels to the first byte of the next row of pixels.
+
+   PNG creates output files with the same number of components as the input.
+   The BMP format expands Y to RGB in the file format and does not
+   output alpha.
+
+   PNG supports writing rectangles of data even when the bytes storing rows of
+   data are not consecutive in memory (e.g. sub-rectangles of a larger image),
+   by supplying the stride between the beginning of adjacent rows. The other
+   formats do not. (Thus you cannot write a native-format BMP through the BMP
+   writer, both because it is in BGR order and because it may have padding
+   at the end of the line.)
+
+   PNG allows you to set the deflate compression level by setting the global
+   variable 'stbi_write_png_compression_level' (it defaults to 8).
+
+   HDR expects linear float data. Since the format is always 32-bit rgb(e)
+   data, alpha (if provided) is discarded, and for monochrome data it is
+   replicated across all three channels.
+
+   TGA supports RLE or non-RLE compressed data. To use non-RLE-compressed
+   data, set the global variable 'stbi_write_tga_with_rle' to 0.
+
+   JPEG does ignore alpha channels in input data; quality is between 1 and 100.
+   Higher quality looks better but results in a bigger image.
+   JPEG baseline (no JPEG progressive).
+
+CREDITS:
+
+
+   Sean Barrett           -    PNG/BMP/TGA
+   Baldur Karlsson        -    HDR
+   Jean-Sebastien Guay    -    TGA monochrome
+   Tim Kelsey             -    misc enhancements
+   Alan Hickman           -    TGA RLE
+   Emmanuel Julien        -    initial file IO callback implementation
+   Jon Olick              -    original jo_jpeg.cpp code
+   Daniel Gibson          -    integrate JPEG, allow external zlib
+   Aarni Koskela          -    allow choosing PNG filter
+
+   bugfixes:
+      github:Chribba
+      Guillaume Chereau
+      github:jry2
+      github:romigrou
+      Sergio Gonzalez
+      Jonas Karlsson
+      Filip Wasil
+      Thatcher Ulrich
+      github:poppolopoppo
+      Patrick Boettcher
+      github:xeekworx
+      Cap Petschulat
+      Simon Rodriguez
+      Ivan Tikhonov
+      github:ignotion
+      Adam Schackart
+
+LICENSE
+
+  See end of file for license information.
+
+*/
+
+#ifndef INCLUDE_STB_IMAGE_WRITE_H
+#define INCLUDE_STB_IMAGE_WRITE_H
+
+#include <stdlib.h>
+
+// if STB_IMAGE_WRITE_STATIC causes problems, try defining STBIWDEF to 'inline' or 'static inline'
+#ifndef STBIWDEF
+#ifdef STB_IMAGE_WRITE_STATIC
+#define STBIWDEF  static
+#else
+#ifdef __cplusplus
+#define STBIWDEF  extern "C"
+#else
+#define STBIWDEF  extern
+#endif
+#endif
+#endif
+
+#ifndef STB_IMAGE_WRITE_STATIC  // C++ forbids static forward declarations
+extern int stbi_write_tga_with_rle;
+extern int stbi_write_png_compression_level;
+extern int stbi_write_force_png_filter;
+#endif
+
+#ifndef STBI_WRITE_NO_STDIO
+STBIWDEF int stbi_write_png(char const *filename, int w, int h, int comp, const void *data, int stride_in_bytes);
+STBIWDEF int stbi_write_bmp(char const *filename, int w, int h, int comp, const void *data);
+STBIWDEF int stbi_write_tga(char const *filename, int w, int h, int comp, const void *data);
+STBIWDEF int stbi_write_hdr(char const *filename, int w, int h, int comp, const float *data);
+STBIWDEF int stbi_write_jpg(char const *filename, int x, int y, int comp, const void *data, int quality);
+
+#ifdef STBI_WINDOWS_UTF8
+STBIWDEF int stbiw_convert_wchar_to_utf8(char *buffer, size_t bufferlen, const wchar_t* input);
+#endif
+#endif
+
+typedef void stbi_write_func(void *context, void *data, int size);
+
+STBIWDEF int stbi_write_png_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const void *data,
+                                    int stride_in_bytes);
+STBIWDEF int stbi_write_bmp_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const void *data);
+STBIWDEF int stbi_write_tga_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const void *data);
+STBIWDEF int stbi_write_hdr_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const float *data);
+STBIWDEF int
+stbi_write_jpg_to_func(stbi_write_func *func, void *context, int x, int y, int comp, const void *data, int quality);
+
+STBIWDEF void stbi_flip_vertically_on_write(int flip_boolean);
+
+#endif//INCLUDE_STB_IMAGE_WRITE_H
+
+#ifdef STB_IMAGE_WRITE_IMPLEMENTATION
+
+#ifdef _WIN32
+#ifndef _CRT_SECURE_NO_WARNINGS
+#define _CRT_SECURE_NO_WARNINGS
+#endif
+#ifndef _CRT_NONSTDC_NO_DEPRECATE
+#define _CRT_NONSTDC_NO_DEPRECATE
+#endif
+#endif
+
+#ifndef STBI_WRITE_NO_STDIO
+#include <stdio.h>
+#endif // STBI_WRITE_NO_STDIO
+
+#include <stdarg.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+#if defined(STBIW_MALLOC) && defined(STBIW_FREE) && (defined(STBIW_REALLOC) || defined(STBIW_REALLOC_SIZED))
+// ok
+#elif !defined(STBIW_MALLOC) && !defined(STBIW_FREE) && !defined(STBIW_REALLOC) && !defined(STBIW_REALLOC_SIZED)
+// ok
+#else
+#error "Must define all or none of STBIW_MALLOC, STBIW_FREE, and STBIW_REALLOC (or STBIW_REALLOC_SIZED)."
+#endif
+
+#ifndef STBIW_MALLOC
+#define STBIW_MALLOC(sz)        malloc(sz)
+#define STBIW_REALLOC(p,newsz)  realloc(p,newsz)
+#define STBIW_FREE(p)           free(p)
+#endif
+
+#ifndef STBIW_REALLOC_SIZED
+#define STBIW_REALLOC_SIZED(p,oldsz,newsz) STBIW_REALLOC(p,newsz)
+#endif
+
+
+#ifndef STBIW_MEMMOVE
+#define STBIW_MEMMOVE(a,b,sz) memmove(a,b,sz)
+#endif
+
+
+#ifndef STBIW_ASSERT
+#include <assert.h>
+#define STBIW_ASSERT(x) assert(x)
+#endif
+
+#define STBIW_UCHAR(x) (unsigned char) ((x) & 0xff)
+
+#ifdef STB_IMAGE_WRITE_STATIC
+static int stbi_write_png_compression_level = 8;
+static int stbi_write_tga_with_rle = 1;
+static int stbi_write_force_png_filter = -1;
+#else
+int stbi_write_png_compression_level = 8;
+int stbi_write_tga_with_rle = 1;
+int stbi_write_force_png_filter = -1;
+#endif
+
+static int stbi__flip_vertically_on_write = 0;
+
+STBIWDEF void stbi_flip_vertically_on_write(int flag)
+{
+   stbi__flip_vertically_on_write = flag;
+}
+
+typedef struct
+{
+   stbi_write_func *func;
+   void *context;
+   unsigned char buffer[64];
+   int buf_used;
+} stbi__write_context;
+
+// initialize a callback-based context
+static void stbi__start_write_callbacks(stbi__write_context *s, stbi_write_func *c, void *context)
+{
+   s->func    = c;
+   s->context = context;
+}
+
+#ifndef STBI_WRITE_NO_STDIO
+
+static void stbi__stdio_write(void *context, void *data, int size)
+{
+   fwrite(data,1,size,(FILE*) context);
+}
+
+#if defined(_MSC_VER) && defined(STBI_WINDOWS_UTF8)
+#ifdef __cplusplus
+#define STBIW_EXTERN extern "C"
+#else
+#define STBIW_EXTERN extern
+#endif
+STBIW_EXTERN __declspec(dllimport) int __stdcall MultiByteToWideChar(unsigned int cp, unsigned long flags, const char *str, int cbmb, wchar_t *widestr, int cchwide);
+STBIW_EXTERN __declspec(dllimport) int __stdcall WideCharToMultiByte(unsigned int cp, unsigned long flags, const wchar_t *widestr, int cchwide, char *str, int cbmb, const char *defchar, int *used_default);
+
+STBIWDEF int stbiw_convert_wchar_to_utf8(char *buffer, size_t bufferlen, const wchar_t* input)
+{
+    return WideCharToMultiByte(65001 /* UTF8 */, 0, input, -1, buffer, (int) bufferlen, NULL, NULL);
+}
+#endif
+
+static FILE *stbiw__fopen(char const *filename, char const *mode)
+{
+   FILE *f;
+#if defined(_MSC_VER) && defined(STBI_WINDOWS_UTF8)
+   wchar_t wMode[64];
+   wchar_t wFilename[1024];
+    if (0 == MultiByteToWideChar(65001 /* UTF8 */, 0, filename, -1, wFilename, sizeof(wFilename)))
+      return 0;
+
+    if (0 == MultiByteToWideChar(65001 /* UTF8 */, 0, mode, -1, wMode, sizeof(wMode)))
+      return 0;
+
+#if _MSC_VER >= 1400
+    if (0 != _wfopen_s(&f, wFilename, wMode))
+        f = 0;
+#else
+   f = _wfopen(wFilename, wMode);
+#endif
+
+#elif defined(_MSC_VER) && _MSC_VER >= 1400
+   if (0 != fopen_s(&f, filename, mode))
+      f=0;
+#else
+   f = fopen(filename, mode);
+#endif
+   return f;
+}
+
+static int stbi__start_write_file(stbi__write_context *s, const char *filename)
+{
+   FILE *f = stbiw__fopen(filename, "wb");
+   stbi__start_write_callbacks(s, stbi__stdio_write, (void *) f);
+   return f != NULL;
+}
+
+static void stbi__end_write_file(stbi__write_context *s)
+{
+   fclose((FILE *)s->context);
+}
+
+#endif // !STBI_WRITE_NO_STDIO
+
+typedef unsigned int stbiw_uint32;
+typedef int stb_image_write_test[sizeof(stbiw_uint32)==4 ? 1 : -1];
+
+static void stbiw__writefv(stbi__write_context *s, const char *fmt, va_list v)
+{
+   while (*fmt) {
+      switch (*fmt++) {
+         case ' ': break;
+         case '1': { unsigned char x = STBIW_UCHAR(va_arg(v, int));
+                     s->func(s->context,&x,1);
+                     break; }
+         case '2': { int x = va_arg(v,int);
+                     unsigned char b[2];
+                     b[0] = STBIW_UCHAR(x);
+                     b[1] = STBIW_UCHAR(x>>8);
+                     s->func(s->context,b,2);
+                     break; }
+         case '4': { stbiw_uint32 x = va_arg(v,int);
+                     unsigned char b[4];
+                     b[0]=STBIW_UCHAR(x);
+                     b[1]=STBIW_UCHAR(x>>8);
+                     b[2]=STBIW_UCHAR(x>>16);
+                     b[3]=STBIW_UCHAR(x>>24);
+                     s->func(s->context,b,4);
+                     break; }
+         default:
+            STBIW_ASSERT(0);
+            return;
+      }
+   }
+}
+
+static void stbiw__writef(stbi__write_context *s, const char *fmt, ...)
+{
+   va_list v;
+   va_start(v, fmt);
+   stbiw__writefv(s, fmt, v);
+   va_end(v);
+}
+
+static void stbiw__write_flush(stbi__write_context *s)
+{
+   if (s->buf_used) {
+      s->func(s->context, &s->buffer, s->buf_used);
+      s->buf_used = 0;
+   }
+}
+
+static void stbiw__putc(stbi__write_context *s, unsigned char c)
+{
+   s->func(s->context, &c, 1);
+}
+
+static void stbiw__write1(stbi__write_context *s, unsigned char a)
+{
+   if (s->buf_used + 1 > sizeof(s->buffer))
+      stbiw__write_flush(s);
+   s->buffer[s->buf_used++] = a;
+}
+
+static void stbiw__write3(stbi__write_context *s, unsigned char a, unsigned char b, unsigned char c)
+{
+   int n;
+   if (s->buf_used + 3 > sizeof(s->buffer))
+      stbiw__write_flush(s);
+   n = s->buf_used;
+   s->buf_used = n+3;
+   s->buffer[n+0] = a;
+   s->buffer[n+1] = b;
+   s->buffer[n+2] = c;
+}
+
+static void stbiw__write_pixel(stbi__write_context *s, int rgb_dir, int comp, int write_alpha, int expand_mono, unsigned char *d)
+{
+   unsigned char bg[3] = { 255, 0, 255}, px[3];
+   int k;
+
+   if (write_alpha < 0)
+      stbiw__write1(s, d[comp - 1]);
+
+   switch (comp) {
+      case 2: // 2 pixels = mono + alpha, alpha is written separately, so same as 1-channel case
+      case 1:
+         if (expand_mono)
+            stbiw__write3(s, d[0], d[0], d[0]); // monochrome bmp
+         else
+            stbiw__write1(s, d[0]);  // monochrome TGA
+         break;
+      case 4:
+         if (!write_alpha) {
+            // composite against pink background
+            for (k = 0; k < 3; ++k)
+               px[k] = bg[k] + ((d[k] - bg[k]) * d[3]) / 255;
+            stbiw__write3(s, px[1 - rgb_dir], px[1], px[1 + rgb_dir]);
+            break;
+         }
+         /* FALLTHROUGH */
+      case 3:
+         stbiw__write3(s, d[1 - rgb_dir], d[1], d[1 + rgb_dir]);
+         break;
+   }
+   if (write_alpha > 0)
+      stbiw__write1(s, d[comp - 1]);
+}
+
+static void stbiw__write_pixels(stbi__write_context *s, int rgb_dir, int vdir, int x, int y, int comp, void *data, int write_alpha, int scanline_pad, int expand_mono)
+{
+   stbiw_uint32 zero = 0;
+   int i,j, j_end;
+
+   if (y <= 0)
+      return;
+
+   if (stbi__flip_vertically_on_write)
+      vdir *= -1;
+
+   if (vdir < 0) {
+      j_end = -1; j = y-1;
+   } else {
+      j_end =  y; j = 0;
+   }
+
+   for (; j != j_end; j += vdir) {
+      for (i=0; i < x; ++i) {
+         unsigned char *d = (unsigned char *) data + (j*x+i)*comp;
+         stbiw__write_pixel(s, rgb_dir, comp, write_alpha, expand_mono, d);
+      }
+      stbiw__write_flush(s);
+      s->func(s->context, &zero, scanline_pad);
+   }
+}
+
+static int stbiw__outfile(stbi__write_context *s, int rgb_dir, int vdir, int x, int y, int comp, int expand_mono, void *data, int alpha, int pad, const char *fmt, ...)
+{
+   if (y < 0 || x < 0) {
+      return 0;
+   } else {
+      va_list v;
+      va_start(v, fmt);
+      stbiw__writefv(s, fmt, v);
+      va_end(v);
+      stbiw__write_pixels(s,rgb_dir,vdir,x,y,comp,data,alpha,pad, expand_mono);
+      return 1;
+   }
+}
+
+static int stbi_write_bmp_core(stbi__write_context *s, int x, int y, int comp, const void *data)
+{
+   int pad = (-x*3) & 3;
+   return stbiw__outfile(s,-1,-1,x,y,comp,1,(void *) data,0,pad,
+           "11 4 22 4" "4 44 22 444444",
+           'B', 'M', 14+40+(x*3+pad)*y, 0,0, 14+40,  // file header
+            40, x,y, 1,24, 0,0,0,0,0,0);             // bitmap header
+}
+
+STBIWDEF int stbi_write_bmp_to_func(stbi_write_func *func, void *context, int x, int y, int comp, const void *data)
+{
+   stbi__write_context s = { 0 };
+   stbi__start_write_callbacks(&s, func, context);
+   return stbi_write_bmp_core(&s, x, y, comp, data);
+}
+
+#ifndef STBI_WRITE_NO_STDIO
+STBIWDEF int stbi_write_bmp(char const *filename, int x, int y, int comp, const void *data)
+{
+   stbi__write_context s = { 0 };
+   if (stbi__start_write_file(&s,filename)) {
+      int r = stbi_write_bmp_core(&s, x, y, comp, data);
+      stbi__end_write_file(&s);
+      return r;
+   } else
+      return 0;
+}
+#endif //!STBI_WRITE_NO_STDIO
+
+static int stbi_write_tga_core(stbi__write_context *s, int x, int y, int comp, void *data)
+{
+   int has_alpha = (comp == 2 || comp == 4);
+   int colorbytes = has_alpha ? comp-1 : comp;
+   int format = colorbytes < 2 ? 3 : 2; // 3 color channels (RGB/RGBA) = 2, 1 color channel (Y/YA) = 3
+
+   if (y < 0 || x < 0)
+      return 0;
+
+   if (!stbi_write_tga_with_rle) {
+      return stbiw__outfile(s, -1, -1, x, y, comp, 0, (void *) data, has_alpha, 0,
+         "111 221 2222 11", 0, 0, format, 0, 0, 0, 0, 0, x, y, (colorbytes + has_alpha) * 8, has_alpha * 8);
+   } else {
+      int i,j,k;
+      int jend, jdir;
+
+      stbiw__writef(s, "111 221 2222 11", 0,0,format+8, 0,0,0, 0,0,x,y, (colorbytes + has_alpha) * 8, has_alpha * 8);
+
+      if (stbi__flip_vertically_on_write) {
+         j = 0;
+         jend = y;
+         jdir = 1;
+      } else {
+         j = y-1;
+         jend = -1;
+         jdir = -1;
+      }
+      for (; j != jend; j += jdir) {
+         unsigned char *row = (unsigned char *) data + j * x * comp;
+         int len;
+
+         for (i = 0; i < x; i += len) {
+            unsigned char *begin = row + i * comp;
+            int diff = 1;
+            len = 1;
+
+            if (i < x - 1) {
+               ++len;
+               diff = memcmp(begin, row + (i + 1) * comp, comp);
+               if (diff) {
+                  const unsigned char *prev = begin;
+                  for (k = i + 2; k < x && len < 128; ++k) {
+                     if (memcmp(prev, row + k * comp, comp)) {
+                        prev += comp;
+                        ++len;
+                     } else {
+                        --len;
+                        break;
+                     }
+                  }
+               } else {
+                  for (k = i + 2; k < x && len < 128; ++k) {
+                     if (!memcmp(begin, row + k * comp, comp)) {
+                        ++len;
+                     } else {
+                        break;
+                     }
+                  }
+               }
+            }
+
+            if (diff) {
+               unsigned char header = STBIW_UCHAR(len - 1);
+               stbiw__write1(s, header);
+               for (k = 0; k < len; ++k) {
+                  stbiw__write_pixel(s, -1, comp, has_alpha, 0, begin + k * comp);
+               }
+            } else {
+               unsigned char header = STBIW_UCHAR(len - 129);
+               stbiw__write1(s, header);
+               stbiw__write_pixel(s, -1, comp, has_alpha, 0, begin);
+            }
+         }
+      }
+      stbiw__write_flush(s);
+   }
+   return 1;
+}
+
+STBIWDEF int stbi_write_tga_to_func(stbi_write_func *func, void *context, int x, int y, int comp, const void *data)
+{
+   stbi__write_context s = { 0 };
+   stbi__start_write_callbacks(&s, func, context);
+   return stbi_write_tga_core(&s, x, y, comp, (void *) data);
+}
+
+#ifndef STBI_WRITE_NO_STDIO
+STBIWDEF int stbi_write_tga(char const *filename, int x, int y, int comp, const void *data)
+{
+   stbi__write_context s = { 0 };
+   if (stbi__start_write_file(&s,filename)) {
+      int r = stbi_write_tga_core(&s, x, y, comp, (void *) data);
+      stbi__end_write_file(&s);
+      return r;
+   } else
+      return 0;
+}
+#endif
+
+// *************************************************************************************************
+// Radiance RGBE HDR writer
+// by Baldur Karlsson
+
+#define stbiw__max(a, b)  ((a) > (b) ? (a) : (b))
+
+static void stbiw__linear_to_rgbe(unsigned char *rgbe, float *linear)
+{
+   int exponent;
+   float maxcomp = stbiw__max(linear[0], stbiw__max(linear[1], linear[2]));
+
+   if (maxcomp < 1e-32f) {
+      rgbe[0] = rgbe[1] = rgbe[2] = rgbe[3] = 0;
+   } else {
+      float normalize = (float) frexp(maxcomp, &exponent) * 256.0f/maxcomp;
+
+      rgbe[0] = (unsigned char)(linear[0] * normalize);
+      rgbe[1] = (unsigned char)(linear[1] * normalize);
+      rgbe[2] = (unsigned char)(linear[2] * normalize);
+      rgbe[3] = (unsigned char)(exponent + 128);
+   }
+}
+
+static void stbiw__write_run_data(stbi__write_context *s, int length, unsigned char databyte)
+{
+   unsigned char lengthbyte = STBIW_UCHAR(length+128);
+   STBIW_ASSERT(length+128 <= 255);
+   s->func(s->context, &lengthbyte, 1);
+   s->func(s->context, &databyte, 1);
+}
+
+static void stbiw__write_dump_data(stbi__write_context *s, int length, unsigned char *data)
+{
+   unsigned char lengthbyte = STBIW_UCHAR(length);
+   STBIW_ASSERT(length <= 128); // inconsistent with spec but consistent with official code
+   s->func(s->context, &lengthbyte, 1);
+   s->func(s->context, data, length);
+}
+
+static void stbiw__write_hdr_scanline(stbi__write_context *s, int width, int ncomp, unsigned char *scratch, float *scanline)
+{
+   unsigned char scanlineheader[4] = { 2, 2, 0, 0 };
+   unsigned char rgbe[4];
+   float linear[3];
+   int x;
+
+   scanlineheader[2] = (width&0xff00)>>8;
+   scanlineheader[3] = (width&0x00ff);
+
+   /* skip RLE for images too small or large */
+   if (width < 8 || width >= 32768) {
+      for (x=0; x < width; x++) {
+         switch (ncomp) {
+            case 4: /* fallthrough */
+            case 3: linear[2] = scanline[x*ncomp + 2];
+                    linear[1] = scanline[x*ncomp + 1];
+                    linear[0] = scanline[x*ncomp + 0];
+                    break;
+            default:
+                    linear[0] = linear[1] = linear[2] = scanline[x*ncomp + 0];
+                    break;
+         }
+         stbiw__linear_to_rgbe(rgbe, linear);
+         s->func(s->context, rgbe, 4);
+      }
+   } else {
+      int c,r;
+      /* encode into scratch buffer */
+      for (x=0; x < width; x++) {
+         switch(ncomp) {
+            case 4: /* fallthrough */
+            case 3: linear[2] = scanline[x*ncomp + 2];
+                    linear[1] = scanline[x*ncomp + 1];
+                    linear[0] = scanline[x*ncomp + 0];
+                    break;
+            default:
+                    linear[0] = linear[1] = linear[2] = scanline[x*ncomp + 0];
+                    break;
+         }
+         stbiw__linear_to_rgbe(rgbe, linear);
+         scratch[x + width*0] = rgbe[0];
+         scratch[x + width*1] = rgbe[1];
+         scratch[x + width*2] = rgbe[2];
+         scratch[x + width*3] = rgbe[3];
+      }
+
+      s->func(s->context, scanlineheader, 4);
+
+      /* RLE each component separately */
+      for (c=0; c < 4; c++) {
+         unsigned char *comp = &scratch[width*c];
+
+         x = 0;
+         while (x < width) {
+            // find first run
+            r = x;
+            while (r+2 < width) {
+               if (comp[r] == comp[r+1] && comp[r] == comp[r+2])
+                  break;
+               ++r;
+            }
+            if (r+2 >= width)
+               r = width;
+            // dump up to first run
+            while (x < r) {
+               int len = r-x;
+               if (len > 128) len = 128;
+               stbiw__write_dump_data(s, len, &comp[x]);
+               x += len;
+            }
+            // if there's a run, output it
+            if (r+2 < width) { // same test as what we break out of in search loop, so only true if we break'd
+               // find next byte after run
+               while (r < width && comp[r] == comp[x])
+                  ++r;
+               // output run up to r
+               while (x < r) {
+                  int len = r-x;
+                  if (len > 127) len = 127;
+                  stbiw__write_run_data(s, len, comp[x]);
+                  x += len;
+               }
+            }
+         }
+      }
+   }
+}
+
+static int stbi_write_hdr_core(stbi__write_context *s, int x, int y, int comp, float *data)
+{
+   if (y <= 0 || x <= 0 || data == NULL)
+      return 0;
+   else {
+      // Each component is stored separately. Allocate scratch space for full output scanline.
+      unsigned char *scratch = (unsigned char *) STBIW_MALLOC(x*4);
+      int i, len;
+      char buffer[128];
+      char header[] = "#?RADIANCE\n# Written by stb_image_write.h\nFORMAT=32-bit_rle_rgbe\n";
+      s->func(s->context, header, sizeof(header)-1);
+
+#ifdef __STDC_WANT_SECURE_LIB__
+      len = sprintf_s(buffer, sizeof(buffer), "EXPOSURE=          1.0000000000000\n\n-Y %d +X %d\n", y, x);
+#else
+      len = sprintf(buffer, "EXPOSURE=          1.0000000000000\n\n-Y %d +X %d\n", y, x);
+#endif
+      s->func(s->context, buffer, len);
+
+      for(i=0; i < y; i++)
+         stbiw__write_hdr_scanline(s, x, comp, scratch, data + comp*x*(stbi__flip_vertically_on_write ? y-1-i : i));
+      STBIW_FREE(scratch);
+      return 1;
+   }
+}
+
+STBIWDEF int stbi_write_hdr_to_func(stbi_write_func *func, void *context, int x, int y, int comp, const float *data)
+{
+   stbi__write_context s = { 0 };
+   stbi__start_write_callbacks(&s, func, context);
+   return stbi_write_hdr_core(&s, x, y, comp, (float *) data);
+}
+
+#ifndef STBI_WRITE_NO_STDIO
+STBIWDEF int stbi_write_hdr(char const *filename, int x, int y, int comp, const float *data)
+{
+   stbi__write_context s = { 0 };
+   if (stbi__start_write_file(&s,filename)) {
+      int r = stbi_write_hdr_core(&s, x, y, comp, (float *) data);
+      stbi__end_write_file(&s);
+      return r;
+   } else
+      return 0;
+}
+#endif // STBI_WRITE_NO_STDIO
+
+
+//////////////////////////////////////////////////////////////////////////////
+//
+// PNG writer
+//
+
+#ifndef STBIW_ZLIB_COMPRESS
+// stretchy buffer; stbiw__sbpush() == vector<>::push_back() -- stbiw__sbcount() == vector<>::size()
+#define stbiw__sbraw(a) ((int *) (void *) (a) - 2)
+#define stbiw__sbm(a)   stbiw__sbraw(a)[0]
+#define stbiw__sbn(a)   stbiw__sbraw(a)[1]
+
+#define stbiw__sbneedgrow(a,n)  ((a)==0 || stbiw__sbn(a)+n >= stbiw__sbm(a))
+#define stbiw__sbmaybegrow(a,n) (stbiw__sbneedgrow(a,(n)) ? stbiw__sbgrow(a,n) : 0)
+#define stbiw__sbgrow(a,n)  stbiw__sbgrowf((void **) &(a), (n), sizeof(*(a)))
+
+#define stbiw__sbpush(a, v)      (stbiw__sbmaybegrow(a,1), (a)[stbiw__sbn(a)++] = (v))
+#define stbiw__sbcount(a)        ((a) ? stbiw__sbn(a) : 0)
+#define stbiw__sbfree(a)         ((a) ? STBIW_FREE(stbiw__sbraw(a)),0 : 0)
+
+static void *stbiw__sbgrowf(void **arr, int increment, int itemsize)
+{
+   int m = *arr ? 2*stbiw__sbm(*arr)+increment : increment+1;
+   void *p = STBIW_REALLOC_SIZED(*arr ? stbiw__sbraw(*arr) : 0, *arr ? (stbiw__sbm(*arr)*itemsize + sizeof(int)*2) : 0, itemsize * m + sizeof(int)*2);
+   STBIW_ASSERT(p);
+   if (p) {
+      if (!*arr) ((int *) p)[1] = 0;
+      *arr = (void *) ((int *) p + 2);
+      stbiw__sbm(*arr) = m;
+   }
+   return *arr;
+}
+
+static unsigned char *stbiw__zlib_flushf(unsigned char *data, unsigned int *bitbuffer, int *bitcount)
+{
+   while (*bitcount >= 8) {
+      stbiw__sbpush(data, STBIW_UCHAR(*bitbuffer));
+      *bitbuffer >>= 8;
+      *bitcount -= 8;
+   }
+   return data;
+}
+
+static int stbiw__zlib_bitrev(int code, int codebits)
+{
+   int res=0;
+   while (codebits--) {
+      res = (res << 1) | (code & 1);
+      code >>= 1;
+   }
+   return res;
+}
+
+static unsigned int stbiw__zlib_countm(unsigned char *a, unsigned char *b, int limit)
+{
+   int i;
+   for (i=0; i < limit && i < 258; ++i)
+      if (a[i] != b[i]) break;
+   return i;
+}
+
+static unsigned int stbiw__zhash(unsigned char *data)
+{
+   stbiw_uint32 hash = data[0] + (data[1] << 8) + (data[2] << 16);
+   hash ^= hash << 3;
+   hash += hash >> 5;
+   hash ^= hash << 4;
+   hash += hash >> 17;
+   hash ^= hash << 25;
+   hash += hash >> 6;
+   return hash;
+}
+
+#define stbiw__zlib_flush() (out = stbiw__zlib_flushf(out, &bitbuf, &bitcount))
+#define stbiw__zlib_add(code,codebits) \
+      (bitbuf |= (code) << bitcount, bitcount += (codebits), stbiw__zlib_flush())
+#define stbiw__zlib_huffa(b,c)  stbiw__zlib_add(stbiw__zlib_bitrev(b,c),c)
+// default huffman tables
+#define stbiw__zlib_huff1(n)  stbiw__zlib_huffa(0x30 + (n), 8)
+#define stbiw__zlib_huff2(n)  stbiw__zlib_huffa(0x190 + (n)-144, 9)
+#define stbiw__zlib_huff3(n)  stbiw__zlib_huffa(0 + (n)-256,7)
+#define stbiw__zlib_huff4(n)  stbiw__zlib_huffa(0xc0 + (n)-280,8)
+#define stbiw__zlib_huff(n)  ((n) <= 143 ? stbiw__zlib_huff1(n) : (n) <= 255 ? stbiw__zlib_huff2(n) : (n) <= 279 ? stbiw__zlib_huff3(n) : stbiw__zlib_huff4(n))
+#define stbiw__zlib_huffb(n) ((n) <= 143 ? stbiw__zlib_huff1(n) : stbiw__zlib_huff2(n))
+
+#define stbiw__ZHASH   16384
+
+#endif // STBIW_ZLIB_COMPRESS
+
+STBIWDEF unsigned char * stbi_zlib_compress(unsigned char *data, int data_len, int *out_len, int quality)
+{
+#ifdef STBIW_ZLIB_COMPRESS
+   // user provided a zlib compress implementation, use that
+   return STBIW_ZLIB_COMPRESS(data, data_len, out_len, quality);
+#else // use builtin
+   static unsigned short lengthc[] = { 3,4,5,6,7,8,9,10,11,13,15,17,19,23,27,31,35,43,51,59,67,83,99,115,131,163,195,227,258, 259 };
+   static unsigned char  lengtheb[]= { 0,0,0,0,0,0,0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4,  4,  5,  5,  5,  5,  0 };
+   static unsigned short distc[]   = { 1,2,3,4,5,7,9,13,17,25,33,49,65,97,129,193,257,385,513,769,1025,1537,2049,3073,4097,6145,8193,12289,16385,24577, 32768 };
+   static unsigned char  disteb[]  = { 0,0,0,0,1,1,2,2,3,3,4,4,5,5,6,6,7,7,8,8,9,9,10,10,11,11,12,12,13,13 };
+   unsigned int bitbuf=0;
+   int i,j, bitcount=0;
+   unsigned char *out = NULL;
+   unsigned char ***hash_table = (unsigned char***) STBIW_MALLOC(stbiw__ZHASH * sizeof(unsigned char**));
+   if (hash_table == NULL)
+      return NULL;
+   if (quality < 5) quality = 5;
+
+   stbiw__sbpush(out, 0x78);   // DEFLATE 32K window
+   stbiw__sbpush(out, 0x5e);   // FLEVEL = 1
+   stbiw__zlib_add(1,1);  // BFINAL = 1
+   stbiw__zlib_add(1,2);  // BTYPE = 1 -- fixed huffman
+
+   for (i=0; i < stbiw__ZHASH; ++i)
+      hash_table[i] = NULL;
+
+   i=0;
+   while (i < data_len-3) {
+      // hash next 3 bytes of data to be compressed
+      int h = stbiw__zhash(data+i)&(stbiw__ZHASH-1), best=3;
+      unsigned char *bestloc = 0;
+      unsigned char **hlist = hash_table[h];
+      int n = stbiw__sbcount(hlist);
+      for (j=0; j < n; ++j) {
+         if (hlist[j]-data > i-32768) { // if entry lies within window
+            int d = stbiw__zlib_countm(hlist[j], data+i, data_len-i);
+            if (d >= best) { best=d; bestloc=hlist[j]; }
+         }
+      }
+      // when hash table entry is too long, delete half the entries
+      if (hash_table[h] && stbiw__sbn(hash_table[h]) == 2*quality) {
+         STBIW_MEMMOVE(hash_table[h], hash_table[h]+quality, sizeof(hash_table[h][0])*quality);
+         stbiw__sbn(hash_table[h]) = quality;
+      }
+      stbiw__sbpush(hash_table[h],data+i);
+
+      if (bestloc) {
+         // "lazy matching" - check match at *next* byte, and if it's better, do cur byte as literal
+         h = stbiw__zhash(data+i+1)&(stbiw__ZHASH-1);
+         hlist = hash_table[h];
+         n = stbiw__sbcount(hlist);
+         for (j=0; j < n; ++j) {
+            if (hlist[j]-data > i-32767) {
+               int e = stbiw__zlib_countm(hlist[j], data+i+1, data_len-i-1);
+               if (e > best) { // if next match is better, bail on current match
+                  bestloc = NULL;
+                  break;
+               }
+            }
+         }
+      }
+
+      if (bestloc) {
+         int d = (int) (data+i - bestloc); // distance back
+         STBIW_ASSERT(d <= 32767 && best <= 258);
+         for (j=0; best > lengthc[j+1]-1; ++j);
+         stbiw__zlib_huff(j+257);
+         if (lengtheb[j]) stbiw__zlib_add(best - lengthc[j], lengtheb[j]);
+         for (j=0; d > distc[j+1]-1; ++j);
+         stbiw__zlib_add(stbiw__zlib_bitrev(j,5),5);
+         if (disteb[j]) stbiw__zlib_add(d - distc[j], disteb[j]);
+         i += best;
+      } else {
+         stbiw__zlib_huffb(data[i]);
+         ++i;
+      }
+   }
+   // write out final bytes
+   for (;i < data_len; ++i)
+      stbiw__zlib_huffb(data[i]);
+   stbiw__zlib_huff(256); // end of block
+   // pad with 0 bits to byte boundary
+   while (bitcount)
+      stbiw__zlib_add(0,1);
+
+   for (i=0; i < stbiw__ZHASH; ++i)
+      (void) stbiw__sbfree(hash_table[i]);
+   STBIW_FREE(hash_table);
+
+   {
+      // compute adler32 on input
+      unsigned int s1=1, s2=0;
+      int blocklen = (int) (data_len % 5552);
+      j=0;
+      while (j < data_len) {
+         for (i=0; i < blocklen; ++i) { s1 += data[j+i]; s2 += s1; }
+         s1 %= 65521; s2 %= 65521;
+         j += blocklen;
+         blocklen = 5552;
+      }
+      stbiw__sbpush(out, STBIW_UCHAR(s2 >> 8));
+      stbiw__sbpush(out, STBIW_UCHAR(s2));
+      stbiw__sbpush(out, STBIW_UCHAR(s1 >> 8));
+      stbiw__sbpush(out, STBIW_UCHAR(s1));
+   }
+   *out_len = stbiw__sbn(out);
+   // make returned pointer freeable
+   STBIW_MEMMOVE(stbiw__sbraw(out), out, *out_len);
+   return (unsigned char *) stbiw__sbraw(out);
+#endif // STBIW_ZLIB_COMPRESS
+}
+
+static unsigned int stbiw__crc32(unsigned char *buffer, int len)
+{
+#ifdef STBIW_CRC32
+    return STBIW_CRC32(buffer, len);
+#else
+   static unsigned int crc_table[256] =
+   {
+      0x00000000, 0x77073096, 0xEE0E612C, 0x990951BA, 0x076DC419, 0x706AF48F, 0xE963A535, 0x9E6495A3,
+      0x0eDB8832, 0x79DCB8A4, 0xE0D5E91E, 0x97D2D988, 0x09B64C2B, 0x7EB17CBD, 0xE7B82D07, 0x90BF1D91,
+      0x1DB71064, 0x6AB020F2, 0xF3B97148, 0x84BE41DE, 0x1ADAD47D, 0x6DDDE4EB, 0xF4D4B551, 0x83D385C7,
+      0x136C9856, 0x646BA8C0, 0xFD62F97A, 0x8A65C9EC, 0x14015C4F, 0x63066CD9, 0xFA0F3D63, 0x8D080DF5,
+      0x3B6E20C8, 0x4C69105E, 0xD56041E4, 0xA2677172, 0x3C03E4D1, 0x4B04D447, 0xD20D85FD, 0xA50AB56B,
+      0x35B5A8FA, 0x42B2986C, 0xDBBBC9D6, 0xACBCF940, 0x32D86CE3, 0x45DF5C75, 0xDCD60DCF, 0xABD13D59,
+      0x26D930AC, 0x51DE003A, 0xC8D75180, 0xBFD06116, 0x21B4F4B5, 0x56B3C423, 0xCFBA9599, 0xB8BDA50F,
+      0x2802B89E, 0x5F058808, 0xC60CD9B2, 0xB10BE924, 0x2F6F7C87, 0x58684C11, 0xC1611DAB, 0xB6662D3D,
+      0x76DC4190, 0x01DB7106, 0x98D220BC, 0xEFD5102A, 0x71B18589, 0x06B6B51F, 0x9FBFE4A5, 0xE8B8D433,
+      0x7807C9A2, 0x0F00F934, 0x9609A88E, 0xE10E9818, 0x7F6A0DBB, 0x086D3D2D, 0x91646C97, 0xE6635C01,
+      0x6B6B51F4, 0x1C6C6162, 0x856530D8, 0xF262004E, 0x6C0695ED, 0x1B01A57B, 0x8208F4C1, 0xF50FC457,
+      0x65B0D9C6, 0x12B7E950, 0x8BBEB8EA, 0xFCB9887C, 0x62DD1DDF, 0x15DA2D49, 0x8CD37CF3, 0xFBD44C65,
+      0x4DB26158, 0x3AB551CE, 0xA3BC0074, 0xD4BB30E2, 0x4ADFA541, 0x3DD895D7, 0xA4D1C46D, 0xD3D6F4FB,
+      0x4369E96A, 0x346ED9FC, 0xAD678846, 0xDA60B8D0, 0x44042D73, 0x33031DE5, 0xAA0A4C5F, 0xDD0D7CC9,
+      0x5005713C, 0x270241AA, 0xBE0B1010, 0xC90C2086, 0x5768B525, 0x206F85B3, 0xB966D409, 0xCE61E49F,
+      0x5EDEF90E, 0x29D9C998, 0xB0D09822, 0xC7D7A8B4, 0x59B33D17, 0x2EB40D81, 0xB7BD5C3B, 0xC0BA6CAD,
+      0xEDB88320, 0x9ABFB3B6, 0x03B6E20C, 0x74B1D29A, 0xEAD54739, 0x9DD277AF, 0x04DB2615, 0x73DC1683,
+      0xE3630B12, 0x94643B84, 0x0D6D6A3E, 0x7A6A5AA8, 0xE40ECF0B, 0x9309FF9D, 0x0A00AE27, 0x7D079EB1,
+      0xF00F9344, 0x8708A3D2, 0x1E01F268, 0x6906C2FE, 0xF762575D, 0x806567CB, 0x196C3671, 0x6E6B06E7,
+      0xFED41B76, 0x89D32BE0, 0x10DA7A5A, 0x67DD4ACC, 0xF9B9DF6F, 0x8EBEEFF9, 0x17B7BE43, 0x60B08ED5,
+      0xD6D6A3E8, 0xA1D1937E, 0x38D8C2C4, 0x4FDFF252, 0xD1BB67F1, 0xA6BC5767, 0x3FB506DD, 0x48B2364B,
+      0xD80D2BDA, 0xAF0A1B4C, 0x36034AF6, 0x41047A60, 0xDF60EFC3, 0xA867DF55, 0x316E8EEF, 0x4669BE79,
+      0xCB61B38C, 0xBC66831A, 0x256FD2A0, 0x5268E236, 0xCC0C7795, 0xBB0B4703, 0x220216B9, 0x5505262F,
+      0xC5BA3BBE, 0xB2BD0B28, 0x2BB45A92, 0x5CB36A04, 0xC2D7FFA7, 0xB5D0CF31, 0x2CD99E8B, 0x5BDEAE1D,
+      0x9B64C2B0, 0xEC63F226, 0x756AA39C, 0x026D930A, 0x9C0906A9, 0xEB0E363F, 0x72076785, 0x05005713,
+      0x95BF4A82, 0xE2B87A14, 0x7BB12BAE, 0x0CB61B38, 0x92D28E9B, 0xE5D5BE0D, 0x7CDCEFB7, 0x0BDBDF21,
+      0x86D3D2D4, 0xF1D4E242, 0x68DDB3F8, 0x1FDA836E, 0x81BE16CD, 0xF6B9265B, 0x6FB077E1, 0x18B74777,
+      0x88085AE6, 0xFF0F6A70, 0x66063BCA, 0x11010B5C, 0x8F659EFF, 0xF862AE69, 0x616BFFD3, 0x166CCF45,
+      0xA00AE278, 0xD70DD2EE, 0x4E048354, 0x3903B3C2, 0xA7672661, 0xD06016F7, 0x4969474D, 0x3E6E77DB,
+      0xAED16A4A, 0xD9D65ADC, 0x40DF0B66, 0x37D83BF0, 0xA9BCAE53, 0xDEBB9EC5, 0x47B2CF7F, 0x30B5FFE9,
+      0xBDBDF21C, 0xCABAC28A, 0x53B39330, 0x24B4A3A6, 0xBAD03605, 0xCDD70693, 0x54DE5729, 0x23D967BF,
+      0xB3667A2E, 0xC4614AB8, 0x5D681B02, 0x2A6F2B94, 0xB40BBE37, 0xC30C8EA1, 0x5A05DF1B, 0x2D02EF8D
+   };
+
+   unsigned int crc = ~0u;
+   int i;
+   for (i=0; i < len; ++i)
+      crc = (crc >> 8) ^ crc_table[buffer[i] ^ (crc & 0xff)];
+   return ~crc;
+#endif
+}
+
+#define stbiw__wpng4(o,a,b,c,d) ((o)[0]=STBIW_UCHAR(a),(o)[1]=STBIW_UCHAR(b),(o)[2]=STBIW_UCHAR(c),(o)[3]=STBIW_UCHAR(d),(o)+=4)
+#define stbiw__wp32(data,v) stbiw__wpng4(data, (v)>>24,(v)>>16,(v)>>8,(v));
+#define stbiw__wptag(data,s) stbiw__wpng4(data, s[0],s[1],s[2],s[3])
+
+static void stbiw__wpcrc(unsigned char **data, int len)
+{
+   unsigned int crc = stbiw__crc32(*data - len - 4, len+4);
+   stbiw__wp32(*data, crc);
+}
+
+static unsigned char stbiw__paeth(int a, int b, int c)
+{
+   int p = a + b - c, pa = abs(p-a), pb = abs(p-b), pc = abs(p-c);
+   if (pa <= pb && pa <= pc) return STBIW_UCHAR(a);
+   if (pb <= pc) return STBIW_UCHAR(b);
+   return STBIW_UCHAR(c);
+}
+
+// @OPTIMIZE: provide an option that always forces left-predict or paeth predict
+static void stbiw__encode_png_line(unsigned char *pixels, int stride_bytes, int width, int height, int y, int n, int filter_type, signed char *line_buffer)
+{
+   static int mapping[] = { 0,1,2,3,4 };
+   static int firstmap[] = { 0,1,0,5,6 };
+   int *mymap = (y != 0) ? mapping : firstmap;
+   int i;
+   int type = mymap[filter_type];
+   unsigned char *z = pixels + stride_bytes * (stbi__flip_vertically_on_write ? height-1-y : y);
+   int signed_stride = stbi__flip_vertically_on_write ? -stride_bytes : stride_bytes;
+
+   if (type==0) {
+      memcpy(line_buffer, z, width*n);
+      return;
+   }
+
+   // first loop isn't optimized since it's just one pixel
+   for (i = 0; i < n; ++i) {
+      switch (type) {
+         case 1: line_buffer[i] = z[i]; break;
+         case 2: line_buffer[i] = z[i] - z[i-signed_stride]; break;
+         case 3: line_buffer[i] = z[i] - (z[i-signed_stride]>>1); break;
+         case 4: line_buffer[i] = (signed char) (z[i] - stbiw__paeth(0,z[i-signed_stride],0)); break;
+         case 5: line_buffer[i] = z[i]; break;
+         case 6: line_buffer[i] = z[i]; break;
+      }
+   }
+   switch (type) {
+      case 1: for (i=n; i < width*n; ++i) line_buffer[i] = z[i] - z[i-n]; break;
+      case 2: for (i=n; i < width*n; ++i) line_buffer[i] = z[i] - z[i-signed_stride]; break;
+      case 3: for (i=n; i < width*n; ++i) line_buffer[i] = z[i] - ((z[i-n] + z[i-signed_stride])>>1); break;
+      case 4: for (i=n; i < width*n; ++i) line_buffer[i] = z[i] - stbiw__paeth(z[i-n], z[i-signed_stride], z[i-signed_stride-n]); break;
+      case 5: for (i=n; i < width*n; ++i) line_buffer[i] = z[i] - (z[i-n]>>1); break;
+      case 6: for (i=n; i < width*n; ++i) line_buffer[i] = z[i] - stbiw__paeth(z[i-n], 0,0); break;
+   }
+}
+
+STBIWDEF unsigned char *stbi_write_png_to_mem(const unsigned char *pixels, int stride_bytes, int x, int y, int n, int *out_len)
+{
+   int force_filter = stbi_write_force_png_filter;
+   int ctype[5] = { -1, 0, 4, 2, 6 };
+   unsigned char sig[8] = { 137,80,78,71,13,10,26,10 };
+   unsigned char *out,*o, *filt, *zlib;
+   signed char *line_buffer;
+   int j,zlen;
+
+   if (stride_bytes == 0)
+      stride_bytes = x * n;
+
+   if (force_filter >= 5) {
+      force_filter = -1;
+   }
+
+   filt = (unsigned char *) STBIW_MALLOC((x*n+1) * y); if (!filt) return 0;
+   line_buffer = (signed char *) STBIW_MALLOC(x * n); if (!line_buffer) { STBIW_FREE(filt); return 0; }
+   for (j=0; j < y; ++j) {
+      int filter_type;
+      if (force_filter > -1) {
+         filter_type = force_filter;
+         stbiw__encode_png_line((unsigned char*)(pixels), stride_bytes, x, y, j, n, force_filter, line_buffer);
+      } else { // Estimate the best filter by running through all of them:
+         int best_filter = 0, best_filter_val = 0x7fffffff, est, i;
+         for (filter_type = 0; filter_type < 5; filter_type++) {
+            stbiw__encode_png_line((unsigned char*)(pixels), stride_bytes, x, y, j, n, filter_type, line_buffer);
+
+            // Estimate the entropy of the line using this filter; the less, the better.
+            est = 0;
+            for (i = 0; i < x*n; ++i) {
+               est += abs((signed char) line_buffer[i]);
+            }
+            if (est < best_filter_val) {
+               best_filter_val = est;
+               best_filter = filter_type;
+            }
+         }
+         if (filter_type != best_filter) {  // If the last iteration already got us the best filter, don't redo it
+            stbiw__encode_png_line((unsigned char*)(pixels), stride_bytes, x, y, j, n, best_filter, line_buffer);
+            filter_type = best_filter;
+         }
+      }
+      // when we get here, filter_type contains the filter type, and line_buffer contains the data
+      filt[j*(x*n+1)] = (unsigned char) filter_type;
+      STBIW_MEMMOVE(filt+j*(x*n+1)+1, line_buffer, x*n);
+   }
+   STBIW_FREE(line_buffer);
+   zlib = stbi_zlib_compress(filt, y*( x*n+1), &zlen, stbi_write_png_compression_level);
+   STBIW_FREE(filt);
+   if (!zlib) return 0;
+
+   // each tag requires 12 bytes of overhead
+   out = (unsigned char *) STBIW_MALLOC(8 + 12+13 + 12+zlen + 12);
+   if (!out) return 0;
+   *out_len = 8 + 12+13 + 12+zlen + 12;
+
+   o=out;
+   STBIW_MEMMOVE(o,sig,8); o+= 8;
+   stbiw__wp32(o, 13); // header length
+   stbiw__wptag(o, "IHDR");
+   stbiw__wp32(o, x);
+   stbiw__wp32(o, y);
+   *o++ = 8;
+   *o++ = STBIW_UCHAR(ctype[n]);
+   *o++ = 0;
+   *o++ = 0;
+   *o++ = 0;
+   stbiw__wpcrc(&o,13);
+
+   stbiw__wp32(o, zlen);
+   stbiw__wptag(o, "IDAT");
+   STBIW_MEMMOVE(o, zlib, zlen);
+   o += zlen;
+   STBIW_FREE(zlib);
+   stbiw__wpcrc(&o, zlen);
+
+   stbiw__wp32(o,0);
+   stbiw__wptag(o, "IEND");
+   stbiw__wpcrc(&o,0);
+
+   STBIW_ASSERT(o == out + *out_len);
+
+   return out;
+}
+
+#ifndef STBI_WRITE_NO_STDIO
+STBIWDEF int stbi_write_png(char const *filename, int x, int y, int comp, const void *data, int stride_bytes)
+{
+   FILE *f;
+   int len;
+   unsigned char *png = stbi_write_png_to_mem((const unsigned char *) data, stride_bytes, x, y, comp, &len);
+   if (png == NULL) return 0;
+
+   f = stbiw__fopen(filename, "wb");
+   if (!f) { STBIW_FREE(png); return 0; }
+   fwrite(png, 1, len, f);
+   fclose(f);
+   STBIW_FREE(png);
+   return 1;
+}
+#endif
+
+STBIWDEF int stbi_write_png_to_func(stbi_write_func *func, void *context, int x, int y, int comp, const void *data, int stride_bytes)
+{
+   int len;
+   unsigned char *png = stbi_write_png_to_mem((const unsigned char *) data, stride_bytes, x, y, comp, &len);
+   if (png == NULL) return 0;
+   func(context, png, len);
+   STBIW_FREE(png);
+   return 1;
+}
+
+
+/* ***************************************************************************
+ *
+ * JPEG writer
+ *
+ * This is based on Jon Olick's jo_jpeg.cpp:
+ * public domain Simple, Minimalistic JPEG writer - http://www.jonolick.com/code.html
+ */
+
+static const unsigned char stbiw__jpg_ZigZag[] = { 0,1,5,6,14,15,27,28,2,4,7,13,16,26,29,42,3,8,12,17,25,30,41,43,9,11,18,
+      24,31,40,44,53,10,19,23,32,39,45,52,54,20,22,33,38,46,51,55,60,21,34,37,47,50,56,59,61,35,36,48,49,57,58,62,63 };
+
+static void stbiw__jpg_writeBits(stbi__write_context *s, int *bitBufP, int *bitCntP, const unsigned short *bs) {
+   int bitBuf = *bitBufP, bitCnt = *bitCntP;
+   bitCnt += bs[1];
+   bitBuf |= bs[0] << (24 - bitCnt);
+   while(bitCnt >= 8) {
+      unsigned char c = (bitBuf >> 16) & 255;
+      stbiw__putc(s, c);
+      if(c == 255) {
+         stbiw__putc(s, 0);
+      }
+      bitBuf <<= 8;
+      bitCnt -= 8;
+   }
+   *bitBufP = bitBuf;
+   *bitCntP = bitCnt;
+}
+
+static void stbiw__jpg_DCT(float *d0p, float *d1p, float *d2p, float *d3p, float *d4p, float *d5p, float *d6p, float *d7p) {
+   float d0 = *d0p, d1 = *d1p, d2 = *d2p, d3 = *d3p, d4 = *d4p, d5 = *d5p, d6 = *d6p, d7 = *d7p;
+   float z1, z2, z3, z4, z5, z11, z13;
+
+   float tmp0 = d0 + d7;
+   float tmp7 = d0 - d7;
+   float tmp1 = d1 + d6;
+   float tmp6 = d1 - d6;
+   float tmp2 = d2 + d5;
+   float tmp5 = d2 - d5;
+   float tmp3 = d3 + d4;
+   float tmp4 = d3 - d4;
+
+   // Even part
+   float tmp10 = tmp0 + tmp3;   // phase 2
+   float tmp13 = tmp0 - tmp3;
+   float tmp11 = tmp1 + tmp2;
+   float tmp12 = tmp1 - tmp2;
+
+   d0 = tmp10 + tmp11;       // phase 3
+   d4 = tmp10 - tmp11;
+
+   z1 = (tmp12 + tmp13) * 0.707106781f; // c4
+   d2 = tmp13 + z1;       // phase 5
+   d6 = tmp13 - z1;
+
+   // Odd part
+   tmp10 = tmp4 + tmp5;       // phase 2
+   tmp11 = tmp5 + tmp6;
+   tmp12 = tmp6 + tmp7;
+
+   // The rotator is modified from fig 4-8 to avoid extra negations.
+   z5 = (tmp10 - tmp12) * 0.382683433f; // c6
+   z2 = tmp10 * 0.541196100f + z5; // c2-c6
+   z4 = tmp12 * 1.306562965f + z5; // c2+c6
+   z3 = tmp11 * 0.707106781f; // c4
+
+   z11 = tmp7 + z3;      // phase 5
+   z13 = tmp7 - z3;
+
+   *d5p = z13 + z2;         // phase 6
+   *d3p = z13 - z2;
+   *d1p = z11 + z4;
+   *d7p = z11 - z4;
+
+   *d0p = d0;  *d2p = d2;  *d4p = d4;  *d6p = d6;
+}
+
+static void stbiw__jpg_calcBits(int val, unsigned short bits[2]) {
+   int tmp1 = val < 0 ? -val : val;
+   val = val < 0 ? val-1 : val;
+   bits[1] = 1;
+   while(tmp1 >>= 1) {
+      ++bits[1];
+   }
+   bits[0] = val & ((1<<bits[1])-1);
+}
+
+static int stbiw__jpg_processDU(stbi__write_context *s, int *bitBuf, int *bitCnt, float *CDU, int du_stride, float *fdtbl, int DC, const unsigned short HTDC[256][2], const unsigned short HTAC[256][2]) {
+   const unsigned short EOB[2] = { HTAC[0x00][0], HTAC[0x00][1] };
+   const unsigned short M16zeroes[2] = { HTAC[0xF0][0], HTAC[0xF0][1] };
+   int dataOff, i, j, n, diff, end0pos, x, y;
+   int DU[64];
+
+   // DCT rows
+   for(dataOff=0, n=du_stride*8; dataOff<n; dataOff+=du_stride) {
+      stbiw__jpg_DCT(&CDU[dataOff], &CDU[dataOff+1], &CDU[dataOff+2], &CDU[dataOff+3], &CDU[dataOff+4], &CDU[dataOff+5], &CDU[dataOff+6], &CDU[dataOff+7]);
+   }
+   // DCT columns
+   for(dataOff=0; dataOff<8; ++dataOff) {
+      stbiw__jpg_DCT(&CDU[dataOff], &CDU[dataOff+du_stride], &CDU[dataOff+du_stride*2], &CDU[dataOff+du_stride*3], &CDU[dataOff+du_stride*4],
+                     &CDU[dataOff+du_stride*5], &CDU[dataOff+du_stride*6], &CDU[dataOff+du_stride*7]);
+   }
+   // Quantize/descale/zigzag the coefficients
+   for(y = 0, j=0; y < 8; ++y) {
+      for(x = 0; x < 8; ++x,++j) {
+         float v;
+         i = y*du_stride+x;
+         v = CDU[i]*fdtbl[j];
+         // DU[stbiw__jpg_ZigZag[j]] = (int)(v < 0 ? ceilf(v - 0.5f) : floorf(v + 0.5f));
+         // ceilf() and floorf() are C99, not C89, but I /think/ they're not needed here anyway?
+         DU[stbiw__jpg_ZigZag[j]] = (int)(v < 0 ? v - 0.5f : v + 0.5f);
+      }
+   }
+
+   // Encode DC
+   diff = DU[0] - DC;
+   if (diff == 0) {
+      stbiw__jpg_writeBits(s, bitBuf, bitCnt, HTDC[0]);
+   } else {
+      unsigned short bits[2];
+      stbiw__jpg_calcBits(diff, bits);
+      stbiw__jpg_writeBits(s, bitBuf, bitCnt, HTDC[bits[1]]);
+      stbiw__jpg_writeBits(s, bitBuf, bitCnt, bits);
+   }
+   // Encode ACs
+   end0pos = 63;
+   for(; (end0pos>0)&&(DU[end0pos]==0); --end0pos) {
+   }
+   // end0pos = first element in reverse order !=0
+   if(end0pos == 0) {
+      stbiw__jpg_writeBits(s, bitBuf, bitCnt, EOB);
+      return DU[0];
+   }
+   for(i = 1; i <= end0pos; ++i) {
+      int startpos = i;
+      int nrzeroes;
+      unsigned short bits[2];
+      for (; DU[i]==0 && i<=end0pos; ++i) {
+      }
+      nrzeroes = i-startpos;
+      if ( nrzeroes >= 16 ) {
+         int lng = nrzeroes>>4;
+         int nrmarker;
+         for (nrmarker=1; nrmarker <= lng; ++nrmarker)
+            stbiw__jpg_writeBits(s, bitBuf, bitCnt, M16zeroes);
+         nrzeroes &= 15;
+      }
+      stbiw__jpg_calcBits(DU[i], bits);
+      stbiw__jpg_writeBits(s, bitBuf, bitCnt, HTAC[(nrzeroes<<4)+bits[1]]);
+      stbiw__jpg_writeBits(s, bitBuf, bitCnt, bits);
+   }
+   if(end0pos != 63) {
+      stbiw__jpg_writeBits(s, bitBuf, bitCnt, EOB);
+   }
+   return DU[0];
+}
+
+static int stbi_write_jpg_core(stbi__write_context *s, int width, int height, int comp, const void* data, int quality) {
+   // Constants that don't pollute global namespace
+   static const unsigned char std_dc_luminance_nrcodes[] = {0,0,1,5,1,1,1,1,1,1,0,0,0,0,0,0,0};
+   static const unsigned char std_dc_luminance_values[] = {0,1,2,3,4,5,6,7,8,9,10,11};
+   static const unsigned char std_ac_luminance_nrcodes[] = {0,0,2,1,3,3,2,4,3,5,5,4,4,0,0,1,0x7d};
+   static const unsigned char std_ac_luminance_values[] = {
+      0x01,0x02,0x03,0x00,0x04,0x11,0x05,0x12,0x21,0x31,0x41,0x06,0x13,0x51,0x61,0x07,0x22,0x71,0x14,0x32,0x81,0x91,0xa1,0x08,
+      0x23,0x42,0xb1,0xc1,0x15,0x52,0xd1,0xf0,0x24,0x33,0x62,0x72,0x82,0x09,0x0a,0x16,0x17,0x18,0x19,0x1a,0x25,0x26,0x27,0x28,
+      0x29,0x2a,0x34,0x35,0x36,0x37,0x38,0x39,0x3a,0x43,0x44,0x45,0x46,0x47,0x48,0x49,0x4a,0x53,0x54,0x55,0x56,0x57,0x58,0x59,
+      0x5a,0x63,0x64,0x65,0x66,0x67,0x68,0x69,0x6a,0x73,0x74,0x75,0x76,0x77,0x78,0x79,0x7a,0x83,0x84,0x85,0x86,0x87,0x88,0x89,
+      0x8a,0x92,0x93,0x94,0x95,0x96,0x97,0x98,0x99,0x9a,0xa2,0xa3,0xa4,0xa5,0xa6,0xa7,0xa8,0xa9,0xaa,0xb2,0xb3,0xb4,0xb5,0xb6,
+      0xb7,0xb8,0xb9,0xba,0xc2,0xc3,0xc4,0xc5,0xc6,0xc7,0xc8,0xc9,0xca,0xd2,0xd3,0xd4,0xd5,0xd6,0xd7,0xd8,0xd9,0xda,0xe1,0xe2,
+      0xe3,0xe4,0xe5,0xe6,0xe7,0xe8,0xe9,0xea,0xf1,0xf2,0xf3,0xf4,0xf5,0xf6,0xf7,0xf8,0xf9,0xfa
+   };
+   static const unsigned char std_dc_chrominance_nrcodes[] = {0,0,3,1,1,1,1,1,1,1,1,1,0,0,0,0,0};
+   static const unsigned char std_dc_chrominance_values[] = {0,1,2,3,4,5,6,7,8,9,10,11};
+   static const unsigned char std_ac_chrominance_nrcodes[] = {0,0,2,1,2,4,4,3,4,7,5,4,4,0,1,2,0x77};
+   static const unsigned char std_ac_chrominance_values[] = {
+      0x00,0x01,0x02,0x03,0x11,0x04,0x05,0x21,0x31,0x06,0x12,0x41,0x51,0x07,0x61,0x71,0x13,0x22,0x32,0x81,0x08,0x14,0x42,0x91,
+      0xa1,0xb1,0xc1,0x09,0x23,0x33,0x52,0xf0,0x15,0x62,0x72,0xd1,0x0a,0x16,0x24,0x34,0xe1,0x25,0xf1,0x17,0x18,0x19,0x1a,0x26,
+      0x27,0x28,0x29,0x2a,0x35,0x36,0x37,0x38,0x39,0x3a,0x43,0x44,0x45,0x46,0x47,0x48,0x49,0x4a,0x53,0x54,0x55,0x56,0x57,0x58,
+      0x59,0x5a,0x63,0x64,0x65,0x66,0x67,0x68,0x69,0x6a,0x73,0x74,0x75,0x76,0x77,0x78,0x79,0x7a,0x82,0x83,0x84,0x85,0x86,0x87,
+      0x88,0x89,0x8a,0x92,0x93,0x94,0x95,0x96,0x97,0x98,0x99,0x9a,0xa2,0xa3,0xa4,0xa5,0xa6,0xa7,0xa8,0xa9,0xaa,0xb2,0xb3,0xb4,
+      0xb5,0xb6,0xb7,0xb8,0xb9,0xba,0xc2,0xc3,0xc4,0xc5,0xc6,0xc7,0xc8,0xc9,0xca,0xd2,0xd3,0xd4,0xd5,0xd6,0xd7,0xd8,0xd9,0xda,
+      0xe2,0xe3,0xe4,0xe5,0xe6,0xe7,0xe8,0xe9,0xea,0xf2,0xf3,0xf4,0xf5,0xf6,0xf7,0xf8,0xf9,0xfa
+   };
+   // Huffman tables
+   static const unsigned short YDC_HT[256][2] = { {0,2},{2,3},{3,3},{4,3},{5,3},{6,3},{14,4},{30,5},{62,6},{126,7},{254,8},{510,9}};
+   static const unsigned short UVDC_HT[256][2] = { {0,2},{1,2},{2,2},{6,3},{14,4},{30,5},{62,6},{126,7},{254,8},{510,9},{1022,10},{2046,11}};
+   static const unsigned short YAC_HT[256][2] = {
+      {10,4},{0,2},{1,2},{4,3},{11,4},{26,5},{120,7},{248,8},{1014,10},{65410,16},{65411,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {12,4},{27,5},{121,7},{502,9},{2038,11},{65412,16},{65413,16},{65414,16},{65415,16},{65416,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {28,5},{249,8},{1015,10},{4084,12},{65417,16},{65418,16},{65419,16},{65420,16},{65421,16},{65422,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {58,6},{503,9},{4085,12},{65423,16},{65424,16},{65425,16},{65426,16},{65427,16},{65428,16},{65429,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {59,6},{1016,10},{65430,16},{65431,16},{65432,16},{65433,16},{65434,16},{65435,16},{65436,16},{65437,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {122,7},{2039,11},{65438,16},{65439,16},{65440,16},{65441,16},{65442,16},{65443,16},{65444,16},{65445,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {123,7},{4086,12},{65446,16},{65447,16},{65448,16},{65449,16},{65450,16},{65451,16},{65452,16},{65453,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {250,8},{4087,12},{65454,16},{65455,16},{65456,16},{65457,16},{65458,16},{65459,16},{65460,16},{65461,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {504,9},{32704,15},{65462,16},{65463,16},{65464,16},{65465,16},{65466,16},{65467,16},{65468,16},{65469,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {505,9},{65470,16},{65471,16},{65472,16},{65473,16},{65474,16},{65475,16},{65476,16},{65477,16},{65478,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {506,9},{65479,16},{65480,16},{65481,16},{65482,16},{65483,16},{65484,16},{65485,16},{65486,16},{65487,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {1017,10},{65488,16},{65489,16},{65490,16},{65491,16},{65492,16},{65493,16},{65494,16},{65495,16},{65496,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {1018,10},{65497,16},{65498,16},{65499,16},{65500,16},{65501,16},{65502,16},{65503,16},{65504,16},{65505,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {2040,11},{65506,16},{65507,16},{65508,16},{65509,16},{65510,16},{65511,16},{65512,16},{65513,16},{65514,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {65515,16},{65516,16},{65517,16},{65518,16},{65519,16},{65520,16},{65521,16},{65522,16},{65523,16},{65524,16},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {2041,11},{65525,16},{65526,16},{65527,16},{65528,16},{65529,16},{65530,16},{65531,16},{65532,16},{65533,16},{65534,16},{0,0},{0,0},{0,0},{0,0},{0,0}
+   };
+   static const unsigned short UVAC_HT[256][2] = {
+      {0,2},{1,2},{4,3},{10,4},{24,5},{25,5},{56,6},{120,7},{500,9},{1014,10},{4084,12},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {11,4},{57,6},{246,8},{501,9},{2038,11},{4085,12},{65416,16},{65417,16},{65418,16},{65419,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {26,5},{247,8},{1015,10},{4086,12},{32706,15},{65420,16},{65421,16},{65422,16},{65423,16},{65424,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {27,5},{248,8},{1016,10},{4087,12},{65425,16},{65426,16},{65427,16},{65428,16},{65429,16},{65430,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {58,6},{502,9},{65431,16},{65432,16},{65433,16},{65434,16},{65435,16},{65436,16},{65437,16},{65438,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {59,6},{1017,10},{65439,16},{65440,16},{65441,16},{65442,16},{65443,16},{65444,16},{65445,16},{65446,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {121,7},{2039,11},{65447,16},{65448,16},{65449,16},{65450,16},{65451,16},{65452,16},{65453,16},{65454,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {122,7},{2040,11},{65455,16},{65456,16},{65457,16},{65458,16},{65459,16},{65460,16},{65461,16},{65462,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {249,8},{65463,16},{65464,16},{65465,16},{65466,16},{65467,16},{65468,16},{65469,16},{65470,16},{65471,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {503,9},{65472,16},{65473,16},{65474,16},{65475,16},{65476,16},{65477,16},{65478,16},{65479,16},{65480,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {504,9},{65481,16},{65482,16},{65483,16},{65484,16},{65485,16},{65486,16},{65487,16},{65488,16},{65489,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {505,9},{65490,16},{65491,16},{65492,16},{65493,16},{65494,16},{65495,16},{65496,16},{65497,16},{65498,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {506,9},{65499,16},{65500,16},{65501,16},{65502,16},{65503,16},{65504,16},{65505,16},{65506,16},{65507,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {2041,11},{65508,16},{65509,16},{65510,16},{65511,16},{65512,16},{65513,16},{65514,16},{65515,16},{65516,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {16352,14},{65517,16},{65518,16},{65519,16},{65520,16},{65521,16},{65522,16},{65523,16},{65524,16},{65525,16},{0,0},{0,0},{0,0},{0,0},{0,0},
+      {1018,10},{32707,15},{65526,16},{65527,16},{65528,16},{65529,16},{65530,16},{65531,16},{65532,16},{65533,16},{65534,16},{0,0},{0,0},{0,0},{0,0},{0,0}
+   };
+   static const int YQT[] = {16,11,10,16,24,40,51,61,12,12,14,19,26,58,60,55,14,13,16,24,40,57,69,56,14,17,22,29,51,87,80,62,18,22,
+                             37,56,68,109,103,77,24,35,55,64,81,104,113,92,49,64,78,87,103,121,120,101,72,92,95,98,112,100,103,99};
+   static const int UVQT[] = {17,18,24,47,99,99,99,99,18,21,26,66,99,99,99,99,24,26,56,99,99,99,99,99,47,66,99,99,99,99,99,99,
+                              99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99};
+   static const float aasf[] = { 1.0f * 2.828427125f, 1.387039845f * 2.828427125f, 1.306562965f * 2.828427125f, 1.175875602f * 2.828427125f,
+                                 1.0f * 2.828427125f, 0.785694958f * 2.828427125f, 0.541196100f * 2.828427125f, 0.275899379f * 2.828427125f };
+
+   int row, col, i, k, subsample;
+   float fdtbl_Y[64], fdtbl_UV[64];
+   unsigned char YTable[64], UVTable[64];
+
+   if(!data || !width || !height || comp > 4 || comp < 1) {
+      return 0;
+   }
+
+   quality = quality ? quality : 90;
+   subsample = quality <= 90 ? 1 : 0;
+   quality = quality < 1 ? 1 : quality > 100 ? 100 : quality;
+   quality = quality < 50 ? 5000 / quality : 200 - quality * 2;
+
+   for(i = 0; i < 64; ++i) {
+      int uvti, yti = (YQT[i]*quality+50)/100;
+      YTable[stbiw__jpg_ZigZag[i]] = (unsigned char) (yti < 1 ? 1 : yti > 255 ? 255 : yti);
+      uvti = (UVQT[i]*quality+50)/100;
+      UVTable[stbiw__jpg_ZigZag[i]] = (unsigned char) (uvti < 1 ? 1 : uvti > 255 ? 255 : uvti);
+   }
+
+   for(row = 0, k = 0; row < 8; ++row) {
+      for(col = 0; col < 8; ++col, ++k) {
+         fdtbl_Y[k]  = 1 / (YTable [stbiw__jpg_ZigZag[k]] * aasf[row] * aasf[col]);
+         fdtbl_UV[k] = 1 / (UVTable[stbiw__jpg_ZigZag[k]] * aasf[row] * aasf[col]);
+      }
+   }
+
+   // Write Headers
+   {
+      static const unsigned char head0[] = { 0xFF,0xD8,0xFF,0xE0,0,0x10,'J','F','I','F',0,1,1,0,0,1,0,1,0,0,0xFF,0xDB,0,0x84,0 };
+      static const unsigned char head2[] = { 0xFF,0xDA,0,0xC,3,1,0,2,0x11,3,0x11,0,0x3F,0 };
+      const unsigned char head1[] = { 0xFF,0xC0,0,0x11,8,(unsigned char)(height>>8),STBIW_UCHAR(height),(unsigned char)(width>>8),STBIW_UCHAR(width),
+                                      3,1,(unsigned char)(subsample?0x22:0x11),0,2,0x11,1,3,0x11,1,0xFF,0xC4,0x01,0xA2,0 };
+      s->func(s->context, (void*)head0, sizeof(head0));
+      s->func(s->context, (void*)YTable, sizeof(YTable));
+      stbiw__putc(s, 1);
+      s->func(s->context, UVTable, sizeof(UVTable));
+      s->func(s->context, (void*)head1, sizeof(head1));
+      s->func(s->context, (void*)(std_dc_luminance_nrcodes+1), sizeof(std_dc_luminance_nrcodes)-1);
+      s->func(s->context, (void*)std_dc_luminance_values, sizeof(std_dc_luminance_values));
+      stbiw__putc(s, 0x10); // HTYACinfo
+      s->func(s->context, (void*)(std_ac_luminance_nrcodes+1), sizeof(std_ac_luminance_nrcodes)-1);
+      s->func(s->context, (void*)std_ac_luminance_values, sizeof(std_ac_luminance_values));
+      stbiw__putc(s, 1); // HTUDCinfo
+      s->func(s->context, (void*)(std_dc_chrominance_nrcodes+1), sizeof(std_dc_chrominance_nrcodes)-1);
+      s->func(s->context, (void*)std_dc_chrominance_values, sizeof(std_dc_chrominance_values));
+      stbiw__putc(s, 0x11); // HTUACinfo
+      s->func(s->context, (void*)(std_ac_chrominance_nrcodes+1), sizeof(std_ac_chrominance_nrcodes)-1);
+      s->func(s->context, (void*)std_ac_chrominance_values, sizeof(std_ac_chrominance_values));
+      s->func(s->context, (void*)head2, sizeof(head2));
+   }
+
+   // Encode 8x8 macroblocks
+   {
+      static const unsigned short fillBits[] = {0x7F, 7};
+      int DCY=0, DCU=0, DCV=0;
+      int bitBuf=0, bitCnt=0;
+      // comp == 2 is grey+alpha (alpha is ignored)
+      int ofsG = comp > 2 ? 1 : 0, ofsB = comp > 2 ? 2 : 0;
+      const unsigned char *dataR = (const unsigned char *)data;
+      const unsigned char *dataG = dataR + ofsG;
+      const unsigned char *dataB = dataR + ofsB;
+      int x, y, pos;
+      if(subsample) {
+         for(y = 0; y < height; y += 16) {
+            for(x = 0; x < width; x += 16) {
+               float Y[256], U[256], V[256];
+               for(row = y, pos = 0; row < y+16; ++row) {
+                  // row >= height => use last input row
+                  int clamped_row = (row < height) ? row : height - 1;
+                  int base_p = (stbi__flip_vertically_on_write ? (height-1-clamped_row) : clamped_row)*width*comp;
+                  for(col = x; col < x+16; ++col, ++pos) {
+                     // if col >= width => use pixel from last input column
+                     int p = base_p + ((col < width) ? col : (width-1))*comp;
+                     float r = dataR[p], g = dataG[p], b = dataB[p];
+                     Y[pos]= +0.29900f*r + 0.58700f*g + 0.11400f*b - 128;
+                     U[pos]= -0.16874f*r - 0.33126f*g + 0.50000f*b;
+                     V[pos]= +0.50000f*r - 0.41869f*g - 0.08131f*b;
+                  }
+               }
+               DCY = stbiw__jpg_processDU(s, &bitBuf, &bitCnt, Y+0,   16, fdtbl_Y, DCY, YDC_HT, YAC_HT);
+               DCY = stbiw__jpg_processDU(s, &bitBuf, &bitCnt, Y+8,   16, fdtbl_Y, DCY, YDC_HT, YAC_HT);
+               DCY = stbiw__jpg_processDU(s, &bitBuf, &bitCnt, Y+128, 16, fdtbl_Y, DCY, YDC_HT, YAC_HT);
+               DCY = stbiw__jpg_processDU(s, &bitBuf, &bitCnt, Y+136, 16, fdtbl_Y, DCY, YDC_HT, YAC_HT);
+
+               // subsample U,V
+               {
+                  float subU[64], subV[64];
+                  int yy, xx;
+                  for(yy = 0, pos = 0; yy < 8; ++yy) {
+                     for(xx = 0; xx < 8; ++xx, ++pos) {
+                        int j = yy*32+xx*2;
+                        subU[pos] = (U[j+0] + U[j+1] + U[j+16] + U[j+17]) * 0.25f;
+                        subV[pos] = (V[j+0] + V[j+1] + V[j+16] + V[j+17]) * 0.25f;
+                     }
+                  }
+                  DCU = stbiw__jpg_processDU(s, &bitBuf, &bitCnt, subU, 8, fdtbl_UV, DCU, UVDC_HT, UVAC_HT);
+                  DCV = stbiw__jpg_processDU(s, &bitBuf, &bitCnt, subV, 8, fdtbl_UV, DCV, UVDC_HT, UVAC_HT);
+               }
+            }
+         }
+      } else {
+         for(y = 0; y < height; y += 8) {
+            for(x = 0; x < width; x += 8) {
+               float Y[64], U[64], V[64];
+               for(row = y, pos = 0; row < y+8; ++row) {
+                  // row >= height => use last input row
+                  int clamped_row = (row < height) ? row : height - 1;
+                  int base_p = (stbi__flip_vertically_on_write ? (height-1-clamped_row) : clamped_row)*width*comp;
+                  for(col = x; col < x+8; ++col, ++pos) {
+                     // if col >= width => use pixel from last input column
+                     int p = base_p + ((col < width) ? col : (width-1))*comp;
+                     float r = dataR[p], g = dataG[p], b = dataB[p];
+                     Y[pos]= +0.29900f*r + 0.58700f*g + 0.11400f*b - 128;
+                     U[pos]= -0.16874f*r - 0.33126f*g + 0.50000f*b;
+                     V[pos]= +0.50000f*r - 0.41869f*g - 0.08131f*b;
+                  }
+               }
+
+               DCY = stbiw__jpg_processDU(s, &bitBuf, &bitCnt, Y, 8, fdtbl_Y,  DCY, YDC_HT, YAC_HT);
+               DCU = stbiw__jpg_processDU(s, &bitBuf, &bitCnt, U, 8, fdtbl_UV, DCU, UVDC_HT, UVAC_HT);
+               DCV = stbiw__jpg_processDU(s, &bitBuf, &bitCnt, V, 8, fdtbl_UV, DCV, UVDC_HT, UVAC_HT);
+            }
+         }
+      }
+
+      // Do the bit alignment of the EOI marker
+      stbiw__jpg_writeBits(s, &bitBuf, &bitCnt, fillBits);
+   }
+
+   // EOI
+   stbiw__putc(s, 0xFF);
+   stbiw__putc(s, 0xD9);
+
+   return 1;
+}
+
+STBIWDEF int stbi_write_jpg_to_func(stbi_write_func *func, void *context, int x, int y, int comp, const void *data, int quality)
+{
+   stbi__write_context s = { 0 };
+   stbi__start_write_callbacks(&s, func, context);
+   return stbi_write_jpg_core(&s, x, y, comp, (void *) data, quality);
+}
+
+
+#ifndef STBI_WRITE_NO_STDIO
+STBIWDEF int stbi_write_jpg(char const *filename, int x, int y, int comp, const void *data, int quality)
+{
+   stbi__write_context s = { 0 };
+   if (stbi__start_write_file(&s,filename)) {
+      int r = stbi_write_jpg_core(&s, x, y, comp, data, quality);
+      stbi__end_write_file(&s);
+      return r;
+   } else
+      return 0;
+}
+#endif
+
+#endif // STB_IMAGE_WRITE_IMPLEMENTATION
+
+/* Revision history
+      1.14  (2020-02-02) updated JPEG writer to downsample chroma channels
+      1.13
+      1.12
+      1.11  (2019-08-11)
+
+      1.10  (2019-02-07)
+             support utf8 filenames in Windows; fix warnings and platform ifdefs
+      1.09  (2018-02-11)
+             fix typo in zlib quality API, improve STB_I_W_STATIC in C++
+      1.08  (2018-01-29)
+             add stbi__flip_vertically_on_write, external zlib, zlib quality, choose PNG filter
+      1.07  (2017-07-24)
+             doc fix
+      1.06 (2017-07-23)
+             writing JPEG (using Jon Olick's code)
+      1.05   ???
+      1.04 (2017-03-03)
+             monochrome BMP expansion
+      1.03   ???
+      1.02 (2016-04-02)
+             avoid allocating large structures on the stack
+      1.01 (2016-01-16)
+             STBIW_REALLOC_SIZED: support allocators with no realloc support
+             avoid race-condition in crc initialization
+             minor compile issues
+      1.00 (2015-09-14)
+             installable file IO function
+      0.99 (2015-09-13)
+             warning fixes; TGA rle support
+      0.98 (2015-04-08)
+             added STBIW_MALLOC, STBIW_ASSERT etc
+      0.97 (2015-01-18)
+             fixed HDR asserts, rewrote HDR rle logic
+      0.96 (2015-01-17)
+             add HDR output
+             fix monochrome BMP
+      0.95 (2014-08-17)
+		       add monochrome TGA output
+      0.94 (2014-05-31)
+             rename private functions to avoid conflicts with stb_image.h
+      0.93 (2014-05-27)
+             warning fixes
+      0.92 (2010-08-01)
+             casts to unsigned char to fix warnings
+      0.91 (2010-07-17)
+             first public release
+      0.90   first internal release
+*/
+
+/*
+------------------------------------------------------------------------------
+This software is available under 2 licenses -- choose whichever you prefer.
+------------------------------------------------------------------------------
+ALTERNATIVE A - MIT License
+Copyright (c) 2017 Sean Barrett
+Permission is hereby granted, free of charge, to any person obtaining a copy of
+this software and associated documentation files (the "Software"), to deal in
+the Software without restriction, including without limitation the rights to
+use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
+of the Software, and to permit persons to whom the Software is furnished to do
+so, subject to the following conditions:
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+------------------------------------------------------------------------------
+ALTERNATIVE B - Public Domain (www.unlicense.org)
+This is free and unencumbered software released into the public domain.
+Anyone is free to copy, modify, publish, use, compile, sell, or distribute this
+software, either in source code form or as a compiled binary, for any purpose,
+commercial or non-commercial, and by any means.
+In jurisdictions that recognize copyright laws, the author or authors of this
+software dedicate any and all copyright interest in the software to the public
+domain. We make this dedication for the benefit of the public at large and to
+the detriment of our heirs and successors. We intend this dedication to be an
+overt act of relinquishment in perpetuity of all present and future rights to
+this software under copyright law.
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
+ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+------------------------------------------------------------------------------
+*/
diff --git a/src/utils/save_image.h b/src/utils/save_image.h
new file mode 100644
index 0000000..87b5c12
--- /dev/null
+++ b/src/utils/save_image.h
@@ -0,0 +1,30 @@
+#ifndef BRAVE2_SAVE_IMAGE_H
+#define BRAVE2_SAVE_IMAGE_H
+
+
+#include <GLFW/glfw3.h>
+#include <vector>
+#include <string>
+
+#define STB_IMAGE_WRITE_IMPLEMENTATION
+
+#include "other/stb_image_write.h"
+
+
+void saveImage(const std::string &path, GLFWwindow *w) {
+    int width, height;
+    glfwGetFramebufferSize(w, &width, &height);
+    GLsizei nrChannels = 3;
+    GLsizei stride = nrChannels * width;
+    stride += (stride % 4) ? (4 - stride % 4) : 0;
+    GLsizei bufferSize = stride * height;
+    std::vector<char> buffer(bufferSize);
+    glPixelStorei(GL_PACK_ALIGNMENT, 4);
+    glReadBuffer(GL_FRONT);
+    glReadPixels(0, 0, width, height, GL_RGB, GL_UNSIGNED_BYTE, buffer.data());
+    stbi_flip_vertically_on_write(true);
+    stbi_write_bmp(path.c_str(), width, height, nrChannels, buffer.data());//, stride);
+}
+
+
+#endif //BRAVE2_SAVE_IMAGE_H
diff --git a/src/utils/util.cpp b/src/utils/util.cpp
new file mode 100644
index 0000000..440547c
--- /dev/null
+++ b/src/utils/util.cpp
@@ -0,0 +1,25 @@
+#include <GL/glew.h>
+#include <GLFW/glfw3.h>
+#include <random>
+#include "util.h"
+
+vec3 util::getRandomRGBColorAround(vec3 base, vec3 offset) {
+    float red = (float) rand() / RAND_MAX * offset.x + base.x;
+    float green = (float) rand() / RAND_MAX * offset.y + base.y;
+    float blue = (float) rand() / RAND_MAX * offset.z + base.z;
+    red /= 255;
+    green /= 255;
+    blue /= 255;
+    return vec3(red, green, blue);
+}
+
+float util::randomOffsetf(float base, float offset) {
+    return (float) rand() / RAND_MAX * offset + base;
+}
+
+float util::randomBetween(float from, float to) {
+    std::random_device rd;
+    std::mt19937 gen(rd());
+    std::uniform_real_distribution<float> dis(from, to);
+    return dis(gen);
+}
\ No newline at end of file
diff --git a/src/utils/util.h b/src/utils/util.h
new file mode 100644
index 0000000..e000354
--- /dev/null
+++ b/src/utils/util.h
@@ -0,0 +1,16 @@
+#ifndef BRAVE2_UTIL_H
+#define BRAVE2_UTIL_H
+
+#include "math.h"
+
+namespace util {
+    vec3 getRandomRGBColorAround(vec3 base, vec3 offset);
+
+    float randomOffsetf(float base, float offset);
+
+    /// random float in interval [from;to)
+    float randomBetween(float from, float to);
+}
+
+
+#endif //BRAVE2_UTIL_H
-- 
GitLab