#!/usr/bin/env python
# coding: utf-8

# In[ ]:


import numpy as np 
import time
from IPython.display import clear_output

from six.moves.urllib.request import urlopen
from contextlib import closing
import json
import k3d

plot = k3d.plot()
plot.display()


# In[ ]:


from k3d.helpers import download
from pyunpack import Archive

filename = download('http://www.semantic3d.net/data/point-clouds/testing1/stgallencathedral_station1_intensity_rgb.7z')
Archive(filename).extractall('./')


# In[ ]:


np.fromfile(filename.replace('.7z', '.txt'), sep=' ', dtype=np.float32).shape


# In[ ]:


import csv

data = None

with open(filename.replace('.7z', '.txt'), mode='r') as csv_file:
    csv_reader = csv.reader(csv_file, delimiter=' ')    
    data = np.array(list(csv_reader), dtype=np.float32)

# compute color in hex format
data[:, 4] = np.sum(data[:, 4:7].astype(np.uint32) * np.array([1, 256, 256 ** 2]), axis=1)    
data = data[:, 0:5]


# In[ ]:


data.shape


# In[ ]:


plot += k3d.points(data[::2, 0:3], data[::2, 4].astype(np.uint32), point_size=0.05, shader="flat") 


# In[ ]:


plot.camera = [5.251483149143791,
 -7.92683507646606,
 3.144285796928443,
 -2.470283607444292,
 3.6558150584160503,
 2.3721091212696286,
 0,
 0,
 1]


# In[ ]:


plot.camera_auto_fit = False
plot -= plot.objects[0]
plot += k3d.points(data[::50, 0:3], data[::50, 4].astype(np.uint32), point_size=0.25, shader="flat") 


# In[ ]:


sigma=10.0
beta=8./3
rho=28.0
def lorenz_deriv(X, sigma=sigma, beta=beta, rho=rho):
    """Compute the time-derivative of a Lorenz system."""
    x, y, z = X.T
    return np.vstack([sigma * (y - x), x * (rho - z) - y, x * y - beta * z]).T


# In[ ]:


plot.camera = [82.36534387751811,
 -119.8210969123126,
 43.968748841328704,
 -0.7272701043451865,
 4.817824060482123,
 35.65948744314234,
 0,
 0,
 1]


# In[ ]:


for p in plot.objects:
    X = p.positions
    for i in range(150):
        X = X + lorenz_deriv(X, sigma=sigma, beta=beta, rho=rho)*0.002
        if i%15==0 and i>0:
            p.positions = X[::1,:]
        #time.sleep(0.1)
        clear_output(wait=True)
        print("iteration:",i)
    p.positions = X


# In[ ]:


for i in range(15):
    for p in plot.objects:
        X = p.positions
        for j in range(15):
            X = X + lorenz_deriv(X, sigma=sigma, beta=beta, rho=rho)*0.001
        p.positions = X[:,:]
        clear_output(wait=True)
        print("iteration:",i)
    


# In[ ]:


p.point_size = .65


# In[ ]: