CODEWITHSUNDEEP
pythonnumpyvtk
fabda01
fabda01

Reputation: 3763

Calculate the distance between two .vtk files in (Python)

Applying a transformation in the following .vtk file:

im1.vtk:

# vtk DataFile Version 3.0
vtk output
ASCII
DATASET POLYDATA
POINTS 10 float
-61.2 40.8 0.0 
-55.3 39.3 0.0 
-49.2 39.3 0.0 
-43.2 40.4 0.0 
-37.3 42.1 0.0 
67.6 44.3 0.0 
63.4 49.8 0.0 
57.7 53.6 0.0 
51.0 55.3 0.0 
44.1 55.2 0.0

Results in another .vtk file:

im2.vtk:

# vtk DataFile Version 4.0
vtk output
ASCII
DATASET POLYDATA
POINTS 10 float
-63.4695 36.4645 0
-57.3647 35.9114 0
-51.1496 36.6507 0
-45.102 38.259 0
-39.2082 40.2851 0
69.7562 40.4176 0
64.6497 45.5255 0
58.1449 49.2956 0
50.8203 51.2899 0
43.4762 51.6839 0

What is the fastest way to find the distance (e.g.: Euclidean) between each point for each file and store it in a list [d1,d2,...,d10] or numpy array? (This is just a sample, the real .vtk file contains 300 points. So the list should be [d1,d2,...,d300]). Is there already a method in VTK which does that?

Example:

The distance d1 between the first point should be:

d1=sqrt((-61.2+63.4695)**2+(40.8-36.4645)**2+(0-0)**2)

EDIT:

The code I could come up with so far is as follows:

import numpy as np

with open('im1.vtk', 'rt') as vtk1:
   vtk_list1 = vtk1.readlines()
with open('im2.vtk', 'rt') as vtk2:
   vtk_list2 = vtk2.readlines()
dist_array = np.array([])
for i in range(5,14):
   landmark1 = np.asarray(vtk_list1[i].split(),dtype=np.float16)
   landmark2 = np.asarray(vtk_list2[i].split(),dtype=np.float16)
   dist = np.sqrt(np.sum((landmark1-landmark2)**2))
   dist_array = np.append(dist_array,dist)

Is there an optimal one that is faster and cleaner?

Upvotes: 1

Views: 571

Answers (1)

Severin Pappadeux
Severin Pappadeux

Reputation: 20130

Is there an optimal one that is faster and cleaner?

Sure - instead of filling rows, fill columns and use full power of numpy

Along the lines

import sys
import numpy as np

def read_vtk(fname):
    """
    """
    with open(fname, 'rt') as f:
        lines = f.readlines()

        l = lines[4]
        s = l.split()
        n = int(s[1])

        x = np.empty(n, dtype=float)
        y = np.empty(n, dtype=float)
        z = np.empty(n, dtype=float)

        for k in range(0, n):
            l = lines[5 + k]
            s = l.split()
            x[k] = float(s[0])
            y[k] = float(s[1])
            z[k] = float(s[2])

        return (x, y, z)

    return None

def main(fa, fb):
    """
    """
    nax, nay, naz = read_vtk(fa)
    nbx, nby, nbz = read_vtk(fb)

    d = np.power(nax-nbx, 2) + np.power(nay-nby, 2) + np.power(naz-nbz, 2)

    return np.sqrt(d)

if __name__ == "__main__":
    d = main("a.vtk", "b.vtk")

    print(d)

    sys.exit(0)

Upvotes: 1

Related Questions