UniformGrid.points setter not stable - need to deprecate

[lifehappenstoyou]

Hi everyone,

I ran into an issue with the UniformGrid's points.setter method giving me the error:

Traceback (most recent call last):
...
...
  File "/Users/Daniel/opt/anaconda3/lib/python3.7/site-packages/pyvista/core/grid.py", line 435, in _from_specs
    self.SetSpacing(xs, ys, zs)

TypeError: SetSpacing argument 1: only size-1 arrays can be converted to Python scalars

I use spacings in my point grid that are float and not integer since I need to have my grid in mm coordinates. It seems that the way the spacings are determined from the point data is not stable in this case. In this case the np.unique call that is supposed to pick one of the spacings determined from the dataset return an array, which then causes the above error.

Interestingly there is a remark in the code:
TODO: this needs to be tested (unique might return a tuple)

dx, dy, dz = np.unique(dx), np.unique(dy), np.unique(dz)

Well, I DID the testing now, is DOES return a tuple :-). In general it does not seem so easy to deal with. I suspect that this is caused by floating point issues. The idea to take the points, reconstruct the grid parameters they have been created with and then recreate the grid seems prone to errors in itself. But I do not have good solution for this. Anyone a good idea?

---

To Reproduce

import numpy as np
import pyvista as pv

du, dv = (0.2, 0.2)
n_u, n_v = (20, 20)
u, v = du * (np.arange(n_u) + 0.5), dv * (np.arange(n_v) + 0.5)
w = 0.0
pixel_grid_uvw = np.vstack(np.meshgrid(u, v, w))

# remove pixel matrix sub-structure by reshaping
pixel_centers = np.reshape(np.dstack(np.split(pixel_grid_uvw, 3, 0)), (-1, 3))
signal = np.ones(pixel_centers.size)

proj_image = pv.UniformGrid()
proj_image.points = pixel_centers
proj_image['signal'] = signal

System Information:

 Date: Wed May 13 12:12:33 2020 CEST

                OS : Darwin
            CPU(s) : 16
           Machine : x86_64
      Architecture : 64bit
               RAM : 32.0 GB
       Environment : Jupyter
        GPU Vendor : ATI Technologies Inc.
      GPU Renderer : AMD Radeon Pro 5500M OpenGL Engine
       GPU Version : 4.1 ATI-3.8.24

  Python 3.7.7 (default, Mar 26 2020, 10:32:53)  [Clang 4.0.1
  (tags/RELEASE_401/final)]

           pyvista : 0.24.2
               vtk : 8.2.0
             numpy : 1.18.1
           imageio : 2.8.0
           appdirs : 1.4.3
            scooby : 0.5.4
        matplotlib : 3.1.3
             PyQt5 : 5.9.2
           IPython : 7.13.0
             scipy : 1.4.1

  Intel(R) Math Kernel Library Version 2019.0.4 Product Build 20190411 for
  Intel(R) 64 architecture applications

[github-actions]

Hi and welcome! Thanks for posting your first issue in the PyVista project! Someone from @pyvista/developers will chime in before too long. If your question is support related, it may be automatically transferred to https://github.com/pyvista/pyvista-support

[GuillaumeFavelier]

Thanks for reporting @lifehappenstoyou ! I'll try to reproduce.

[GuillaumeFavelier]

As expected, I encounter the same fate:

Traceback (most recent call last):
  File "/tmp/bug.py", line 15, in <module>
    proj_image.points = pixel_centers
  File "/home/guillaume/source/pyvista/pyvista/core/grid.py", line 472, in points
    self._from_specs((nx,ny,nz), (dx,dy,dz), (ox,oy,oz))
  File "/home/guillaume/source/pyvista/pyvista/core/grid.py", line 435, in _from_specs
    self.SetSpacing(xs, ys, zs)
TypeError: SetSpacing argument 1: only size-1 arrays can be converted to Python scalars

The idea to take the points, reconstruct the grid parameters they have been created with and then recreate the grid seems prone to errors in itself.

If you know those informations, there is still the option to set dimensions, spacing and origin instead of setting the points directly (as a workaround).

[GuillaumeFavelier]

Using a version of unique that supports a tolerance factor might help in this situation:

def _unique(x):
    def eclose(a,b,rtol=1.0000000000000001e-05, atol=1e-08):
        return np.abs(a - b) <= (atol + rtol * np.abs(b))

    y = x.flat.copy()
    y.sort()
    ci = 0
    U = np.empty((0,),dtype=y.dtype)

    while ci < y.size:
        ii = eclose(y[ci],y)
        mi = np.max(ii.nonzero())
        U = np.concatenate((U,[y[mi]])) 
        ci = mi + 1
    return U

From: https://stackoverflow.com/questions/5426908/find-unique-elements-of-floating-point-array-in-numpy-with-comparison-using-a-d

Also, in your code snippet, you use:

signal = np.ones(pixel_centers.size)

I think you mean:

signal = np.ones(len(pixel_centers))

size will return the total number of entries and len the effective number of rows (corresponding to the number of points).

[lifehappenstoyou]

Hi Guillaume,

yes, I think introducing a tolerance in unique might be a way to solve the issue. I will for now try to build the UniformGrid from dimensions and spacings. For my immediate use there should be no problem. And propably doing it this way anyways might be best and safest for a UniformGrid. The question would be if this constructor is needed then in the first place for a UniformGrid? I originally tried a StructuredGrid which is a different story then. There are slight floating point differences do not matter when setting it up as you have to specify all points anyways.

Thanks for pointing out my mistake with pixel_centers.size, you are right, it must be len(pixel_centers) instead. Since the bug appeared before I did not care to much in the minimum example.

best Daniel

[GuillaumeFavelier]

What is your opinion on this @banesullivan ?

[banesullivan]

IMO (as the person who originally implemented all of this code for UniformGrid), we should take out the points setter on UniformGrid meshes (honestly, I thought I did this a long time ago because of a lot of issues with it).

UniformGrids are implicitly defined by three parameters: 1) origin, 2) spacing, 3) dimensions/extent. Thus, we should not support letting users explicitly set the spatial reference. The issue outlined here is one such reason: there's no great way to guarantee that the points passed by the user make a perfect uniform grid of points - and thus we'd have to write TON of checks to properly figure out the spacing, origin, extent all while also gaurunteeing the points do not have any bit of a curvilinear structure. They could be structured and curvilinear and this setter will fail miserably.

If you want to be able to explicitly set the points of a grid mesh, then you should use the StructuredGrid class, not the UniformGrid class.

Also, having the points setter is misleading: you cannot dynamically adjust the points of a UniformGrid like you can other PyVista mesh types. and when you do, the actual points generated by the points property getter will not match the points passed by the user to the property setter.

[banesullivan]

SO my recommendation to you, @lifehappenstoyou (great username!), is to switch to using the StructuredGrid class unless you have a compelling reason to use the UnifromGrid class (the only thing I can think of is volume rendering support, in which cases we have work-arounds)

[lifehappenstoyou]

Thanks for your input on this, @banesullivan. After discussing it above, not providing a point setter method for this case at all seemed to me a reasonable conclusion as well. Setting the points explicitely does not really make sense. I just ran into this problem, because I had the point coordinates at hand anyways, so I thought, why not use it? But personally I can work around this, so no real issue for me and there should not be any issue for anybody else for the reasons you pointed out once this method is taken out.

Thanks to all contributors for developing pyvista. I just started using it for my project and it really saved me from digging through all the rather clumsy vtk code in most cases! Awesome work!

[banesullivan]

I'm going to keep this open as a ticket to deprecate that setter and try to do it this weekend.

I just ran into this problem, because I had the point coordinates at hand anyways, so I thought, why not use it?

That's exactly the logic I had when implementing as I was using UniformGirds on a weird little project and just threw this feature into PyVista when implementing the UniformGrid class... in. hindsight, not a great idea.

Otherwise, it's awesome to hear you are enjoying and benefiting from PyVista in your work!