Commit 89d54913 authored by philip's avatar philip
Browse files

Merge branch 'sabrina' into 'main'

merge sabrina's branch for computing and visualizing voronoi diagrams with geometry shader

See merge request philipclaude/avro!88
parents 81004560 1623df6e
......@@ -89,6 +89,12 @@ ShaderProgram::ShaderProgram( const std::string& name , bool with_tess , const s
std::string geo_src = get_shader_src( base + "-geo.glsl" );
avro_assert_msg( compile(name_.c_str(),vtx_src,frg_src,geo_src) , "error compiling raytracer shader" );
}
else if (name == "voronoi") {
std::string vtx_src = get_shader_src( base + "-vtx.glsl" );
std::string frg_src = get_shader_src( base + "-frg.glsl" );
std::string geo_src = get_shader_src( base + "-geo.glsl" );
avro_assert_msg( compile(name_.c_str(),vtx_src,frg_src,geo_src) , "error compiling voronoi shader" );
}
else {
printf("unknown shader %s\n",name.c_str());
avro_assert_not_reached;
......
// the version should be specified by a macro at compile time
// e.g. #version 410 or #version 330
#version 410
layout( location = 0 ) out vec4 fragColor;
in vec3 v_Position;
in vec3 v_Normal;
in vec3 v_Parameter;
const int ncolor = 256;
uniform float u_umin;
uniform float u_umax;
uniform samplerBuffer colormap;
#define DEBUG 0
void
get_color( float u , out vec3 color ) {
float umin = 0.;//u_umin;
float umax = 7.;//u_umax;
// int indx = int(ncolor*(u - umin)/(umax - umin));
int indx = int(ncolor*(u - umin)/(umax - umin));
if (indx < 0) indx = 0;
if (indx > 255) indx = 255;
float r0 = texelFetch( colormap , 3*(indx) + 0 ).x;
float g0 = texelFetch( colormap , 3*(indx) + 1 ).x;
float b0 = texelFetch( colormap , 3*(indx) + 2 ).x;
color = vec3(r0,g0,b0);
//if (indx == 36) color = vec3(1,0,1);
}
void main() {
vec3 color = vec3(.8, .8, .8);//vec3(0.8,0.8,0.2);
int idx = gl_PrimitiveID;
float f = idx;
#if DEBUG
if (idx == 0) color = vec3(1,0,0);
if (idx == 1) color = vec3(0,1,0);
if (idx == 2) color = vec3(0,0,1);
if (idx == 3) color = vec3(1,1,0);
if (idx == 4) color = vec3(0,1,1);
if (idx == 5) color = vec3(1,0,1);
if (idx == 6) color = vec3(.5,.5,.5);
#else
get_color(idx, color);
#endif
fragColor = vec4(color, 1);
}
\ No newline at end of file
// the version should be specified by a macro at compile time
// e.g. #version 410 or #version 330
#version 410
layout (points) in;
// uniform mat4 u_ModelViewProjectionMatrix;
// uniform mat4 u_NormalMatrix;
// uniform mat4 u_ModelViewMatrix;
// out vec3 v_Position;
// out vec3 v_Normal;
// out vec3 v_Parameter;
// layout (triangle_strip, max_vertices = 18) out;
// uniform usamplerBuffer connectivity;
// uniform samplerBuffer coordinates;
//test
// the version should be specified by a macro at compile time
// e.g. #version 410 or #version 330
// the version should be specified by a macro at compile time
// e.g. #version 410 or #version 330
// layout (points) in;
uniform mat4 u_ModelViewProjectionMatrix;
// uniform mat4 u_NormalMatrix;
// uniform mat4 u_ModelViewMatrix;
// uniform int u_clip;
// uniform vec3 u_clip_center;
// uniform vec3 u_clip_normal;
// uniform int u_render_full;
// out vec3 v_Position;
// out vec3 v_Normal;
// out vec3 v_Parameter;
layout (triangle_strip, max_vertices = 16) out;
//uniform usamplerBuffer connectivity;
uniform samplerBuffer seeds;
// render full square
void render_square(vec3 pt) {
//define all four points of the square
vec3 bl = pt + vec3(-.02, -.02, 0);
vec3 br = pt + vec3(.02, -.02, 0);
vec3 tl = pt + vec3(-.02, .02, 0);
vec3 tr = pt + vec3(.02, .02, 0);
//top triangle
gl_Position = u_ModelViewProjectionMatrix*vec4(bl, 1.0);
EmitVertex();
gl_Position = u_ModelViewProjectionMatrix*vec4(tl, 1.0);
EmitVertex();
gl_Position = u_ModelViewProjectionMatrix*vec4(tr, 1.0);
//EmitVertex();
gl_PrimitiveID = gl_PrimitiveIDIn;
EmitVertex();
EndPrimitive();
gl_Position = u_ModelViewProjectionMatrix*vec4(bl, 1.0);
EmitVertex();
gl_Position = u_ModelViewProjectionMatrix*vec4(tr, 1.0);
EmitVertex();
gl_Position = u_ModelViewProjectionMatrix*vec4(br, 1.0);
//EmitVertex();
gl_PrimitiveID = gl_PrimitiveIDIn;
EmitVertex();
EndPrimitive();
}
void
main() {
// lookup the seed points, taking in glPrimitiveIDIn
vec3 p0 = vec3(texelFetch( seeds , gl_PrimitiveIDIn ).xy, 0);
render_square(p0);
}
// the version should be specified by a macro at compile time
// e.g. #version 410 or #version 330
#version 410
// this code taken from example of tetrahedron/volume toy by Philip
layout (location = 0 ) in vec3 a_Position;
// nothing is done here because the geometry shader produces and processes the actual vertices
// but we still need to set gl_Position
void main() {
gl_Position = vec4(0.0,0.0,0.0,0.0);
}
// the version should be specified by a macro at compile time
// e.g. #version 410 or #version 330
//#version 410
layout( location = 0 ) out vec4 fragColor;
in vec3 v_Position;
in vec3 v_Normal;
in vec3 v_Parameter;
in float v_Index;
const int ncolor = 256;
uniform float u_umin;
uniform float u_umax;
uniform int u_nb_points;
uniform samplerBuffer colormap;
#define DEBUG 0
float random (vec2 uv) {
return fract(sin(dot(uv,vec2(12.9898,78.233)))*43758.5453123);
}
void
get_color( float u , out vec3 color ) {
// This line turns on the random shading (for structured points)
//u = random( vec2(v_Index,v_Index) ) * (u_nb_points+1);
float umin = 0.;
float umax = u_nb_points +1;
// int indx = int(ncolor*(u - umin)/(umax - umin));
int indx = int(ncolor*(u - umin)/(umax - umin));
if (indx < 0) indx = 0;
if (indx > 255) indx = 255;
float r0 = texelFetch( colormap , 3*(indx) + 0 ).x;
float g0 = texelFetch( colormap , 3*(indx) + 1 ).x;
float b0 = texelFetch( colormap , 3*(indx) + 2 ).x;
color = vec3(r0,g0,b0);
}
void main() {
vec3 color = vec3(.2, .2, .2);
float idx = v_Index; //gl_PrimitiveID;
get_color(idx, color);
fragColor = vec4(color, 1);
}
\ No newline at end of file
// the version should be specified by a macro at compile time
// e.g. #version 410 or #version 330
#extension GL_ARB_gpu_shader_fp64 : enable
layout (points) in;
precision highp float;
uniform int u_nb_neighbors;
uniform mat4 u_ModelViewProjectionMatrix;
const int MAX_VERTS = 30;
flat in int[] instance_ID;
out float v_Index;
layout (triangle_strip, max_vertices = 40) out;
uniform samplerBuffer seeds;
uniform usamplerBuffer nn;
void render_poly(vec3[MAX_VERTS] poly, vec3 center, int space){
vec3 p0 = poly[0];
for (int i = 0; i < space -2; i++){
gl_Position = u_ModelViewProjectionMatrix*vec4(p0, 1.0);
v_Index = float(instance_ID[0]);
EmitVertex();
gl_Position = u_ModelViewProjectionMatrix*vec4(poly[i+1], 1.0);
v_Index = float(instance_ID[0]);
EmitVertex();
gl_Position = u_ModelViewProjectionMatrix*vec4(poly[i+2], 1.0);
//gl_PrimitiveID = instance_ID[0];
v_Index = float(instance_ID[0]);
EmitVertex();
EndPrimitive();
}
}
vec3 calc_intersect(vec3 p1, vec3 p2, vec3 q1, vec3 q2){ // p1 and p2 are seeds, q1 and q2 are borders
/*
* Originally from geogram (https://github.com/BrunoLevy/geogram)
* in src/lib/geogram/voronoi/generic_RVD_vertex.h
*/
//initialize variables
float d = 0;
float l1 = 0;
float l2 = 0;
float n = 0;
//deal with the x
n = p1.x - p2.x;
d -= n* (p2.x + p1.x);
l1 += q2.x * n;
l2 += q1.x * n;
//deal with the y
n = p1.y - p2.y;
d -= n* (p2.y + p1.y);
l1 += q2.y * n;
l2 += q1.y * n;
//
d = .5 * d;
l1 = abs(l1 + d);
l2 = abs(l2 + d);
float l12 = l1 + l2;
if (l12 > 1e-30) { // was 1e-30, changed to 1e-6
l1 /= l12;
l2 /= l12;
}
else{
l1 = .5;
l2 = .5;
}
float x0 = l1 * q1.x + l2 * q2.x;
float y0 = l1 * q1.y + l2 * q2.y;
return vec3(x0, y0, 0);
}
void calc_poly(vec3 zi, vec3 zj, inout vec3[MAX_VERTS] curr_poly, inout int cspace){
//cspace is an argument, which is inout as it goes along with curr_poly and represents the actual length of the polygon information in curr_poly at any given moment (ideally)
int nspace = 0; // nspace is 'new space,' the variable which acts as a pointer to keep track of the next available spot to add a vertex to new_poly
//where zi is the current site pt and zj is the current nearest neighbor
vec3 new_poly[MAX_VERTS];
for (int i = 0; i < cspace; i++) {
vec3 b0;
vec3 b1;
if (i == cspace-1){
b0 = curr_poly[cspace -1];
b1 = curr_poly[0];
}
else{
b0 = curr_poly[i];
b1 = curr_poly[i+1];
}
int side1 = 0;
int side2 = 0;
//calculate side 1: determine if the first border point is inside or outside
if (distance(b0, zi) < distance(b0, zj)) {
side1 = 1;
}
else {
side1 = -1;
}
//calculate side 2: determine if the second border point is inside or outside
if (distance(b1, zi) < distance(b1, zj)) {
side2 = 1;
}
else {
side2 = -1;
}
vec3 intersect;
if (side1 != side2){
intersect = calc_intersect(zi, zj, b0, b1); //b0 and b1 are the current border points
if (side1 == 1){
new_poly[nspace] = b0;
new_poly[nspace+1] = intersect;
nspace += 2;
}
else{
new_poly[nspace] = intersect;
nspace += 1;
}
}
else if (side1 == 1){ // meaning both sides are inside
new_poly[nspace] = b0; // we add in the first border point since we will add in the next point at the next step of the loop
nspace += 1;
}
// otherwise both sides are -1 and we take neither
}
curr_poly = new_poly;
cspace = nspace; // setting the current space to be the new space so that it is updated for the new polygon.
}
void
main() {
// Lookup the seed points, taking in instance_ID which we output from the vertex shader. This should be equivalent to glPrimitiveID_In but is not
int site_idx = instance_ID[0];
vec3 zi = vec3(texelFetch( seeds , site_idx ).xy, 0);
// This is the initial polygon, we start by initializing it to the polygon which covers the total area of the diagram
vec3[MAX_VERTS] curr_poly;
curr_poly[0] = vec3(0, 0, 0);
curr_poly[1] = vec3(1, 0, 0);
curr_poly[2] = vec3(1, 1, 0);
curr_poly[3] = vec3(0, 1, 0);
int space = 4;
for (int i = 1; i < u_nb_neighbors; i ++){
int nn0 = int(texelFetch( nn, i + (site_idx * u_nb_neighbors)).x); // i + (site_idx * 7)
vec3 zj= vec3(texelFetch( seeds, nn0).xy, 0); // nn0 is just the index so we need to look up the corresponding point in seeds
calc_poly(zi, zj, curr_poly, space);
// Below is the radius of security theorem
float radius = 0;
// loop over the vertices of the polygon to find the furthest one, the distance from zi to that point becomes the radius.
for (int j = 0; j < space; j++ ) {
vec3 curr_pt = curr_poly[j];
if(distance(curr_pt, zi) > radius) {
radius = distance(curr_pt, zi);
}
}
// If the distance between the current nearest neighbor and the current site is greater than twice the radius, stop clipping
if (distance(zj, zi) > (2.1 * radius)){
break;
}
}
render_poly(curr_poly, zi, space);
}
// the version should be specified by a macro at compile time
// e.g. #version 410 or #version 330
//#version 410
// this code taken from example of tetrahedron/volume toy by Philip
flat out int instance_ID;
layout (location = 0 ) in vec3 a_Position;
// nothing is done here because the geometry shader produces and processes the actual vertices
// but we still need to set gl_Position
void main() {
gl_Position = vec4(0.0,0.0,0.0,0.0);
instance_ID = gl_VertexID;
}
......@@ -93,11 +93,13 @@ PowerDiagram::compute() {
// compute the power diagram
//printf("computing the power diagram in %ud\n",sites_.dim());
clock_t t0 = clock();
index_t nb_thread = 1;
#if PARALLEL
typedef PowerDiagram thisclass;
ProcessCPU::parallel_for(
parallel_for_member_callback( this , &thisclass::compute ), 0,cell_.size()
);
nb_thread = ProcessCPU::maximum_concurrent_threads();
#else
for (index_t k = 0; k < cell_.size(); k++) {
compute(k);
......@@ -105,7 +107,7 @@ PowerDiagram::compute() {
#endif
// record the time it took to compute the voronoi diagram
time_voronoi_ = real_t(clock()-t0)/real_t(CLOCKS_PER_SEC)/ProcessCPU::maximum_concurrent_threads();
time_voronoi_ = real_t(clock()-t0)/real_t(CLOCKS_PER_SEC)/nb_thread;
}
......
......@@ -17,6 +17,7 @@ add_test_files( TEST_FILES "ut" api )
# sandbox files should have _toy.cpp
add_test_files( SBX_FILES "toy" sandbox/philip )
add_test_files( SBX_FILES "toy" sandbox/hiro )
add_test_files( SBX_FILES "toy" sandbox/sabrina )
set( unit_skip mesh/delaunay/vertex_ut.cpp )
......
#include "unit_tester.hpp"
#include "graphics/application.h"
#include "library/ckf.h"
#include "avro_config.h" // we need AVRO_SOURCE_DIR to find shaders
// the following two header files are only part of the avro-extensions repository
// #include "gl.h"
#include "graphics/shader.h"
#include "graphics/colormap.h"
#include "common/tools.h"
#include <geogram/nn_search.h>
#include <time.h>
#include <iostream>
#include <fstream>
#include "common/process.h"
#include "common/parallel_for.h"
using namespace avro;
UT_TEST_SUITE( voronoi_toy )
// this is all borrowed from geometry toy
// set up the view
int width = 500, height = width;
graphics::vec3 eye = {0.5,0.5,1.5};
graphics::vec3 up = {0,1,0};
graphics::vec3 center = {0.5,0.5,0.0};
float fov = M_PI/4.0;
float aspect = width/height;
float znear = 0.001;
float zfar = 100;
// define the transformation matrices, and compute MVP matrix
graphics::mat4 perspective_matrix = graphics::glm::perspective( fov , aspect , znear , zfar );
graphics::mat4 view_matrix = graphics::glm::lookAt( eye , center , up );
graphics::mat4 model_matrix = graphics::glm::identity();
graphics::mat4 mv = view_matrix * model_matrix;
graphics::mat4 mvp = perspective_matrix * mv;
struct NearestNeighboursCalculator {
typedef NearestNeighboursCalculator thisclass;
NearestNeighboursCalculator( coord_t d , const std::vector<real_t>& coords , index_t max_n ) :
dim(d),
nb_sites( coords.size() / dim ),
max_neighbours(max_n),
sites(coords),
neighbours(max_n*nb_sites) {
search = GEO::NearestNeighborSearch::create(dim,"ANN");
search->set_points( nb_sites , sites.data() );
}
~NearestNeighboursCalculator() {
delete search;
}
void compute( index_t k ) {
std::vector<double> distance2( max_neighbours ); // +1 because first neighbour is the point itself
std::vector<index_t> neighbours0( max_neighbours );
search->get_nearest_neighbors( max_neighbours, &sites[k*dim] , neighbours0.data() , distance2.data() );
for (int i = 0; i < max_neighbours; i++)
neighbours[k*max_neighbours+i] = neighbours0[i];
}
void calculate() {
ProcessCPU::parallel_for( parallel_for_member_callback( this , &thisclass::compute ) , 0 , nb_sites );
}
coord_t dim;
index_t nb_sites;
index_t max_neighbours;
const std::vector<real_t>& sites;
GEO::NearestNeighborSearch* search = nullptr;
std::vector<GLuint> neighbours;
};
UT_TEST_CASE( test1 )
{
//also taken from geometry toy
// initialize OpenGL
avro_assert_msg( glfwInit() , "problem initializing OpenGL!" );
glfwWindowHint(GLFW_CONTEXT_VERSION_MAJOR, 4);
glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 1);
glfwWindowHint(GLFW_OPENGL_FORWARD_COMPAT, GL_TRUE);
glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE);
GLFWwindow* window = glfwCreateWindow( width , height , "voronoi" , NULL, NULL);
if (!window) {
const char* description;
int code = glfwGetError(&description);
if (description)
printf("GLFW error (%d): %s\n",code,description);
glfwTerminate();
}
glfwMakeContextCurrent(window);
// load GL functions
gladLoadGL();
graphics::dumpGLInfo(); // this will be dumpGLInfo() in avro (it isn't necessary but prints some info about OpenGL, GLSL and the GPU)
//initialize the shader:
graphics::ShaderProgram shader("voronoi",false,{"#version 410"});