Commit ae14fe01 authored by Sabrina Templeton's avatar Sabrina Templeton
Browse files

added in nearest neighbors calculation and shading for square diagram

parent c69b572f
......@@ -20,10 +20,20 @@ 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 ) {
u = random( vec2(v_Index,v_Index) ) * (u_nb_points+1);
float umin = 0.;//u_umin;
float umax = u_nb_points +1;//u_umax;
......
......@@ -5,41 +5,15 @@
layout (points) in;
precision highp float;
// uniform mat4 u_ModelViewProjectionMatrix;
// uniform mat4 u_NormalMatrix;
// uniform mat4 u_ModelViewMatrix;
uniform int u_nb_neighbors;
// 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
uniform int u_nb_neighbors;
// 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;
const int MAX_VERTS = 30;
// out vec3 v_Position;
// out vec3 v_Normal;
// out vec3 v_Parameter;
flat in int[] instance_ID;
out float v_Index;
......@@ -59,15 +33,6 @@ void render_poly(vec3[MAX_VERTS] poly, vec3 center, int space){
//vec3 p1;
for (int i = 0; i < space -2; i++){
// if (i == space-1){
// p0 = poly[space -1];
// p1 = poly[0];
// }
// else{
// p0 = poly[i];
// p1 = poly[i+1];
// }
gl_Position = u_ModelViewProjectionMatrix*vec4(p0, 1.0);
v_Index = float(instance_ID[0]);
......@@ -88,54 +53,6 @@ void render_poly(vec3[MAX_VERTS] poly, vec3 center, int space){
}
}
//render "points" as squares
void render_square(vec3 pt, float size) {
//define all four points of the square
vec3 bl = pt + vec3(-size, -size, 0);
vec3 br = pt + vec3(size, -size, 0);
vec3 tl = pt + vec3(-size, size, 0);
vec3 tr = pt + vec3(size, size, 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 = instance_ID[0];
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 = instance_ID[0];
EmitVertex();
EndPrimitive();
}
float dist(vec3 a, vec3 b){
//return length(b -a);
float comp1 = pow(b.x - a.x, 2);
float comp2 = pow(b.y - a.y, 2);
return comp1 + comp2;
}
vec3 calc_intersect(vec3 p1, vec3 p2, vec3 q1, vec3 q2){ // p1 and p2 are seeds, q1 and q2 are borders
//initialize variables
float d = 0;
......@@ -186,13 +103,6 @@ void calc_poly(vec3 zi, vec3 zj, inout vec3[MAX_VERTS] curr_poly, inout int cspa
vec3 new_poly[MAX_VERTS];//= vec3[4](vec3(0, 0, 0), vec3(0, 0, 0), vec3(0, 0, 0), vec3(0, 0, 0));
for (int i = 0; i < cspace; i++) {
//curr_poly[i] = vec3( 0.1*i , 0, 0.0 );
//render_square(new_poly[i],0.05);
//continue;
vec3 b0;
vec3 b1;// initialize b0 and b1
if (i == cspace-1){
......@@ -239,20 +149,12 @@ void calc_poly(vec3 zi, vec3 zj, inout vec3[MAX_VERTS] curr_poly, inout int cspa
}
//otherwise both sides are -1 and we take neither
}
// note to self: this is definitely outside of the loop ( and should remain to be)
curr_poly = new_poly; // does this need to be a deep copy? I don't think so.
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 glPrimitiveIDIn
......@@ -270,11 +172,7 @@ main() {
vec3 curr_poly3 = vec3(1, 1, 0);
vec3 curr_poly4 = vec3(0, 1, 0);
//render it here
// render_square(curr_poly1, .02);
// render_square(curr_poly2, .02);
// render_square(curr_poly3, .02);
// render_square(curr_poly4, .02);
vec3 poly[MAX_VERTS];
vec3[MAX_VERTS] curr_poly;
......@@ -285,39 +183,6 @@ main() {
curr_poly[3] = curr_poly4;
int space = 4; //this is cspace
//************************************************************************************************
//This is a test of the the intersect function. It works!!
// vec3 test_pt = calc_intersect(vec3(.3, .5, 0), vec3(.3, -.5, 0), vec3(0, .5, 0), vec3(0, -.5, 0));
// render_square(test_pt, .05);
//************************************************************************************************
//nearest neighbor hardcoding
// int neighbor_hardcoded[49] = int[49]( 0, 1, 2, 3, 4, 5, 6,
// 1, 0, 4, 3, 2, 5, 6,
// 2, 0, 1, 3, 4, 5, 6,
// 3, 0, 1, 2, 4, 5, 6,
// 4, 1, 0, 2, 3, 5, 6,
// 5, 6, 3, 2, 0, 1, 4,
// 6, 5, 3, 2, 0, 1, 4 );
//nearest neighbor test
// for (int i = 1; i < 49; i ++){
// int nn_test = int(texelFetch( nn, i).x);
// // if (nn_test == neighbor_hardcoded[i]) {
// // render_square(vec3(-.5, -.5, 0), .05);
// // }
// if (nn_test > 1) {
// render_square(vec3(-.8, -.8, 0), .05);
// }
// }
// this will obviously change once we are actually fetching from the nearest neighbor instead of the seed
float d = 0;
for (int i = 1; i < u_nb_neighbors; i ++){
......@@ -331,12 +196,6 @@ main() {
float d1 = d;
d = distance(nearest_neighbor, p0);
// This is a test to make sure that nearest neighbors are sorted in increaasing order. it passes
// if (d1 > d) {
// render_square(vec3(0, 1, 0), .5);
// }
calc_poly(p0, nearest_neighbor, curr_poly, space);
//calculate radius of security
......@@ -351,31 +210,13 @@ main() {
furthest_pt = curr_pt;
}
}
//render_square(furthest_pt, .02);
// if (site_idx == 0){
// render_square(furthest_pt, i * .001);
// }
// if (radius != distance(furthest_pt, p0)){
// render_square(vec3(0, 0, 0), .5);
// };
if (distance(nearest_neighbor, p0) > (2.1 * radius)){
//render_square(vec3(0, 0, 0), .1);
break;
}
}
//if (site_idx == 0) {
render_poly(curr_poly, p0, space);
//render_square(vec3(1, 1, 0), .01);
//}
//render_square(p0, .03);
render_poly(curr_poly, p0, space);
}
......
......@@ -24,6 +24,10 @@
#include <fstream>
#include "common/process.h"
#include "common/parallel_for.h"
using namespace avro;
UT_TEST_SUITE( voronoi_toy )
......@@ -47,6 +51,39 @@ 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
......@@ -81,17 +118,16 @@ UT_TEST_CASE( test1 )
// graphics::ShaderProgram points_shader("points",false,{"#version 410"});
// points_shader.use();
graphics::gl_index vertex_array;
GL_CALL( glGenVertexArrays( 1, &vertex_array ) );
GL_CALL( glBindVertexArray(vertex_array) );
//generate random seeds desired times 3
int num_points = 1e6;
int num_points = 1000;
int size = num_points*3;
int DEBUG = 0;
int DEBUG = 1;
std::vector<graphics::gl_float> seeds;
......@@ -127,6 +163,7 @@ else{
}
}
num_points = seeds.size()/3;
}
......@@ -173,7 +210,7 @@ else{
nns->set_points( num_points , seeds_d.data() );
nns->set_points( num_points , seeds_d.data() );
index_t nb_neighbors = 30; // should be 50 or 30
std::vector<index_t> neighbors(nb_neighbors);
......@@ -201,8 +238,16 @@ else{
}
}
NearestNeighboursCalculator nnc( 3 , seeds_d , nb_neighbors );
nnc.calculate();
std::vector<GLuint>& neighbors_f = nnc.neighbours;
clock_t TIME1 = clock();
printf("--> nearest neighbor calculation time = %g seconds\n",real_t(TIME1-TIME0)/real_t(CLOCKS_PER_SEC));
int nb_thread = ProcessCPU::maximum_concurrent_threads();
printf("--> nearest neighbor calculation time = %g seconds\n",real_t(TIME1-TIME0)/real_t(CLOCKS_PER_SEC)/nb_thread);
// std::vector<graphics::gl_float> seeds = {
......@@ -227,7 +272,7 @@ else{
// };
// convert nearest neighbor type
std::vector<GLuint> neighbors_f( total_neighbors.begin() , total_neighbors.end() );
//std::vector<GLuint> neighbors_f( total_neighbors.begin() , total_neighbors.end() );
graphics::gl_index points_buffer;
......@@ -360,8 +405,9 @@ else{
// call some OpenGL commands
//code borrowed from https://www.lighthouse3d.com/tutorials/opengl-timer-query/
//glEnable(GL_RASTERIZER_DISCARD);
GL_CALL( glDrawArrays(GL_POINTS, 0, num_points) );
//glDisable(GL_RASTERIZER_DISCARD);
glEndQuery(GL_TIME_ELAPSED);
GLuint64 time;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment