//vodafone call visualisation
//
//Tim Redfern Jan 2012
//
//scale of wall 10.57x2.8m - some parts are obscured
//
//416 x 110 in. ~ 62400x16500 ~ 1029 mpx ~ 2945 MB ram to open
//
//CousinMarriageWorld.svg, nominally 940 × 470 pixels, file size: 1.99 MB
//World_map_(Miller_cylindrical_projection,_blank).svg‎, nominally 634 × 477 pixels, file size: 1.84 MB
//wikimedia
//
//http://www.vectortemplates.com/vector-world-map.php

//requirements - 
//must be able to transform points to a new projection
// maybe not through normal shape library
//must be able to find a random point within the shape
// 

//250112 so far so good, can iterate children.
//safest thing is to establish that I can place a point as being inside or outside of a shape
//(are they closed properly!)

//If we can find points in them then we can start the main task:
//connecting parts of the SVG with rows in the spreadsheet

//maybe next first I should look at projections (to keep them happy)
//I have a map with rectangular coords,
//is easy to enough to play with the projection

// public boolean contains(RGeomElem shp)
// public geomerative.RRectangle getBounds()
// --> this is kind of hidden. why? anyway, it should be possible to use it
//another option is to use the OUTLINE of the country (as a purely geometric thing).

/*
 http://en.wikipedia.org/wiki/3D_projection
 
 a- the 3D position of a point A that is to be projected.
 c- the 3D position of a point C representing the camera.
 O- The orientation of the camera (represented, for instance, by Tait–Bryan angles).
 e- the viewer's position relative to the display surface.[1]


http://williams.best.vwh.net/avform.htm#Crs
==good stuff

Intermediate points on a great circle
Here we find points (lat,lon) a given fraction of the distance (d) between them.
Suppose the starting point is (lat1,lon1) and the final point (lat2,lon2) and we want the point
a fraction f along the great circle route. f=0 is point 1. f=1 is point 2. The two points cannot
be antipodal ( i.e. lat1+lat2=0 and abs(lon1-lon2)=pi) because then the route is undefined.
The intermediate latitude and longitude is then given by:

        A=sin((1-f)*d)/sin(d)
        B=sin(f*d)/sin(d)
        x = A*cos(lat1)*cos(lon1) +  B*cos(lat2)*cos(lon2)
        y = A*cos(lat1)*sin(lon1) +  B*cos(lat2)*sin(lon2)
        z = A*sin(lat1)           +  B*sin(lat2)
        lat=atan2(z,sqrt(x^2+y^2))
        lon=atan2(y,x)
        
        //find 
        
        RPoint sp,ep,mp;
        float f=0.8; //fraction along path;
        float d=acos(sin(sp.y)*sin(ep.y)+cos(sp.y)*cos(ep.y)*cos(sp.x-ep.x));
        float A=sin((1-f)*d)/sin(d);
        float B=sin(f*d)/sin(d);
        float x = A*cos(lat1)*cos(lon1) +  B*cos(lat2)*cos(lon2);
        float y = A*cos(lat1)*sin(lon1) +  B*cos(lat2)*sin(lon2);
        float z = A*sin(lat1)           +  B*sin(lat2);
        mp.y=atan2(z,sqrt(x^2+y^2));
        mp.x=atan2(y,x);
        
        
        
        drawing - circles (fixed diameter), closed shape made of 2 bezier curves

paths above 3d globe
fix missing countries
render above map
colours
sort lines according to destination ?

match 3D map or render texture out of processing
or export 3d lines - dx

cache points!

 */

import processing.pdf.*;
import geomerative.*;

float log10 (float x) {
  return (log(x) / log(10));
}
class pointNormalise {
  //take pixel coords and turn into lat/lng radians
  float xo,xs,yo,ys;
  pointNormalise(float _xo,float _xs,float _yo,float _ys) {
    xo=_xo;
    xs=_xs;
    yo=_yo;
    ys=_ys;
  }
  RPoint alise(RPoint p) {
        //normalise
        float px=(((p.x-xo)/xs)-0.5)*PI*2; //26)/736);
        float py=(((p.y-yo)/ys)-0.5)*PI; //90)/390);
        
        return new RPoint(px,py);
  }
}

class pointTransform  {
  
  RPoint form(RPoint in) {
    //transform
    return new RPoint (((in.x*0.053)+(PI/2)-.01),((in.y)*0.14)+(PI/2)-.1); //front half of sphere
  }
  
}

class sphereMap {
        
    RPoint per(RPoint p,float _r) {
        

        //map to 3D sphere
        float r=getHeight()*_r;
        float x=r*cos(p.x)*(sin(p.y)) *4;
        float z=r*sin(p.x)*(sin(p.y)) -985;
        float y=r*cos(p.y)+(getHeight()*1)-330;  
        
        //camera at 0,0,0 
        //screen plane at 0,0,100
        //(dx-ex)(ez/dz)
        //(dy-ey)(ez/dz)
        
        //rotate camera
        
        //transform into 2D plane @100

        return new RPoint(x*(100/z),y*(100/z));
    }
}

RPoint screenMapper(RPoint p) {
    p.x=((p.x*getWidth())/(PI*2))+(getWidth()/2);
    p.y=((p.y*getHeight())/PI)+(getHeight()/2);
    return p;
}

RPoint plerp(RPoint s,RPoint e,float a) {
  return new RPoint(lerp(s.x,e.x,a),lerp(s.y,e.y,a));
}

float GSphereDist(RPoint p1,RPoint p2) {
  return acos(sin(p1.y)*sin(p2.y)+cos(p1.y)*cos(p2.y)*cos(p1.x-p2.x));
}

RPoint GCircFract(RPoint sp,RPoint ep,float f) {
        float d=GSphereDist(sp,ep);
        float A=sin((1-f)*d)/sin(d);
        float B=sin(f*d)/sin(d);
        float x = A*cos(sp.y)*cos(sp.x) +  B*cos(ep.y)*cos(ep.x);
        float y = A*cos(sp.y)*sin(sp.x) +  B*cos(ep.y)*sin(ep.x);
        float z = A*sin(sp.y)           +  B*sin(ep.y);
        return new RPoint(atan2(y,x),atan2(z,sqrt(pow(x,2)+pow(y,2))));
}
void ellipse3d(RPoint pos,float s,float r){
  //take ellipse coord in UV space
  //return in 3d globe space
  float er=s*0.5;
  curveTightness(-0.67);
  RPoint s1=smap.per(new RPoint(pos.x-er,pos.y),r);
  RPoint s2=smap.per(new RPoint(pos.x,pos.y-er),r);
  RPoint s3=smap.per(new RPoint(pos.x+er,pos.y),r);
  RPoint s4=smap.per(new RPoint(pos.x,pos.y+er),r);
  beginShape();
    curveVertex(s1.x+hw, s1.y+hh);
    curveVertex(s2.x+hw, s2.y+hh);
    curveVertex(s3.x+hw, s3.y+hh);
    curveVertex(s4.x+hw, s4.y+hh);
    curveVertex(s1.x+hw, s1.y+hh);
    curveVertex(s2.x+hw, s2.y+hh);
    curveVertex(s3.x+hw, s3.y+hh);
  endShape();
}

RShape shp;
pointNormalise pnorm;
pointTransform ptrans;
sphereMap smap;

String mode;

csvloader data;
calldata calls;

bitmapcountry Ireland;
bezierstroke bstroke;
gradientstroke3D gstroke;

float hw,hh;

PImage lightmap;
color bg,fg;

void setup(){
  println("vodaviz v0.21");
  RG.init(this);

  mode="PDF";
  if (mode=="PDF") size(832,220,PDF, "vodaviz_circles_190212.pdf"); //P3D); //832,220); //nb pdf is 800x600
  else size(832,220); //,PDF, "testoutput.pdf"); //P3D); //832,220); //nb pdf is 800x600
  
  // C. 33  -  M. 3  -  Y. 0  -  K. 0
  bg=color(164,215,244);
  fg=color(#D2131D);
  
  hw=width/2;
  hh=height/2;
  
  smooth();
  float m = millis();
  
  shp = RG.loadShape("countries_named_mercator.svg"); //test_drawing.svg"); //world_countries_outlines_split.svg");
  
  pnorm = new pointNormalise(18.279,746.302,109,374.293); 
  ptrans = new pointTransform();
  smap = new sphereMap();
  
  lightmap=loadImage("vodaviz_bg_160212.png");
  
  float startsize=0.2;
  float endsize=0.1;
  float linewidth=0.05;
  float mpfract=0.5;
  float raisefract=0.1;
  float bezierfract=0.25;
  bstroke = new bezierstroke(startsize,endsize,linewidth,mpfract,raisefract,bezierfract);
  
  startsize=2;
  endsize=.2;
  float[] transpos={0.0,0.25,0.75,1.0};
  float[] transamt={0.5,0.05,0.05,0.5};
  color col=color(0,0,0);
  gstroke=new gradientstroke3D(startsize,endsize,transpos,transamt,col,smap);
  
  RG.ignoreStyles();
  println("loaded svg in "+((millis()-m)*.001)+" seconds");
  Ireland=new bitmapcountry("Ireland",0,0,0,0,0,0,shp.children[0]);
  Ireland.analyse(18.279,746.302,109,374.293,5);
  data=new csvloader("calls.csv");
  calls=new calldata(data.data,shp,"countries5.xml",5,this);
 
  background(255,255,255);
  noFill();
  stroke(255);
  strokeWeight(.002);
  
  if (false) { //check worked example http://williams.best.vwh.net/avform.htm#Example
    RPoint LAX=new RPoint(2.066470,0.592539);
    RPoint JFK=new RPoint(1.287762,0.709186);
    println("LAX to JFK:"+GSphereDist(LAX,JFK)+" radians");
    RPoint LmJ=GCircFract(LAX,JFK,0.4);
    println("40% of LAX to JFK:"+LmJ.y+","+LmJ.x+" radians");
  }
  

  
}

int i=0;
int j=0;

void draw() {
  if (calls.loaded) { // wait for 1st xmlevent
    noStroke();
    
    float hw=getWidth()/2;
    float hh=getHeight()/2;
    
   if (i==0) {
     //image(lightmap,0,0,getWidth(),getHeight());
     //fill(bg);
     //rect(0,0,width,height);
     //background(0);
     
     if (false) { //render lightmap in 3D
       float pw=0.5/lightmap.width;
       float ph=0.5/lightmap.height;
          for (int j=0;j<lightmap.width;j++) {
            for (int k=0;k<lightmap.height;k++) {
              fill(lightmap.get(j,k));
              float u=(((float)j)/lightmap.width)-0.5;
              float v=(((float)k)/lightmap.height)-0.5;
              RPoint p0=smap.per(ptrans.form(new RPoint((u-pw)*TWO_PI,(v+ph)*PI)),4);
              RPoint p1=smap.per(ptrans.form(new RPoint((u+pw)*TWO_PI,(v+ph)*PI)),4);
              RPoint p2=smap.per(ptrans.form(new RPoint((u-pw)*TWO_PI,(v-ph)*PI)),4);
              RPoint p3=smap.per(ptrans.form(new RPoint((u+pw)*TWO_PI,(v-ph)*PI)),4);
              quad(p0.x+hw,p0.y+hh,p1.x+hw,p1.y+hh,p3.x+hw,p3.y+hh,p2.x+hw,p2.y+hh); //clockwise
              //quad(p0.x,p0.y,p1.x,p1.y,p3.x,p3.y,p2.x,p2.y); //clockwise
            }
          }
      }
   }
   //
    
    //pick a random colour
    fg=color(random(150)+10,random(150)+10,random(150)+10);
    //stroke(fg);
    noStroke();
    fill(fg);
    //stroke(0,0,0,50);
  
    
    if (true) //draw globe 1 country per frame
    {
         if (false) { //draw countries
         //noStroke();
         //fill(fg);
          println("drawing "+calls.getcountry(i).name+": "+calls.getcountry(i).outline.paths.length+" paths");
          
          for (int l=0;l<calls.getcountry(i).outline.paths.length;l++) {
            beginShape();
            for (int k=0;k<calls.getcountry(i).outline.paths[l].commands.length;k++) {
              RPoint sp=calls.getcountry(i).outline.paths[l].commands[k].startPoint;
              RPoint dp=smap.per(ptrans.form(pnorm.alise(sp)),4);
              
              vertex(dp.x+(getWidth()/2),dp.y+(getHeight()/2)); //,z); //z);
            }
            endShape();
            }
          
         }
         
         if (false) { //draw centre lines from countries to ireland
            RPoint Sp=ptrans.form(pnorm.alise(calls.countries.get(i).getcentre()));
            RPoint Ep=ptrans.form(pnorm.alise(Ireland.getcentre()));
            gstroke.drawstroke(Sp,Ep,fg);
         }       
        
        if (false) { //draw great circle lines
          for (int j=0;j<calls.countries.get(i).calls*.001;j++) {
            RPoint s=calls.countries.get(i).getpoint();
            RPoint e=Ireland.getpoint();
             if (s.x>0&&e.x>0){ //point found
              RPoint sp=ptrans.form(pnorm.alise(s));
              RPoint ep=ptrans.form(pnorm.alise(e));
              RPoint Sp=smap.per(sp,4);
              RPoint Ep=smap.per(ep,4);
              RPoint Mp=smap.per(GCircFract(sp,ep,0.8),4.1);
                beginShape();
                  line(Sp.x+(getWidth()/2),Sp.y+(getHeight()/2),Mp.x+(getWidth()/2),Mp.y+(getHeight()/2));
                  line(Mp.x+(getWidth()/2),Mp.y+(getHeight()/2),Ep.x+(getWidth()/2),Ep.y+(getHeight()/2));
                endShape();
            }
          }
            //println("plotting "+calls.countries.get(i).name+": "+calls.countries.get(i).calls+" calls");
          }
         
          if (true) { //draw gradient 3D lines
          if (calls.countries.get(i).points.size()>0) {
            for (int j=0;j<log10(calls.countries.get(i).calls)*15;j++) {
              weightedpixel px=calls.countries.get(i).getpixel();
              RPoint Sp=ptrans.form(pnorm.alise(px.randompt(5.0)));
              RPoint Ep=ptrans.form(pnorm.alise(Ireland.getpoint()));
              //gstroke.drawstroke(Sp,Ep,fg);
              //RPoint Dp=smap.per(Sp,4);
              //ellipse(Dp.x+hw,Dp.y+hh,2,2);
              ellipse3d(Sp,(px.bright*.00001)+.000002,4);
            }
          }
       }   
    }
    else //draw 2D
    {
       if (false) { //draw countries
          beginShape();
          for (int k=0;k<calls.getcountry(i).outline.paths[0].commands.length;k++) {
            RPoint sp=calls.getcountry(i).outline.paths[0].commands[k].startPoint;
            RPoint dp=screenMapper(pnorm.alise(sp));
            
            vertex(dp.x,dp.y); //,z); //z);
          }
          endShape();
          //println("drawing "+calls.getcountry(i).name+": "+calls.getcountry(i).outline.paths[0].commands.length+" points");
        }
       if (true) { //draw lines
          for (int j=0;j<log10(calls.countries.get(i).calls)*1;j++) {
            RPoint s=calls.countries.get(i).getpoint();
            RPoint e=Ireland.getpoint();
             if (s.x>0&&e.x>0){ //point found
              RPoint Sp=screenMapper(pnorm.alise(s));
              RPoint Ep=screenMapper(pnorm.alise(e));
              RPoint Mp=plerp(Sp,Ep,0.75);
                //beginShape();
                  //line(Sp.x+(getWidth()/2),Sp.y+(getHeight()/2),Mp.x+(getWidth()/2),Mp.y+(getHeight()/2));
                  //line(Mp.x+(getWidth()/2),Mp.y+(getHeight()/2),Ep.x+(getWidth()/2),Ep.y+(getHeight()/2));
                  //bezier(Sp.x, Sp.y, Sp.x, Sp.y-((getHeight()-Sp.y)*.1), Ep.x, Ep.y-((getHeight()-Ep.y)*.1), Ep.x, Ep.y);
                //endShape();
                bstroke.drawstroke(Sp,Ep);
             }
          }
            //println("plotting "+calls.countries.get(i).name+": "+calls.countries.get(i).calls+" calls");
          }
          if (false) { //fill country
          for (int j=0;j<calls.countries.get(i).points.size();j++) {
              RPoint Sp=calls.countries.get(i).points.get(j).pt;
              ellipse(Sp.x,Sp.y,1,1);
          }
        }
    }
  
         
    i++;
  
    //}
    if (i==calls.countries.size()-1) {
          println("finished");
         noLoop();
         if (mode=="PDF") exit();
         else save("160212_test.png");
    }
  }
}

void mousePressed() {
  
}