lvpyio manual

This package is a reader and writer for DaVis image sets, vector sets and particle fields (.set), as well as vector fields (.vc7, .vec) and images (.im7, .imx).

Supported Python versions are 3.7 to 3.11, CPython only, 64-bit only. OS: Windows and Linux.

Table of Contents

Difference between Lvpyio and Lvpyio_wrapped

There are two packages to choose from: lvpyio and lvpyio_wrapped. Their documented interface is identical. Most users should use lvpyio: it's faster, consumes less memory, and is easier to call from other applications, like MATLAB and ParaView. However, using lvpyio and Qt5 in the same process, e.g. when you embed lvpyio inside a PySide2 application or when you select Qt5Agg as the matplotlib backend, can cause a Qt-related ImportError. If you encounter this error, install and use lvpyio_wrapped as a drop-in replacement. It starts a separate process for reading and writing, so you will see a second Python process appear in the task manager after importing it.

Dependencies

Lvpyio

lvpyio has the following necessary dependencies:

Optional dependencies:

On very minimal Linux installations (e.g. Arch base install), you may have to install libdb as well. On most Linux systems it's already present.

Lvpyio_wrapped

lvpyio_wrapped has the following necessary dependencies:

Optional dependencies:

Installation

Create and activate a virtual environment with a tool of your choice, then install the package with pip:

pip install lvpyio

When in doubt, see the pip documentation.

Notes for Conda Users

If you use conda, please install (or update, if already installed) the necessary dependencies with conda before installing lvpyio with pip. Otherwise the missing dependencies will be installed with pip as well and this is known to cause conflicts. If unsure, do it in a fresh environment or call pip install --no-deps instead of pip install.

Notes for Users of Older Pip Versions

If your pip version is older than 19.3 (you can check it with pip --version), you need to update pip first.

If you use an environment created with venv:

python -m pip install -U pip

If you use a conda environment:

conda update pip

Usage

Reading Sets and Multi-sets

Although you can work with .im7 and .vc7 files directly (see below), the preferred way to read and write DaVis data is via sets. DaVis recordings often don't consist of separate .im7 files, so using the set interface becomes absolutely necessary.

The function read_set accepts a path and returns a Set object or, if the path points to a multi-set, a Python list of Set objects. If you want to check in advance if your set is a multi-set, you can do it with is_multiset.

>>> import lvpyio as lv
>>> path = "path/to/data/Cam_Date=041014_Time=173435"
>>> lv.is_multiset(path)
False
>>> s = lv.read_set(path)
>>> s
<Image set 'Cam_Date=041014_Time=173435' with 10 buffers>
>>> s.title
'Cam_Date=041014_Time=173435'
>>> s.type_id
'SET_TYPE_ID_RECORDING'
>>> len(s)  # number of buffers
10

A set consists of one or more buffers, and a buffer consists of one or more frames (also called volumes).

Buffers and Frames

A buffer can be read from a Set object by indexing it:

# reading the first buffer from a set
buffer = s[0]

# iterating over all buffers
for buffer in s:
    ... # do something

If you read a set, a buffer isn't loaded from disk until you access it, and will be reloaded every time you access it:

import lvpyio as lv
s = lv.read_set("path/to/data") # does not load any buffer
copy = s[0] # loads a buffer from disk
another_copy = s[0] # loads the same buffer again

In particular, modifying the loaded buffer will never modify the set it came from. To force loading all buffers and keep their copies in memory, which isn't recommended for large sets, you can use all_buffers = list(s). For only loading a range of buffers, comprehensions can be used, e.g. [s[i] for i in range(2, 6)] copies the specified buffers from disk to memory immediately, and (s[i] for i in range(2, 6)) creates a lazy iterator without copying anything yet.

Buffer attributes can be accessed with buffer.attributes. It's a Python dict whose keys are str:

>>> from pprint import pprint
>>> pprint(buffer.attributes)
{'AttributeDisplayFrameInfo': 'SetupDialogBufferAttrFrameInfo(-1)',
 'ExpPathValue_Date': '041014',
 'ExpPathValue_Time': '173435',
 'FastStorageStartTime': '22463.467724',
 'LoadFile': 'E:\\examples\\Example PIV '
             'Stereo\\Cam_Date=041014_Time=173435\\B00001.im7',
 'LoadSet': 'E:/examples/Example PIV Stereo/Cam_Date=041014_Time=173435',
 'LoadSetIndex': '1',
 '_Date': '14.10.04',
 '_FrameScale': '1\n0\nframe\n',
 '_Header_PackType': '1',
 '_Time': '17:34:36.906'}

Device data won't be displayed separately like in DaVis. You will find it in the attributes whose name start with DevData.

Frames can be accessed with buffer.frames. It's a Python list of ImageFrame or VectorFrame objects. A Buffer object is sized, subscriptable and iterable:

>>> len(buffer) # number of frames
4
>>> [frame.scales.x.unit for frame in buffer]
['mm', 'mm', 'mm', 'mm']
>>> frame = buffer[0]
>>> frame.type_id
'ImageFrame'

Frame attributes can be accessed just like buffer attributes:

>>> pprint(frame.attributes)
{'AcqTime': '500',
 'AcqTimeSeries': '1530749 µs',
 'CCDExposureTime': '30000 µs',
 'CamPixelSize': '6.45 µm',
 'CameraName': '1A: Imager Intense (1376 x 1040 pixel)',
 'FrameProcessing': '0'}

For example, to investigate the acquisition time series of a set, you can do:

>>> import lvpyio as lv
>>> s = lv.read_set("path/to/data/Cam_Date=041014_Time=173435")
>>> [b[0].attributes["AcqTimeSeries"] for b in s]
['1530749 µs', '1948585 µs', '2366377 µs', '2782542 µs', '3200737 µs',
'3617083 µs', '4035520 µs', '4452483 µs', '4868068 µs', '5286472 µs']

Easy Data Access

>>> arr = buffer.as_masked_array()
>>> type(arr)
<class 'numpy.ma.core.MaskedArray'>

A frame, conceptually, consists of planes. A plane of an image frame is an image, a plane of a vector frame is a vector field. Three-dimensional data has multiple planes, two-dimensional data only one plane. The number of planes can be obtained with len(frame), the shape of each plane with frame.shape. For an image, shape is measured in pixel, for a vector field, in vectors.

>>> len(frame)
1 # one image
>>> frame.shape
(1024, 1376) # height x width, numpy convention

Every frame has a method called as_masked_array which accepts a plane index and returns a masked numpy array. If no plane index is provided, 0 will be assumed. Calling it with different indices is useful when working with 3D data. The values will be scaled like in DaVis. To avoid scaling errors, the dtype of the result will be 'float64' for image frames regardless of the data type stored in the original buffer, [('u', '<f8'), ('v', '<f8')] for vector frames with two components and [('u', '<f8'), ('v', '<f8'), ('w', '<f8')] for vector frames with three components, where f8 again stands for float64.

buffer.as_masked_array(i, j) is shorthand for buffer[i].as_masked_array(j). If no frame and plane indices are provided, 0 will be assumed.

You may modify the resulting arrays as you wish. The modifications won't affect the original buffer object or frame object.

Working with Structured and Masked Numpy Arrays: Examples

The following two examples are runnable. We encourage you to open a Python interpreter and try them out with the accompanying data.

Image buffer example:

import lvpyio as lv
b = lv.read_set("data/word_masked")[0]
# extract first frame, first plane as array:
arr = b.as_masked_array()

print(arr)
# output:
# [[554.0 364.0 2262.0 ... -- -- --]
#  [2453.0 423.0 428.0 ... -- -- --]
#  [1689.0 508.0 382.0 ... -- -- --]
#  ...
#  [-- -- -- ... -- -- --]
#  [-- -- -- ... -- -- --]
#  [-- -- -- ... -- -- --]]

# get unmasked data:
print(arr.data)
# output:
# [[ 554.  364. 2262. ...   67.   77.   76.]
#  [2453.  423.  428. ...   93.   76.   80.]
#  [1689.  508.  382. ...   79.   90.   75.]
#  ...
#  [  65.   59.   69. ...   75.   78.   80.]
#  [  63.   66.   59. ...   71.   84.   83.]
#  [  70.   66.   64. ...   68.   79.   85.]]

# get mask:
print(arr.mask)
# output:
# [[False False False ...  True  True  True]
#  [False False False ...  True  True  True]
#  [False False False ...  True  True  True]
#  ...
#  [ True  True  True ...  True  True  True]
#  [ True  True  True ...  True  True  True]
#  [ True  True  True ...  True  True  True]]

# first row of the masked array:
print(arr[0])
# output:
# [554.0 364.0 2262.0 ... -- -- --]

Vector buffer example:

import lvpyio as lv
b = lv.read_set("data/2d3c")[0]
arr = b.as_masked_array()

# get the topmost leftmost vector (it's masked in this example):
print(arr[0,0])
# output:
# (--, --, --)

# get the first five vectors from the second row:
print(arr[1,:5])
# output:
# [(--, --, --) (--, --, --) (--, --, --)
# (-0.16445358684289454, 0.15648226495242118, 0.0090000756944865)
# (-0.17152579677414895, 0.15407613824236394, -0.005028278277654201)]

# extract the unmasked data and the mask:
arr.data
arr.mask
# extract the vector components as masked arrays:
arr["u"]
arr["v"]
arr["w"]  # this component only exists in 3C vector fields

More details about structured arrays and masked arrays can be found in the numpy documentation:

Structured arrays

Masked arrays

Plotting (Requires matplotlib)

Every frame has a plot method:

>>> frame.plot()  # a plot should appear

An image will be displayed with masked pixels whited out.

A vector field will be displayed as arrows over a background map that shows the absolute velocity (or, more generally, the vector length). If vectors have a third component, it's being ignored by the arrows but included in the background image calculation. Masked and disabled vectors will be whited out.

Just like as_masked_array, plot accepts a plane index. If no plane index is provided, 0 will be assumed. buffer.plot(i, j) is shorthand for buffer[i].plot(j). If no frame and plane indices are provided, 0 will be assumed.

A keyword argument show (default is True) can be used if you are using a non-interactive backend or want to customize the plot. With show=False, the method doesnt display a plot and returns a matplotlib.figure.Figure object instead which you can modify and display as you need.

Advanced Data Access

The easy interface is limited. It also assumes well-formed data (the same shape for every plane, mask having the same shape as the data itself, no vector frames in an image set etc.).

A frame has a components attribute which behaves like a dict and can be used for advanced data access.

>>> vector_frame.components
<17 components: U0, V0, U1, V1, U2, V2, U3, V3, ACTIVE_CHOICE, ENABLED,
TS:Peak ratio, MASK, TS:Correlation value, TS:Particle size, TS:Uncertainty V,
TS:Uncertainty Vx, TS:Uncertainty Vy>

Its console representation already contains the names of the components. This is useful to interactively explore attached scalar fields, vector components, vector choices etc. Alternatively, you can obtain the names of all components with frame.components.keys().

Scalar fields' names are prefixed with "TS:". For example, if DaVis displays a field called Particle size, it can be obtained with vector_frame.components["TS:Particle size"].

For example, in an image frame, you can get the raw unscaled image data of the first plane with frame.components["PIXEL"] to do the appropriate scaling yourself. The scales can be accessed with frame.scales. If k is the plane index, j the row index and i the column index, you can get the coordinates of the corresponding point like this:

x = frame.scales.x.slope * i + frame.scales.x.offset
y = frame.scales.y.slope * j + frame.scales.y.offset
z = frame.scales.z.slope * k + frame.scales.z.offset

If frame is a vector frame, the grid should be taken into account while scaling coordinates. Use frame.scales.x.slope * grid.x * i + frame.scales.x.offset etc.

The data itself can be scaled with frame.scales.i in the same manner (i means intensity, just like in DaVis). Unlike with coordinates, taking the grid into account is not needed. When scaling vector data, both DaVis and as_masked_array a) scale every vector component with the intensity scale, and then b) multiply the u component with the sign of the slope of the x scale, the v component with the sign of the slope of the y scale, and the w component with the sign of the slope of the z scale. You should do the same if you wish to obtain the same results manually. Step b usually only makes a difference for the y component: in most DaVis data, x and z scales have positive slope and the y scale has negative slope.

Scalar fields, like all components, have scales as well:

vector_frame.components["TS:Particle size"].scale

In a vector frame, as_masked_array doesn't differentiate between disabled vectors and masked vectors. You can access the raw data with frame.components["MASK"] and frame.components["ENABLED"]. Be aware that those components aren't always present. Be also aware that DaVis and numpy use different conventions for masks. Both use an array of ones and zeros as a mask, but in DaVis, 0 means invalid value and 1 means valid value, and in numpy it's the other way around. frame.masks and frame.disabled return the corresponding components as they are stored (0 means invalid value), but frame.as_masked_array uses numpy masked arrays and, therefore, converts masks to numpy format.

To sum up the differences between easy access and advanced access:

# Assuming we have a vector buffer:
array = buffer.as_masked_array(frame=m, plane=n)
# Accessing the mask (no distinction between disabled and masked,
# 1 means invalid as usual in numpy):
array["u"].mask
# Accessing the u component of the data (scaled, choices taken into account):
array["u"].data

# U component of the data, choice 0, unscaled:
buffer[m].components['U0'].planes[n]
# Mask (1 means valid as usual in DaVis):
buffer[m].masks[n]
# Which vectors are enabled:
buffer[m].enabled[n]

The advanced interface does not copy the original data, so it can be used to modify it.

Example: Reading the Uncertainties

>>> s = lv.read_set("data/2d3c")
>>> buffer = s[0]
>>> frame = buffer[0]
>>> frame.components
<25 components: U0, V0, U1, V1, U2, V2, U3, V3, ACTIVE_CHOICE, ENABLED, W0, W1,
W2, W3, TS:Peak ratio, MASK, TS:Correlation value, TS:Stereo reconstruction 
error, TS:Uncertainty C_xy, TS:Uncertainty C_xz, TS:Uncertainty C_yz, 
TS:Uncertainty V, TS:Uncertainty Vx, TS:Uncertainty Vy, TS:Uncertainty Vz>
>>> uncertainty_V = frame.components["TS:Uncertainty V"]
>>> uncertainty_V
<1 plane, scale=Scale(slope=0.213246, offset=0.0, unit='m/s', description=
'Uncertainty V')>
>>> uncertainty_V[0] # get the plane with the actual data
array([[0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       [0.        , 0.        , 0.05930465, ..., 0.        , 0.        ,
        0.        ],
       ...,
       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ]], dtype=float32)
# we can ignore the offset because it's 0.0:
>>> uncertainty_V_scaled = uncertainty_V[0] * uncertainty_V.scale.slope
>>> uncertainty_V_scaled
array([[0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       [0.        , 0.        , 0.01264648, ..., 0.        , 0.        ,
        0.        ],
       ...,
       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ]], dtype=float32)

Other Useful Methods and Properties

frame.masks is a list of masks (0 means masked). The number of arrays in this list is the number of planes in this frame. frame.masks[i] is identical to frame.components["MASK"][i] (where i is a plane index), but will be available (and filled with 1) if the data doesn't have a "MASK" component for some reason. Image frames additionally have frame.images, a list of unscaled images, for similar access to the "PIXEL" component. Vector frames have frame.enabled and frame.active_choice for access to the "ENABLED" and "ACTIVE_CHOICE" pixel components. These properties exist to make working with the advanced interface more convenient.

Vector frames have the following additional properties:

A frame can be indexed with plane indices. The result is a dict that maps component names to numpy arrays. No copying is performed, any changes in the array will be propagated to the original frame object and buffer object:

# multiply whole image buffer with a factor
for frame in buffer:
    for plane in frame:
        plane["PIXEL"] *= 2

Closing Sets

Like many other Python resources, e.g. file handles, Set objects can't be properly managed by the Python garbage collector. If you are opening and discarding large amounts of sets, it's a good idea to close them manually (or to use a context manager, but this won't work with multi-sets as they are simply lists).

# opening and closing a set manually:
s = lv.read_set("path/to/ImageSet")  # opening set
buffer = s[0]  # copying buffer to memory
s.close()  # set is closed, but the buffer is still available

# opening a multi-set and closing its subsets:
ms = lv.read_set("path/to/MultiSet")
s = ms[0]
buffer = s[0]
for s in ms:
    s.close()

# context manager:
with lv.read_set("path/to/ImageSet") as s:
    buffer = s[0]

You can check with s.closed if a set has already been closed. You should not read from a closed set.

For more examples on reading and plotting image and buffer data, please check the accompanying Jupyter notebooks, plots.ipynb and corr_map.ipynb.

Writing Sets

For writing a set to disk, use write_set. The first argument can be a Set object (use it if you want to convert a set produced in the latest DaVis to a set that can be read by an earlier DaVis version), a Python list of Buffer objects (use it if you want to read existing buffers, modify them and save them to a new set), or a Python list of numpy arrays (use them to create completely new data). Putting a new set in a project can be done by simply putting it in the project folder with your file manager.

Using write_set with Set Objects

lv.write_set(lv.read_set("data/2d3c"), "my_new_set")

Using write_set with Lists of Buffer Objects

When working with a multi-frame image, make sure that you modify all frames in a compatible way.

Example: adding noise to an image.

import numpy as np
import lvpyio as lv

s = lv.read_set("data/word_masked")

result = []
for buffer in s:
    for frame in buffer:
        image = frame.components["PIXEL"]
        plane = image[0]  # this is a 2D image, so only one plane
        noise = np.random.choice(5000, image.shape).astype(image.dtype)
        plane += noise
    result.append(buffer)

lv.write_set(result, "noise")

Example: modifying only the U component of a vector buffer.

import lvpyio as lv

buffer = lv.read_set("data/2d3c")[0]

for frame in buffer:
    for i in range(frame.choices):
        # U0, U1, U2, U3 because there are 4 choices:
        component = frame.components[f"U{i}"]
        for plane in component:
            plane *= 2  # multiply in-place

lv.write_set([buffer], "multiplied")

Using write_set with Lists of Numpy Arrays

Buffer can be represented by 2D arrays for 2D data and 3D arrays for 3D data. Buffers with multiple frames can be represented by lists of arrays.

Input Output
list of 2D arrays each buffer in set has one frame with one plane
list of 3D arrays each buffer in set has one frame with multiple planes
list of lists of 2D arrays each buffer in set has multiple frames, each frame has one plane
list of lists of 3D arrays each buffer in set has multiple frames, each frame has multiple planes

Allowed dtypes are:

For example, a list of lists with four 2D arrays of dtype numpy.dtype("uint8") will be written as a set of word images with four frames each, and a list with a single 3D-array with dtype lvpyio.vec3c will be written as a set with one 3C/3D vector volume.

Masked arrays are allowed. Their masks will be converted as well.

Example usage:

import numpy as np
import lvpyio as lv

# A set with two 3x3 images, each pixel contains a double-precision number
first = np.array([[1,2,3],[4,5,6],[7,8,9]], dtype=np.double)
second = np.array([[10,11,12],[13,14,15],[16,17,18]], dtype=np.double)
lv.write_set([first, second], "image")

# A set with one buffer, the buffer has 5 frames
frames = [first * i for i in range(5)]
lv.write_set([frames], "frames")

# Ten 2x2 vector fields with 2C vectors
vectors = np.array([[(1,2),(3,4)],[(5,6),(7,8)]], dtype=lv.vec2c)
s = []
for i in range(10):
    copy = vectors.copy()
    copy["u"] *= i
    copy["v"] *= i
    s.append(copy)

lv.write_set(s, "vectors")

# One 20 x 20 x 20 masked vector volume with 3C vectors
masked_vectors = np.ma.zeros((20, 20, 20), dtype=lv.vec3c)
masked_vectors.mask = True
lv.write_set([masked_vectors], "masked_vectors")

Reading and Writing Buffers Directly to Disk

read_buffer and write_buffer can be used to read a buffer directly from an .im7 or .vc7 file or write buffers to disk as .im7 or .vc7 files. write_buffer can be called with Buffer objects (use it if you want to read an existing buffer and modify it) or with numpy arrays (use it if you want to create a completely new buffer).

DaVis will ignore buffers that don't belong to a set, so you can't just put them in the project folder. Use the Import dialog for placing a new buffer into the project.

Using read_buffer

# reading a buffer from a path
buffer = lv.read_buffer("some/path/B00001.im7")

Using write_buffer

Just like with write_set, use 2D arrays for 2D data and 3D arrays for 3D data. Use lists of arrays for data with multiple frames.

Input Output
2D array one frame with one plane
3D array one frame with multiple planes
list of 2D arrays multiple frames, each has one plane
list of 3D arrays multiple frames, each has multiple planes

Allowed dtypes are the same as for write_set.

Reading Particle Fields

Particle fields (e.g. Shake-the-Box results) are a special kind of sets that should be read with the read_particles function:

>>> tr = lv.read_particles("data/particles/SingleTrack")
>>> tr
<Time-resolved particle field with 1001 time steps and 1 track>

Like image and vector sets, they consist of a .set file and a folder with the same name. Trying to read a particle field with read_set, or an image or vector set with read_particles, raises RuntimeError.

A particle field object has the following fields and methods:

Just like buffer-based sets, a particle field object can be used as a context manager, and has closed and close for manual management as well.

Time Steps

Particle fields don't consist of buffers. Instead, they consist of time steps. Time step indices are 1-based in DaVis but 0-based in Python, just like buffers in buffer-based sets have 1-based indices in DaVis and 0-based indices in Python. The number of time steps can be obtained with len:

>>> len(tr)
1001
>>> tr[0]
TimeStep(1 particle with track, 0 particles without track, 0 scalars)

Particle fields can be time-resolved or multi-pulse. In multi-pulse particle fields, time steps are called pulses and grouped together into blocks of equal length:

>>> mp = lv.read_particles("data/particles/SegmentationTrackingMp")
>>> mp
<2-pulse particle field with 5 blocks and 741 tracks>
>>> len(mp)
5
>>> mp[0]
MultiPulse((170 particles with track, 0 particles without track, 1 scalar),
(170 particles with track, 0 particles without track, 1 scalar))
>>> mp[0][0]
Pulse(170 particles with track, 0 particles without track, 1 scalar)

Because a particle field can be indexed and has a length, it can be used like any other sequence type:

for pulses in mp:
    for pulse in pulses:
        print(len(pulse.particles))

Timestamps of a multi-pulse particle field (as returned by the times method) are grouped together as well, so it's a 1D array for a time-resolved particle field but a 2D array for a multi-pulse particle field:

>>> tr.times()
array([0.   , 0.001, 0.002, ..., 0.998, 0.999, 1.   ])
>>> mp.times()
array([[0.   , 0.001],
       [0.005, 0.006],
       [0.01 , 0.011],
       [0.015, 0.016],
       [0.02 , 0.021]])
>>> tr.times().shape
(1001,)
>>> mp.times().shape
(5, 2)

For convenience, the type_id of a particle field contains the length of a pulse:

>>> tr.type_id
<Type.TIME_RESOLVED: 1>
>>> tr.type_id.value
1
>>> mp.type_id
<Type.DOUBLE_PULSE: 2>
>>> mp.type_id.value
2

Particles and Tracks

A time step contains particles. In the following example, the first time step contains one particle:

>>> tr[0].particles
array([(0, (0., 4000., 0., 1.))],
      dtype=[('track_id', '<u8'),
             ('particle', [('x', '<f4'), ('y', '<f4'), ('z', '<f4'),
                           ('intensity', '<f4')])])

In the other example, the first pulse of the first multi-pulse contains 170 particles:

>>> mp[0][0]
Pulse(170 particles with track, 0 particles without track, 1 scalar)
>>> len(mp[0][0].particles)
170

Each particle has coordinates (x, y, z) and intensity.

A particle can be part of a track. The number of tracks in the particle field can be accessed with track_count. Each track can be accessed by its ID. Track IDs are numbers from 0 to track_count - 1. Additionally, one can load all tracks at once with tracks. For large particle fields, loading all tracks can be both very slow and very memory-intensive, so if you only need a few tracks whose IDs you already know, use single_track.

>>> tr.track_count
1
>>> tr.tracks()
[Track(start=0, 1001 particles, 0 scalars)]
>>> track = tr.single_track(0)
>>> track
Track(start=0, 1001 particles, 0 scalars)
>>> track.start
0
>>> track.particles
array([(0.0000000e+00, 4000.      , 0.0000000e+00, 1.000e+00),
       (2.0000000e-06, 4000.      , 1.0000000e-06, 2.000e+00),
       (1.6000000e-05, 4000.      , 8.0000000e-06, 3.000e+00), ...,
       (1.9880239e+03,   23.952032, 9.9401196e+02, 9.990e+02),
       (1.9940060e+03,   11.988004, 9.9700299e+02, 1.000e+03),
       (2.0000000e+03,   -0.      , 1.0000000e+03, 1.001e+03)],
      dtype=[('x', '<f4'), ('y', '<f4'), ('z', '<f4'), ('intensity', '<f4')])
>>> track.particles.shape
(1001,)
>>> particle = track.particles[5] # fifth particle from the track
>>> particle
(0.00025, 3999.9995, 0.000125, 6.)
>>> particle["x"] # values can be accessed by name
0.00025
>>> x, y, z, i = particle # a particle can be destructured
>>> i
6.0

The same particle can be accessed via its time step:

>>> step = tr[5]
>>> step
TimeStep(1 particle with track, 0 particles without track, 0 scalars)
>>> step.particles
array([(0, (0.00025, 3999.9995, 0.000125, 6.))],
      dtype=[('track_id', '<u8'), ('particle', [('x', '<f4'), ('y', '<f4'),
      ('z', '<f4'), ('intensity', '<f4')])])
>>> track_id, particle = step.particles[0]
>>> track_id
0
>>> particle
(0.00025, 3999.9995, 0.000125, 6.)

If you know the time step and the track ID, you can access the particle directly by the track ID:

>>> step.particle_for_track(0)
(0.00025, 3999.9995, 0.000125, 6.)

Here is another example from a particle field with multiple tracks:

>>> mp.track_count
741
>>> track = mp.single_track(100)
>>> track
Track(start=0, 2 particles, 1 scalar)
>>> track.start
0
>>> len(track)
2

As we can see, track with ID 100 starts in the first block and has two particles.

For time-resolved data, start of a track is the index of its first time step (that is, the time step containing the first particle of the track). For 4-pulse data, start is a tuple of two indices: the index of its first multi-pulse block and the index of the pulse in the block. For 2-pulse data, start is the index of its first multi-pulse block (the index of the pulse in the block is not necessary: tracks for double-pulse data always start in the first of the two pulses).

IPR (iterative particle reconstruction, see FlowMaster Shake-the-Box product manual for details) data consists of particles without tracks. In this case, the particles array doesn't contain track IDs:

>>> ipr = lv.read_particles("data/particles/IPR")
>>> ipr
<2-pulse particle field with 3 blocks and 0 tracks>
>>> ipr[0][0]
Pulse(0 particles with track, 29 particles without track, 0 scalars)
>>> ipr[0][0].particles[0]
(2.3992462, 1.7948761, 8.602752, 0.5882383)

Scaling

Neither the particle coordinates nor the intensities are scaled. You will have to scale them manually using the corresponding scales stored in the particle field object. The scales can be obtained like this:

>>> tr.scales
Scales(x=Scale(slope=0.1, offset=-100.0, unit='mm', description='X'),
y=Scale(slope=-0.1, offset=200.0, unit='mm', description='Y'),
z=Scale(slope=0.1, offset=-50.0, unit='mm', description='Z'),
i=Scale(slope=1.0, offset=0.0, unit='', description=''))

The bounds (see above) are unscaled as well. They have the form ((min_x, min_y, min_z), (max_x, max_y, max_z)), where min_x and max_x should be scaled with the x scale, min_y and max_y with the y scale and so on.

Scalars

Particles may have scalar values attached to them (this can be the case in e.g. Tomographic Sizing experiments). For data with tracks, the scalar values may be read both from tracks (in which case they have a corresponding track_id value) and from the time steps / pulses:

>>> mp = lv.read_particles("data/particles/SegmentationTrackingMp")
>>> mp.scalar_scales
{'Particle diameter': Scale(slope=1.0, offset=0.0, unit='mm', description='')}
>>> pulse0 = mp[0][0]
>>> pulse0
Pulse(170 particles with track, 0 particles without track, 1 scalar)
>>> len(pulse0.scalars["Particle diameter"])
170
>>> pulse0.scalars["Particle diameter"][:2]
array([(0, 7.1812015), (1, 0.8418105)],
      dtype=[('track_id', '<u8'), ('scalar', '<f4')])
>>> pulse1 = mp[0][1]
>>> pulse1.scalars["Particle diameter"][:2]
array([(0, 7.249796), (1, 0.77885 )],
      dtype=[('track_id', '<u8'), ('scalar', '<f4')])
>>> track0 = mp.single_track(0)
>>> track0.scalars["Particle diameter"]
array([7.1812015, 7.249796 ], dtype=float32)
>>> track1 = mp.single_track(1)
>>> track1.scalars["Particle diameter"]
array([0.8418105, 0.77885  ], dtype=float32)

As you can see, the values are indeed the same.

For IPR data, the scalar values may only be read from time steps / pulses and don't have track IDs corresponding to them.

The names of all scalars can be obtained with .scalar_scales.keys().

Scalars are not to be confused with quantities like velocity components that one can see in DaVis. Velocities are not stored in the the particle field, so it's impossible for lvpyio to extract them. They are calculated by DaVis on the fly every time the particle field is displayed.

Attributes

A particle field has global attributes:

>>> mp.attributes["Timestamp"]
'2021-09-21T09:53:53,262Z'

In addition, each time step (for time-resolved data) or each multi-pulse block (for multi-pulse data) has attributes as well, usually copies of the buffer attributes of the image set that the particle set was calculated from:

>>> mp[0].attributes["TomoVolumeSubSetName"]
'Data/S00001'

Writing Particle Fields

You can read a particle field and write it back to disk:

pf = lv.read_particles("path/to/data")
lv.write_particles(pf, "copy")

This is useful for converting old particle fields created with DaVis 8 to the new, much faster format used by DaVis 10 (but it also means that particle fields written with lvpyio can't be read with DaVis 8). Note the important difference between write_particles, which converts old data to new format, and write_buffer/write_set, which convert new data to old format for maximum compatibility, so image and vector sets written with lvpyio can be read with older DaVis versions (however, some newer time-related attributes might not be understood).

Custom particle fields may be written as well by using a dictionary with the following keys:

Examples

Writing time-resolved tracks:

times = np.array([0.0,0.1,0.2,0.3,0.4])
tracks = [lv.Track(3, np.array([(0,1,2,3),(4,5,6,0)], dtype=lv.particle_t)),
          lv.Track(2, np.array([(1,2,3,4),(5,6,0,1)], dtype=lv.particle_t)),
          lv.Track(2, np.array([(2,3,4,5),(6,0,1,2),
                                (3,4,5,6)], dtype=lv.particle_t)),
          lv.Track(1, np.array([(0,1,2,3),(4,5,6,0)], dtype=lv.particle_t)),
          lv.Track(1, np.array([(1,2,3,4),(5,6,0,1),
                                (2,3,4,5)], dtype=lv.particle_t)),
          lv.Track(1, np.array([(6,0,1,2),(3,4,5,6),
                                (0,1,2,3),(4,5,6,0)], dtype=lv.particle_t)),
          lv.Track(0, np.array([(1,2,3,4),(5,6,0,1)], dtype=lv.particle_t)),
          lv.Track(0, np.array([(2,3,4,5),(6,0,1,2),(3,4,5,6)], dtype=lv.particle_t)),
          lv.Track(0, np.array([(0,1,2,3),(4,5,6,0),
                                (1,2,3,4),(5,6,0,1)], dtype=lv.particle_t)),
          lv.Track(0, np.array([(2,3,4,5),(6,0,1,2),(3,4,5,6),
                                (0,1,2,3),(4,5,6,0)], dtype=lv.particle_t))]

lv.write_particles(dict(times=times, tracks=tracks), "tracks_tr")

Writing 2-pulse IPR data:

times = np.array([(0.0,0.1),(0.2,0.3),(0.4,0.5)])
particles = [(np.array([(1,2,3,1),(2,3,1,2)], dtype=lv.particle_t),
              np.array([(3,1,2,3)], dtype=lv.particle_t)),
             (np.array([], dtype=lv.particle_t),
              np.array([(4,5,1,2),(3,4,5,1),(2,3,4,5),
                        (1,2,3,4),(5,1,2,3)], dtype=lv.particle_t)),
             (np.array([(-1,-2,-3,1)], dtype=lv.particle_t),
              np.array([(4,5,1,2),(3,4,5,1),(2,3,4,5),(1,2,3,4),
                        (5,1,2,3),(4,5,1,2),(3,4,5,1)], dtype=lv.particle_t))]
lv.write_particles(dict(times=times, particles=particles), "ipr_2p")

Writing 4-pulse tracks with scalars:

times = np.array([(0.0,0.1,0.2,0.3),(1.0.,1.1,1.2,1.3)])
tracks = [lv.Track(start=(0,0),
                   particles=np.array((2,1,2,1),
                                      (1,1,1,1),
                                      (4,3,2,1),
                                      (0,1,0,1), dtype=lv.particle_t),
                   scalars={"Particle diameter": np.array([1,2,3,2])}),
          lv.Track(start=(0,1),
                   particles=np.array((0,0,0,1),
                                      (0,1,1,1),
                                      (1,0,1,1), dtype=lv.particle_t),
                   scalars={"Particle diameter": np.array([3,2,1])}),
          lv.Track(start=(1,0),
                   particles=np.array((0,1,2,3),
                                      (4,5,6,7),
                                      (8,9,8,7),
                                      (6,5,4,3), dtype=lv.particle_t),
                   scalars={"Particle diameter": np.array([2,1,2,3])})]
scalar_scales = {"Particle diameter": lv.Scale(1, 0, "mm", "")}
d = dict(times=times, tracks=tracks, scalar_scales=scalar_scales)
lv.write_particles(d, "tracks_4p")

Absolute times are not guaranteed to be preserved during writing of custom particle fields, only the time differences between them are. We recommend to use times starting with 0.0 to avoid surprises:

>>> times = np.array([1,2,4,8], dtype=np.float64)
>>> lv.write_particles(dict(times=times, tracks=[]), "output")
>>> lv.read_particles("output").times()
array([0., 1., 3., 7.])

Scaling

All data (particle coordinates, particle intentities, and bounds) has to be unscaled, that is, all spatial data should be in voxel units. To make using data in other units more convenient, lvpyio provides two functions, unscale_particles and unscale_tracks, which take a particles array (or a list of tracks, respectively) and a Scales object and unscale the particles in-place (if you wish to preserve the scaled data, you'll have to copy it before applying the function). The same Scales object should be then put into the custom particle field dictionary that will be provided to the write_particles function.

Many DaVis operations, e.g. binning, expect ca. 500 .. 2000 voxels in each spatial direction. The scales for unscaling should be chosen accordingly.

Code Examples for Common Tasks

The examples in this section are runnable. We encourage you to try them out with the data that was included with the manual.

Image Boundaries in Scaled Units

import lvpyio as lv
s = lv.read_set("data/word_masked")
buffer = s[0]
frame = buffer[0]
height, width = frame.shape # height and width in pixel
x, y = frame.scales.x, frame.scales.y
print("left boundary:", x.offset, x.unit)
print("right boundary:", x.offset + width * x.slope, x.unit)
print("top boundary:", y.offset, y.unit)
print("bottom boundary:", y.offset + height * y.slope, y.unit)

# output:
# left boundary: -28.7315 mm
# right boundary: 33.93842 mm
# top boundary: 17.2664 mm
# bottom boundary: -29.37168 mm

Convert .png to DaVis Format

import matplotlib.pyplot as plt
import numpy as np
import lvpyio as lv

image = plt.imread("data/logo.png") # jpg or gif works, too
# converting RGBA to grayscale with weighted luminosity:
grayscale = np.dot(image[...,:3], [0.299, 0.587, 0.114])
buffer = grayscale.astype(np.float32)
lv.write_set([buffer], "logo_set") # a set with one buffer
lv.write_buffer(buffer, "logo_buffer.im7") # a single .im7 file

Convert Single Buffers to Sets

import lvpyio as lv
lv.write_set([lv.read_buffer("data/2d3c/B00001.vc7")], "my_new_set")

Quick Preview of 3D Vector Volumes

import matplotlib.pyplot as plt
import numpy as np
import lvpyio as lv

s = lv.read_set("data/3d3c")
buffer = s[0]
frame = buffer[0]

depth = len(frame)
height, width = frame.shape

planes = [frame.as_masked_array(plane=i) for i in range(depth)]

u = np.ma.array([plane["u"] for plane in planes])
v = np.ma.array([plane["v"] for plane in planes])
w = np.ma.array([plane["w"] for plane in planes])

# https://numpy.org/doc/stable/reference/generated/numpy.meshgrid.html
x, y, z = np.meshgrid(
    np.linspace(0, 1, depth),
    np.linspace(0, 1, height),
    np.linspace(0, 1, width),
    indexing="ij"
)

assert u.shape == x.shape == (depth, height, width)

fig = plt.figure()
ax = fig.gca(projection='3d')

# https://matplotlib.org/mpl_toolkits/mplot3d/tutorial.html
ax.quiver(x, y, z, u, v, w, length=0.02)

plt.show()

Find Out the Camera Name of Each Frame

import lvpyio as lv
s = lv.read_set("data/word_masked")
buffer = s [0]
names = [frame.attributes["CameraName"] for frame in buffer]
print(names)

Customize Plots

from matplotlib.ticker import AutoLocator
import matplotlib.pyplot as plt
import lvpyio as lv

s = lv.read_set("data/word_masked")
buffer = s[0]
figure = buffer.plot(show=False)
ax = figure.axes[0]

# restoring default ticks:
ax.xaxis.set_major_locator(AutoLocator())
ax.yaxis.set_major_locator(AutoLocator())
# adding a grid and a title:
ax.grid()
ax.set_title("Customized plot")
# cropping the image:
ax.set_xlim(500, 1000)
ax.set_ylim(300, 800)

# more can be found at https://matplotlib.org/api/axes_api.html

plt.show()

Convert Vector Data to HDF5

A script for converting vector sets or vector buffers to HDF5 (requires h5py) can be found in scripts/vector_to_hdf5.py. The script supports single-frame, single-plane (1D), single-choice vector data (so it's suitable for Strain results, but not for PIV results).

Limitations

As of now, only images, vector fields, and particle fields are supported. It can happen that a set of another format (e.g. profile data) can be read without errors, but this is an implementation detail you shouldn't rely on.

Subsets, e.g. calculation results or additional DVC data, have to be loaded manually by the user. For example, if you have read a DVC set called DVC_set and need access to its Data folder, you'll have to load DVC_set/Data/S00001.set, DVC_set/Data/S00002.set etc. as separate sets. An example (a script to resize a DVC recording to half size along all three axes, requires opencv-python and tqdm) can be found in scripts/dvc_resize.py.

Bayer RGB can be loaded, but has to be interpreted by the user. Other RGB image formats are not supported.

Writing multisets is not supported.

Writing custom particle fields with block attributes is not supported. The block attributes are supposed to be copies of the buffer attributes of the image set that the particle set was calculated from. If there is no such set, it doesn't make sense to have block attributes.

To facilitate support, all data written with lvpyio is automatically marked with a special _lvpyioVersion attribute. Writing data without this attribute is not possible.

Troubleshooting

Problem: on Windows, read_buffer throws file doesn't exist, but the file definitely exists.

Cause: Windows paths use backslash as delimiter, but Python treats it as an escape character. For example, in a path like "E:\2020\B00001.vc7", Python interprets \2 as a single character.

Solution: Possible workarounds are:

Problem: importing lvpyio in IPython or in a Jupyter notebook fails with ModuleNotFoundError.

Cause: this happens if you work inside a virtual environment but IPython / Jupyter is installed globally.

Solution: install IPython or Jupyter in your virtual environment.

Problem: lvpyio can't open .ims files.

Cause: .ims is not a standalone format like .im7. An .ims file can only be read as part of a set. Not even DaVis can import a separate .ims file.

Solution: let's assume you have a recording set named MyRecording and it contains .ims files, e.g. MyRecording/Frame0-0.ims. Don't try to open them. Instead, open MyRecording with lvpyio.read_set.

Miscellaneous

Before version 1.3.0, lvpyio was called lvreader.

Feedback

Please report bugs and confusing error messages to python@lavision.de.