Normalizing the scale of zoomable signal plots?

[__NullPointer __exception]

I've taken to adapt the realtime_signals demo to my datasets. The data I work with is non-stationary and usually exhibits a trend of sorts. As a consequence, zooming into the data usually doesn't work well since drifting values fall out of the window.

I'm trying to add min-max normalization along the y-axis to the example code. That way, whatever my zoom level is, data should remain centered on the window. However, my solution produces glitches I can't explain.

from vispy import gloo
from vispy import app
import numpy as np
import math
import numpy.ma as ma



#Data
num_samples = 1000
num_features = 3
df_raw = np.reshape(1+(np.random.normal(size=num_samples*num_features, loc=[0], scale=[0.01])), [num_samples, num_features]).cumprod(0).astype('float32')
df_raw_std = 1+((df_raw-np.min(df_raw,0))*2)/(np.min(df_raw,0)-np.max(df_raw,0))

# Generate the signals as a (num_features, num_samples) array. 320x1000
y = df_raw_std.transpose()

# Signal 2D index of each vertex (row and col) and x-index (sample index
# within each signal).
index_col = np.c_[np.tile(np.arange(num_samples), num_features),np.repeat(np.arange(num_features), num_samples)].astype(np.float32)
y_flat = np.reshape(y,y.shape[0]*y.shape[1])
index_y_scaled_orig = np.c_[-1 + 2*(index_col[:,0] / num_samples),y_flat].astype(np.float32)
index_y_scaled = index_y_scaled_orig.copy()

index_min = np.c_[np.tile(np.arange(num_samples), num_features),np.repeat(1, num_samples*num_features)].astype(np.float32)
index_max = np.c_[np.tile(np.arange(num_samples), num_features),np.repeat(1, num_samples*num_features)].astype(np.float32)


#This is called once for each vertex
VERT_SHADER = """
#version 120
//scaling. Running minimum and maximum of visible time series
attribute float index_min;
attribute float index_max;

// y coordinate of the position.
attribute float y;

// row, and time index.
attribute vec2 index_col;

// 2D scaling factor (zooming).
uniform vec2 scale;
uniform vec2 num_features;

// Number of samples per signal.
uniform float num_samples;

// for fragment shader
varying vec2 v_index;

// Varying variables used for clipping in the fragment shader.
varying vec2 v_position;
varying vec4 v_ab;

void main() {
    float nrows = num_features.x;

    // Compute the x coordinate from the time index
    float x = -1 + 2 * (index_col.x / (num_samples-1));

    //0 is zoom from center. 1 is zoom on the right. -1 is zoom on the left. WE should map mouse x pos to this.
    float zoom_x_pos = 0.0;

    // RELATIVE LINE POSITION
    // =============================
    // Manipulate x/y position here?
    // =============================
    // vec2 position = vec2(x - (1 - 1 / scale.x)*zoom_x_pos,y); // DEACTIVATED SCALING, NICE PLOTS EMERGE
    vec2 position = vec2(x - (1 - 1 / scale.x)*zoom_x_pos,(y-index_min)/(index_max-index_min)); //SCALING, GLITCHY

    // SPREAD
    //does not scale the x pos, just the y pos by an equal amount per row
    float spread = 1;
    vec2 a = vec2(spread, spread/nrows);

    // LOCATION
    vec2 b = vec2(0, -1 + 2*(index_col.y+.5) / nrows);

    // COMBINE RELATIVE LINE POSITION + SPREAD + LOCATION
    gl_Position = vec4(a*scale*position+b, 0.0, 1.0);

    // WRAP UP
    v_index = index_col;
    // For clipping test in the fragment shader.
    v_position = gl_Position.xy;
    v_ab = vec4(a, b);
}
"""

FRAG_SHADER = """
#version 120
varying vec2 v_index;
varying vec2 v_position;
varying vec4 v_ab;
void main() {
    gl_FragColor = vec4(1., 1., 1., 1.);
    // Discard the fragments between the signals (emulate glMultiDrawArrays).

    if (fract(v_index.y) > 0.)
        discard;

    // Clipping test.
    vec2 test = abs((v_position.xy-v_ab.zw)/v_ab.xy);
    if ((test.x > 1) || (test.y > 1))
        discard;
}
"""


class Canvas(app.Canvas):
    def __init__(self):
        app.Canvas.__init__(self, title='Use your wheel to zoom!',
                            keys='interactive')
        self.program = gloo.Program(VERT_SHADER, FRAG_SHADER)
        self.program['y'] = y.reshape(-1, 1)
        self.program['index_col'] = index_col
        self.program['scale'] = (1., 1.)
        self.program['num_features'] = (num_features, 1)
        self.program['num_samples'] = num_samples
        self.program['index_min'] = index_min[:,0].reshape(-1, 1)
        self.program['index_max'] = index_max[:,0].reshape(-1, 1)
        gloo.set_viewport(0, 0, *self.physical_size)

        self._timer = app.Timer('auto', connect=self.on_timer, start=True)

        gloo.set_state(clear_color='black', blend=True,
                       blend_func=('src_alpha', 'one_minus_src_alpha'))

        self.show()

    def on_resize(self, event):
        gloo.set_viewport(0, 0, *event.physical_size)

    def on_mouse_wheel(self, event):
        dx = np.sign(event.delta[1]) * .05
        scale_x, scale_y = self.program['scale']

        index_y_scaled[:,0] = index_y_scaled_orig[:,0] * scale_x
        index_y_scaled[:, 1] = index_y_scaled_orig[:, 1] * scale_x
        valid = ((index_y_scaled[:,0]>-1)*(index_y_scaled[:,0]<1))

        index_y_scaled_reshaped = (np.reshape(index_y_scaled[:, 1],[num_features,num_samples]))
        shown = ma.masked_array(index_y_scaled_reshaped, mask=np.logical_not(valid))
        runmin = np.array(np.min(shown, 1))
        runmax = np.array(np.max(shown, 1))
        index_min[:, 1] = np.repeat(runmin, num_samples)
        index_max[:, 1] = np.repeat(runmax, num_samples)

        print(scale_x)
        print(runmin)
        print(runmax)

        self.program['index_min'] = index_min[:,1].reshape(-1, 1)
        self.program['index_max'] = index_max[:,1].reshape(-1, 1)
        #print(self.program['print_position'])

        scale_x_new, scale_y_new = (scale_x * math.exp(1.0*dx),
                                    scale_y * math.exp(1.0*dx))
        #print(scale_x_new)
        self.program['scale'] = (max(1, scale_x_new), max(1, scale_y_new))
        self.update()

    def on_timer(self, event):
        """Add some data at the end of each signal (real-time signals)."""
        self.program['y'].set_data(y.ravel().astype(np.float32)) #(10920,)
        self.update()

    def on_draw(self, event):
        gloo.clear()
        self.program.draw('line_strip')

if __name__ == '__main__':
    c = Canvas()
    app.run()

1. What am I doing wrong? i am "estimating" the visible line selection outside opengl and passing scale correction parameters to the opengl pipeline. Yet I get obvious visual glitches as well as distorted lines
2. Is there a smarter way of approaching this problem in in vispy? Perhaps a way to solve the normalization in the fragment shader or via camera tricks?

[David Hoese]

Hi,

This is a difficult one. There is a lot of code here. Are you saying that if you swap the two lines in your vertex shader things look mostly as you expect (at least that's what I got from your code comments)?

My initial guess is that there is something odd going on with the limits. If I hard code index_min to -2 and index_max to 2 I get a slightly better plot, but hard to say.

Dave

[__NullPointer __exception]

Hello David

Apologies for the delayed response.

turns out you were indeed correct about the limits, the glitches were produced by a bug in the min max computation. I've posted the correct code below.

Still, I am not happy with my current solution. It requires the normalization parameters to still be computed outside opengl, potentially massively slowing down execution once I move to higher data amounts and many lineplots to visualize.

I understand that it is not possible to do this in the vertex shader since it only operates over individual vertices - normalization requires knowledge of the relative position of all other vertices in a given lineplot.

I wonder however if it is possible to achieve normalization - which is just a linear affine scale and transform of a line object - for example within the geometry shader. I am still new to opengl, but if I understand the pipeline correctly it will require each lineplot (3 in my example code) to be defined in its own vertex buffer primitive. I can then consume these primitives in a geometry shader, iterate over the vertices of a given lineplot using a loop over gl_in, compute the overall min and max and then translate and scale the position of each vertex in gl_in.

Is this possible with vispy? I know I can define a geometry shader, but I have trouble actually creating individual primitives for each line plot. In my example code, i think I am treating all 3 line plots as a single primitive.

Greetings,

Esteban

_____________

Corrected code

______________

from vispy import gloo
from vispy import app
import numpy as np
import math
import numpy.ma as ma

import matplotlib.pyplot as plt
import pandas as pd


#Data
num_samples = 10000

num_features = 3
df_raw = np.reshape(1+(np.random.normal(size=num_samples*num_features, loc=[0], scale=[0.01])), [num_samples, num_features]).cumprod(0).astype('float32')
df_raw_std = 1+((df_raw-np.min(df_raw,0))*2)/(np.min(df_raw,0)-np.max(df_raw,0))

# Generate the signals as a (num_features, num_samples) array. 320x1000
y = df_raw_std.transpose()

# Signal 2D index of each vertex (row and col) and x-index (sample index
# within each signal).
index_col = np.c_[np.tile(np.arange(num_samples), num_features),np.repeat(np.arange(num_features), num_samples)].astype(np.float32)

y_flat = y.flatten()

"""
#version 120
//scaling. Running minimum and maximum of visible time series
attribute float index_min;
attribute float index_max;

// y coordinate of the position.
attribute float y;

// row, and time index.
attribute vec2 index_col;

// 2D scaling factor (zooming).
uniform vec2 scale;
uniform vec2 num_features;

// Number of samples per signal.
uniform float num_samples;

// for fragment shader
varying vec2 v_index;

// Varying variables used for clipping in the fragment shader.
varying vec2 v_position;
varying vec4 v_ab;

void main() {
    float nrows = num_features.x;

    // Compute the x coordinate from the time index
    float x = -1 + 2 * (index_col.x / (num_samples-1));

    //0 is zoom from center. 1 is zoom on the right. -1 is zoom on the left. WE should map mouse x pos to this.
    float zoom_x_pos = 0.0;

    // RELATIVE LINE POSITION
    // =============================
    // Manipulate x/y position here?
    // =============================
    // vec2 position = vec2(x - (1 - 1 / scale.x)*zoom_x_pos,y); // DEACTIVATED SCALING, NICE PLOTS EMERGE
    vec2 position = vec2(x - (1 - 1 / scale.x)*zoom_x_pos,y); //SCALING, GLITCHY

    vec2 yscale_a = vec2(0., index_min);
    vec2 yscale_b = vec2(1., 2/(index_max-index_min));
    vec2 yscale_c = vec2(0., -1/nrows);


    // SPREAD
    //does not scale the x pos, just the y pos by an equal amount per row
    float spread = 1;
    vec2 a = vec2(spread, spread/nrows);

    // LOCATION
    vec2 b = vec2(0, -1 + 2*(index_col.y+.5) / nrows);

    // COMBINE RELATIVE LINE POSITION + SPREAD + LOCATION
    // gl_Position = vec4(a*(scale*position-yscale_a)*yscale_b+b, 0.0, 1.0);
    gl_Position = vec4(a*(scale*position-yscale_a)*yscale_b+b+yscale_c, 0.0, 1.0);



    // WRAP UP
    v_index = index_col;
    // For clipping test in the fragment shader.
    v_position = gl_Position.xy;
    v_ab = vec4(a, b);
}
"""

FRAG_SHADER = """
#version 120
varying vec2 v_index;
varying vec2 v_position;
varying vec4 v_ab;
void main() {
    gl_FragColor = vec4(1., 1., 1., 1.);
    // Discard the fragments between the signals (emulate glMultiDrawArrays).

    if (fract(v_index.y) > 0.)
        discard;

    // Clipping test.
    vec2 test = abs((v_position.xy-v_ab.zw)/v_ab.xy);
    if ((test.x > 1) || (test.y > 1))
        discard;
}
"""


class Canvas(app.Canvas):
    def __init__(self):
        app.Canvas.__init__(self, title='Use your wheel to zoom!',
                            keys='interactive')
        self.program = gloo.Program(VERT_SHADER, FRAG_SHADER)

        self.program['y'] = y_flat
        

self.program['index_col'] = index_col
        self.program['scale'] = (1., 1.)
        self.program['num_features'] = (num_features, 1)
        self.program['num_samples'] =

 num_samples
        self.program['index_min'] = index_min[:,0].flatten()
        self.program['index_max'] = index_max[:,0].flatten()

        gloo.set_viewport(0, 0, *self.physical_size)

        self._timer = app.Timer('auto', connect=self.on_timer, start=True)

        gloo.set_state(clear_color='black', blend=True,
                       blend_func=('src_alpha', 'one_minus_src_alpha'))

        self.show()

    def on_resize(self, event):
        gloo.set_viewport(0, 0, *event.physical_size)

    def on_mouse_wheel(self, event):
        dx = np.sign(event.delta[1]) * .05
        scale_x, scale_y = self.program['scale']

        index_y_scaled[:,0] = index_y_scaled_orig[:,0] *

 scale_x
        valid = ((index_y_scaled[:,0]>-1)*(index_y_scaled[:,0]<1))
        y_flat_scaled = y_flat * scale_x
        shown = ma.masked_array(y_flat_scaled, mask=np.logical_not(valid))
        shown_reshaped = (shown.reshape(num_features,num_samples))

        runmin = np.array(np.min(shown_reshaped, 1))
        runmax = np.array(np.max(shown_reshaped, 1))

        index_min[:, 1] = np.repeat(runmin, num_samples)
        index_max[:, 1] = np.repeat(runmax, num_samples)

        # scale_x = 10
        # scale_y = 10
        # print(scale_x)
        # print(scale_y)
        # print(runmin)
        # print(runmax)
        # forplot=(y_flat_scaled*valid).reshape(num_features,num_samples).transpose()
        # pd.DataFrame(forplot).plot(subplots=True)
        # forplot2 = (((y_flat_scaled * valid - index_min[:, 1])).flatten()).reshape([num_features, num_samples]).transpose()
        # pd.DataFrame(forplot2).plot(subplots=True)
        # forplot3 = (((y_flat_scaled * valid - index_min[:, 1]) / (index_max[:, 1] - index_min[:, 1])).flatten()).reshape([num_features, num_samples]).transpose()
        # pd.DataFrame(forplot3).plot(subplots=True)

        self.program['index_min'] = index_min[:,1].flatten()
        self.program['index_max'] = index_max[:,1].flatten()

        #print(self.program['print_position'])

        scale_x_new, scale_y_new = (scale_x * math.exp(1.0*dx),
                                    scale_y * math.exp(1.0*dx))
        #print(scale_x_new)
        self.program['scale'] = (max(1, scale_x_new), max(1, scale_y_new))
        self.update()

    def on_timer(self, event):
        """Add some data at the end of each signal (real-time signals)."""

        # y[:, :-1] = y[:, 1:]
        # y[:, -1:] = np.random.normal(size=3, loc=[0], scale=[0.01]).reshape(3, 1)

        self.program['y'].set_data(y.flatten().astype(np.float32)) #(10920,)

[__NullPointer __exception]

i believe I've reached the line limit - allow me to re-post the corrected code:

[David Hoese]

I've never thought of using geometry shaders that way, but I also haven't coded my own geometry shader before. It is technically possible to use geometry shaders in newer versions of vispy but I've had mixed results getting them to work on all platforms (especially OSX).

How much data (numpy array shape of each line) do you plan on working with? If/when data is updated, are new points added to the end of the line and old points removed from the other end? Or do you replace the entire line with new data for every update? Sorry if you've answered this before.

In my opinion using a geometry shader or trying to do min/max calculations in the GPU is overkill. I'd be worried that you are prematurely optimizing your code. There are also similar options that might perform well for you and produce a similar end result visualization-wise. For example, you could hold on to a min/max value for each line in the python CPU-side code, then when new data comes in you compare that to the min/max of the new data (assuming the entire line isn't replaced). Or only computing the min/max every X updates and provided a larger Y-axis buffer around the data. So instead of doing min - 0.1 to max + 0.1 (I'm thinking of this like a matplotlib line plot), you could do a percentage of the difference like min + (max - min) * 0.4 and max + (max - min) * 0.4. This way you aren't computing the min/max every update and the data probably fits on the screen for quite a while. There are downsides to all of this, but it also saves you from having to write a geometry shader.

Just brainstorming.

Dave

[__NullPointer __exception]

Hi David,

Indeed, the displayed data array is a stream and updated in FIFO style, oldest row to have entered is the first row to be removed with incoming data. With this in mind, good point regarding the running min / max calculation, as the running comparison is indeed computationally negligible then. 

The problem or downside comes more with the concept of user interaction and zooming. The visualization is used in a GUI and the data load is heavy and may well succeed a million points per column at any given moment. Because of the granularity of the data, zooming is essential. The user should be able to zoom in to view tiny details of the time series and subsequently select these. As detailed in my first post, normalization has to be tied to the zoom factor in order to be able to see anything at close zoom levels and also do a selection. Hence, every time a user interacts with zoom, minmax has to be recalculated accross the entire +visible" portion of the number array of potentially millions datapoints per column.

I agree with you that this may be a bit premature in terms of optimization though - still, I think of it as a nice challenge to get to learn opengl and use the shader pipelines. If not in the geometry shader, am I correct in that already the fragment shader can operate over an entire priitive? By re-defining my code to have one line-strip per column of data, I might be able to try something in the fragment shader.

My trouble is at this stage in how to re-work the code so that it builds one line-strip per column of data. 

[David Hoese]

I think I understand what you're saying. As the user zooms in on a particular point in the X-dimension, the line should be zoomed in to a particular range in the Y dimension? So let's say at time T, your line plot from T-1 to T+1 spans from 50 to 75 (for example). If the user's mouse is at time T, then as they zoom in with the scroll wheel the line will be scaled to show the data between 50 and 75. As they continue zooming in we aren't looking at T-1 to T+1 anymore, but rather T-0.5 to T+0.5 so our y limits are now 60 and 68 (for example). I'm not sure that description makes sense to you, but I think it matches what you want and makes sense to me (for now) :D. This is different from how most plotting libraries work, isn't it? Usually the zoom follows the user's mouse in both the X and Y dimensions, right?

The fragment shader operates on...fragments. So even less "big picture" than the vertex shader. I'm not sure you'll be able to determine what you want there.

Dave

[__NullPointer __exception]

Exactly :) The last code I posted simulates the behavior correctly, albeit outside OpenGl.

It is indeed different than what regular plotting usually does. The reason for this approach being that the user must then "range-select" a granular piece of line from a very large dataset for further downstream computation. Because the data has lots of structural breaks, there will be locations with extreme growth rates which will become tedious to select if the zooming also affects the y-axis. I want the user to always have a complete view of the slice he is selecting while preserving the granularity of selection flexibility on the X axis (user may want to select a slice 3 points wide in 1 million points).

I guess then it really comes down to perhaps giving geometry shaders a go. I understand the concerns with stability, but I consider it a fun little project for during my downtimes..! Sorry to ask so upfront, but do you perhaps have an example of drawing a line-strip primitive per column, given an input 2D-numpy matrix? Most examples I find for vispy and line drawing use the high-level interface which I'd like to avoid in order to have full flexibility later on.

 

[David Hoese]

Regarding user behavior: sounds good. Good luck.

For lines, I would start by looking at the LineVisual docstring if you haven't already (https://github.com/vispy/vispy/blob/master/vispy/visuals/line/line.py#L54-L91). However, because of your need for custom shaders and stuff this might not be great. The basic GL line drawing method provides the most flexibility for how you provide the data to the GPU and might give you ideas for using index buffers or something similar to draw each line segment separately. Otherwise, no I don't think I have a good example. There might be some examples in the vispy repository that show the different uses of the LineVisual object which might also give you some ideas of what the docstring is trying to say.

Dave