/* -*- mia-c++ -*-
 *
 * This file is part of viewitgui - a library and program for the
 * visualization of 3D data sets. 
 *
 * Copyright (c) Leipzig, Madrid 1999-2013 Mirco Hellmann, Gert Wollny
 *
 * viewitgui is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or
 * (at your option) any later version.
 *
 * viewitgui is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with viewitgui; if not, see <http://www.gnu.org/licenses/>.
 */


#include <viewit/trajects.hh>

#include <iostream>
#include <queue>

using namespace std;
const string type_TTrajectories("TTrajectories");

TColor high_color(1.0,.5,0.0,0.9);
TColor low_color(0.0,1.0,1.,.8);

TTrajectories::TTrajectories(const string& name, P3DTransformation _field, shared_ptr<GL2DHTexture> _texture):
	TDrawable(name),
	field(_field),
	texture(_texture)
{
	toggle_transparent();
}

TTrajectories::TTrajectories(const string& name,const TVertexList& start_list,
			     P3DTransformation _field, shared_ptr<GL2DHTexture> _texture):
	TDrawable(name),
	field(_field),
	texture(_texture)
{
	TVertexList slist(start_list);
	toggle_transparent();
	generate_trace_list(&slist);
}

void TTrajectories::add(const TVertexList& start_list)
{
	TVertexList slist; 
	for (TVertexList::const_iterator i = start_list.begin(); 
	     i != start_list.end(); i++) {
		C3DFVector loc = *i;
		slist.push_back(loc + C3DFVector(-2.0,-2.0,-2.0));
		slist.push_back(loc + C3DFVector(-2.0,-2.0, 2.0));
		slist.push_back(loc + C3DFVector(-2.0, 2.0,-2.0));
		slist.push_back(loc + C3DFVector(-2.0, 2.0, 2.0));		
		slist.push_back(loc + C3DFVector( 2.0,-2.0,-2.0));
		slist.push_back(loc + C3DFVector( 2.0,-2.0, 2.0));
		slist.push_back(loc + C3DFVector( 2.0, 2.0,-2.0));
		slist.push_back(loc + C3DFVector( 2.0, 2.0, 2.0));		
	}
	generate_trace_list(&slist);
}

void TTrajectories::add(const TSpeedPixelQueue & start_list,size_t max_size, float min_distance)
{
	min_distance *= min_distance; 
	
	TSpeedPixelQueue  starts(start_list);
	TVertexList list; 
	while ( ( ! starts.empty() ) && ( list.size() < max_size ) ) {
		TSpeedPixel p = starts.top(); 
		starts.pop();
		
		bool usefull = true; 
		for (TVertexList::const_iterator i = list.begin(); 
		     i != list.end(); i++) {
			if (((*i) - p.get_loc()).norm2() < min_distance) {
				usefull = false; 
				break;
			}
		}
		if (usefull) 
			list.push_back(p.get_loc());
	}
	generate_trace_list(&list);
}

void TTrajectories::clear()
{
	flow.clear();
}

const string& TTrajectories::get_type() const
{
	return type_TTrajectories;
}


void TTrajectories::do_gl_draw(const TCamera& c)const
{
	glEnable(GL_LINE_SMOOTH);
	if (texture) {
		texture->gl_texture_matrix_on();
		glEnable(GL_TEXTURE_2D);
		texture->bind();
		
		for (TFlowfield::const_iterator i = flow.begin(); 
		     i != flow.end(); i++) {
			glBegin(GL_LINE_STRIP);
			
			for (TStreamline::const_iterator j = (*i).begin(); 
			     j != (*i).end(); j++) {
				glTexCoord4fv(&(*j).tex.x);
				glColor4fv(&(*j).col.x);
				glVertex3fv(&(*j).loc.x);
			}
			glEnd();
		}
		texture->gl_texture_matrix_off();
		glDisable(GL_TEXTURE_2D);

	}else{
		for (TFlowfield::const_iterator i = flow.begin(); 
		     i != flow.end(); i++) {
			glBegin(GL_LINE_STRIP);
			
			for (TStreamline::const_iterator j = (*i).begin(); 
			     j != (*i).end(); j++) {
				glColor4fv(&(*j).col.x);
				glVertex3fv(&(*j).loc.x);
			}
			glEnd();
		}
	}
	glDisable(GL_LINE_SMOOTH);
}


void TTrajectories::generate_trace_list(TVertexList *start)
{
	typedef std::list< C3DTrajectory > CTrajectoryList; 
	CTrajectoryList traces; 
	
	float max_tn = 0; 
	
	cerr << "have " << start->size() <<" start points\n";
	while (!start->empty()) {
		C3DTrajectory trace;
		float tn = get_trajectory_from(start->front(), 1.0, 40, 0.1, trace);
		
		if (trace.size() > 5) {
			traces.push_back(trace);
			if (tn > max_tn) 
				max_tn = tn; 
		}else{
			fprintf(stderr,"d");
		}
		start->pop_front();
	}
	
	for (auto trace = traces.begin(); trace != traces.end(); trace++) {
		TStreamline t;
		int pos = 0; 
		for (auto i = trace->begin(); i != trace->end(); i++) {
			stream_element_t p; 
			float v = i->weight / max_tn;
			p.col =  high_color * v +  low_color * (1-v);
			p.col.a *= alpha_shape[pos++];
			p.tex.x = i->tangent.x;
			p.tex.y = i->tangent.y;
			p.tex.z = i->tangent.z;
			p.tex.a = 1.0f;			
			
			
			pos %= NUM_ALPHA_STAGES; 
			p.loc = i->location;
			t.push_front(p);
		}
		flow.push_back(t);
	}
	cerr << "flow has "<< flow.size() <<" streamlines\n";
}

float TTrajectories::get_trajectory_from(const C3DFVector& Start,
				     float step, size_t sizelimit,
				     float normlimit, 
				     C3DTrajectory& trace)
{
	// Rough idea, needs some mor thinking
	C3DTracepoint Current; 
	auto direction = field->apply(Start);
	Current.location = Start;
	const C3DFVector transform_range(field->get_size()); 

	auto inside = [&transform_range](const C3DFVector& v) {
		return v.x >= 0 && v.x < transform_range.x && 
		v.y >= 0 && v.y < transform_range.y && 
		v.z >= 0 && v.z < transform_range.z;  
	}; 
	
	Current.weight = direction.norm();
	float maxTraceNorm = Current.weight;
	trace.push_back(Current);
	
	while (Current.weight > normlimit && trace.size() <= sizelimit) {
		Current.tangent = direction  /  Current.weight;

		Current.location += Current.tangent * step;
			
		if (!inside(Current.location)) 
			break;
		
		trace.push_back(Current);
		direction = field->apply(Current.location);
		
		float dnorm = direction.norm();
		if (dot(Current.tangent,  direction) < dnorm * (- 0.9))
			break;
		
		Current.weight = dnorm;
		if (maxTraceNorm < Current.weight) {
			maxTraceNorm = Current.weight;
		}
	}  
	
	direction = field->apply(Start);
	Current.location = Start;
	Current.weight = direction.norm();

	while (Current.weight > normlimit && trace.size() <= sizelimit) {
		
		Current.tangent = direction  /  (-Current.weight);

		Current.location = Current.tangent * step;

		if (!inside(Current.location)) 
			break;
		direction = field->apply(Current.location);
		trace.push_front(Current);
		
		float dnorm = direction.norm();
		if (dot(Current.tangent, direction) < dnorm * (- 0.9))
			break;
		
		Current.weight = dnorm;
		if (maxTraceNorm < Current.weight) {
			maxTraceNorm = Current.weight;
		}
	}
	return maxTraceNorm;
}

void TTrajectories::povray_output_object(FILE *f)const
{
	if (!is_visible()) 
		return; 
	
	
	fprintf(f,"#declare %s_line = object {\n",get_name().c_str());
	for (TFlowfield::const_iterator i = flow.begin(); 
	     i != flow.end(); i++) {
		
		TStreamline::const_iterator j = i->begin();
		
		while (j != i->end()) {
		
			const C3DFVector& last_loc   = (*j).loc; 
			const TColor& last_color = (*j).col;

			fprintf(f,"  cone { <%f, %f, %f>, %f, ",last_loc.x,last_loc.y,last_loc.z,last_color.a);
			j++; 
		
			
			const C3DFVector& loc   = (*j).loc; 
			const TColor& color = (*j).col;
			fprintf(f," <%f, %f, %f>, %f open ",loc.x,loc.y,loc.z,color.a);
			
			fprintf(f,"   texture { pigment {color rgb <%f, %f, %f> }}",color.x,color.y,color.z);
		}
		fprintf(f," finish { phong 0.75 }\n");
	}
	fprintf(f,"}\n");
}

void TTrajectories::povray_output_description(FILE *f)const
{
	
	cerr << get_name() << " povray_output_description not yet implemented " << endl; 
}

float TTrajectories::alpha_shape[NUM_ALPHA_STAGES] = 
        {1., .95,.9,.85,.8,.75,.7,.65,.6,.55,.5,.45,.4,.35,.3,.25};

//	{.25,.3,.35,.4,.45,.5,.55,.6,.65,.7,.75,.8,.85,.9,.95,1.};


struct TTraceStart {
	float norm; 
 	C3DFVector loc;
};

struct TArrowElement {
	C3DFVector loc;
	C3DFVector mid;
	float r2 ;
};

struct TFieldTracerData {
	P3DTransformation field_;
	float min_norm_;
	float min_dist_;
	
	typedef shared_ptr<priority_queue<TTraceStart> > PTraceStart; 
	PTraceStart start_list;
	list<TArrowElement> result_collector; 
	
	TFieldTracerData(P3DTransformation field, float min_norm, float min_dist);
	void init_lists();
	shared_ptr<TVertexList> get_start_vertexes();
	void clean(const TArrowElement& elm); 
	bool check(const C3DFVector& v)const;
	void collect_result(); 
	void suggest(const TVertexList& suggest);
};

bool operator < (const TTraceStart& a, const TTraceStart& b) {
	return a.norm < b.norm; 
}

TFieldTracer::TFieldTracer(P3DTransformation field, float min_norm, float min_dist):
	data(new TFieldTracerData(field, min_norm, min_dist))
{
}

TFieldTracer::~TFieldTracer() 
{
	delete data;
}

shared_ptr<TVertexList> TFieldTracer::get_start_vertexes()
{
	data->init_lists();
	return data->get_start_vertexes();
}

shared_ptr<TVertexList> TFieldTracer::get_start_vertexes(const TVertexList& suggest)
{
	data->suggest(suggest);
	return data->get_start_vertexes();
}

TFieldTracerData::TFieldTracerData(P3DTransformation field, float min_norm, float min_dist):
	field_(field),
	min_norm_(min_norm), 
	min_dist_(min_dist),
	start_list(new priority_queue<TTraceStart>())
{
}

void TFieldTracerData::init_lists() 
{
	auto src = field_->begin();

	
	TTraceStart ts; 
	for (ts.loc.z =0; ts.loc.z < field_->get_size().z; ++ts.loc.z)
		for (ts.loc.y =0; ts.loc.y < field_->get_size().y; ++ts.loc.y)
			for (ts.loc.x =0; ts.loc.x < field_->get_size().x; ++ts.loc.x,++src) {
				ts.norm = (*src).norm();
				if (ts.norm > min_norm_)
					start_list->push(ts);
			}
}

void TFieldTracerData::suggest(const TVertexList& suggest)
{
	auto end = suggest.end();
	
	TTraceStart ts;
	for (auto i = suggest.begin();  i != end; ++i) {
		
		ts.norm = i->norm();
		
		if ( ts.norm > min_norm_ ) {
			ts.loc = *i;
			start_list->push(ts);
		}else {
			cerr << ts.norm << endl; 
		}
	}
}

shared_ptr<TVertexList> TFieldTracerData::get_start_vertexes()
{

	collect_result();
		
	shared_ptr<TVertexList> result(new TVertexList);

	auto i = result_collector.begin();
	auto end = result_collector.end();
	
	while (i != end) {
		result->push_back(i->loc);
		++i; 
	}
	return result;
}

void TFieldTracerData::collect_result()
{
	int rejects = 0; 
	while ( ! start_list->empty() ) {
		// get next starting candidate
		TTraceStart ts = start_list->top();
		start_list->pop();
		
		cerr << result_collector.size() << " of " <<
			start_list->size() << "\r"; 
		// check whether we will really use it
		
		if ( ts.norm < min_norm_ ) {
			rejects++; 
			continue; 
		}
		
		TArrowElement elm; 
		C3DFVector dir = field_->apply(ts.loc);
		elm.loc = ts.loc; 
		elm.mid = ts.loc -  dir * 0.5f; 
		elm.r2 = ( dir.norm() + min_dist_ ); 
		elm.r2 *= elm.r2; 
		
		// use it.
		clean(elm);

		// add starting and end point to list
		TTraceStart l;
		l.loc = ts.loc - (2.0f * min_dist_ / ts.norm) * dir;
		if (check(l.loc)) {
			l.norm = field_->apply(l.loc).norm();
			if ( l.norm > min_norm_ )
				start_list->push(l);
		}
		
		
		l.loc = ts.loc + (1.0f + 2.0f * min_dist_ / ts.norm) * dir;
		if (check(l.loc)) {
			l.norm = field_->apply(l.loc).norm();
			if ( l.norm > min_norm_ )
				start_list->push(l);
		}
		
		result_collector.push_back(elm); 
	}
}

template <class T> void swap(T *a, T *b)
{
	T help = *a; 
	*a = *b; 
	*b = help; 
}

void TFieldTracerData::clean(const TArrowElement& elm)
{
	PTraceStart help(new priority_queue<TTraceStart>()); 
	while ( ! start_list->empty() ) {
		// get next starting candidate
		TTraceStart ts = start_list->top();
		start_list->pop();
		
		if ( (ts.loc - elm.mid).norm2() > elm.r2)
			help->push(ts);
	}		
	start_list = help; 
}

bool TFieldTracerData::check(const C3DFVector& v)const 
{
	auto i = result_collector.begin();
	auto end = result_collector.end();
	
	while (i != end) {
		if ( ( i->mid - v).norm2() < i->r2)
			return false; 
		++i; 
	}
	return true; 
}