Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Voxelization routine #44

Closed
Kresa96 opened this issue Sep 8, 2021 · 17 comments · Fixed by #45, #46 or #47
Closed

Voxelization routine #44

Kresa96 opened this issue Sep 8, 2021 · 17 comments · Fixed by #45, #46 or #47
Assignees

Comments

@Kresa96
Copy link

Kresa96 commented Sep 8, 2021

Very nice work!
We try to use your package MicroStructPy in our research group to built polycristalline electrolytes for impedance simulation. We usually use numpy arrays in which each grain is represented by several array entries (very similar to voxels, or images in 2D). The grains are distinguishable by different integer values in the matrix. In this way it's easy to modify the structure and use it for several different purposes (as an example for our simulations).

Example for a 2D structure consisting of 2 grains:
1111112222
1111112222
1111222222
1112222222

I think it would be very nice if you could implement a function which can transform the 2d and 3d polygon mesh into numpy arrays.

Looking forward to your reply!

@kip-hart kip-hart self-assigned this Sep 8, 2021
@kip-hart
Copy link
Owner

kip-hart commented Sep 8, 2021

Hi thank you! I'm glad to hear you're using MicroStructPy. I understand what you're asking for, where there is a grid of points and the value assigned to each point is the grain number, then these values are stored in a matrix.

I'm not familiar with voxelization approaches. My first thought is to create a list of points then test if point "i" is within grain "j", but that would be very slow. My second thought is to put bounding boxes around each grain, then only test if point "i" is within grain "j" if the point is within the bounding box. This is better, but probably still slow.

Do you have suggestions for voxelization routines?

@Kresa96
Copy link
Author

Kresa96 commented Sep 8, 2021

Thanks for the fast anwer!
No I am also not very familiar with those methods. But I saw that a voxelization routine for ellipsoidal grains is already implemented in 'Kanapy', a Python package from 'mrgprasad', which was also made to create microstructures in Python. Maybe all you need is in their documentation.

In 2D I was able to translate the polymesh.plot in a very 'dirty' way into a numpy array:

polymesh.plot(facecolors = "white", edgecolor = "k")
plt.axis('off')
buf = io.BytesIO()
plt.savefig(buf, format="png", bbox_inches='tight' , pad_inches = 0, dpi=500)
buf.seek(0)
img_arr = np.frombuffer(buf.getvalue(), dtype=np.uint8)
buf.close()
img = cv2.imdecode(img_arr, 0)

s = [[1,1,1],[1,1,1],[1,1,1]]
labeled_mask, num_labels = ndimage.label(img, structure = s)
labeled_mask = expand_labels(labeled_mask, distance = 20)

But in 3D this dirty procedure will not work.
If you are able to implement something like this, we will probably use your package for lots of stuff.

@kip-hart
Copy link
Owner

kip-hart commented Sep 8, 2021

Ok got it, I'll check out Kanapy. I've started a branch to add a RasterMesh class. The pseudo-code for creating the mesh is

# 1. Create node and element grids
# 2. Compute element centers
# 3. Remove elements outside domain
# 4. For each region:
# A. Create a bounding box
# B. Isolate element centers with box
# C. For each facet, remove centers on the wrong side
# D. Assign remaining centers to region
# 5. Combine regions of the same seed number (remove voids)
# 6. Define remaining facets, inherit their attributes

This mesh class would have the same attributes as TriMesh. To create the array of grain numbers, my pseudo-code is

# 1. Create array full of -1 values
# 2. Convert 1st node of each element into array indices
# 3. For each element: populate array with element attributes

I'd also add functions to write these meshes out in VTK and Abaqus formats.

I can keep you posted with updates. How urgent is the need?

@Kresa96
Copy link
Author

Kresa96 commented Sep 8, 2021

Nice!
Yeah I think it is a nice idea to implement the output of these files!
Keep me updated!
The need is very urgent 😁

@kip-hart kip-hart mentioned this issue Sep 9, 2021
3 tasks
@kip-hart
Copy link
Owner

kip-hart commented Sep 9, 2021

I've added details in #45 for how this is implemented. You can get the desired array with a script that first creates/opens a PolyMesh, then calls RasterMesh.from_polymesh. The raster mesh created by this function has a method as_array. Array indices are 0 in the lower-left corner. The functions work in both 2D and 3D.

Upgrading to v1.5.0 will have these features. Please let me know if you encounter any issues.

kip-hart added a commit that referenced this issue Sep 9, 2021
@Kresa96
Copy link
Author

Kresa96 commented Sep 9, 2021

Thank you very much for the fast implementation!!
I got a view problems:
In 2D the plot() method does not work in my case. If I plot the matrix with plt.imshow() for several values of the argument 'mesh_size' I get anomalies in the matrix:
image

In 3D voxelization does not work at all in my case:
Input:
image
Output:
image

@kip-hart kip-hart reopened this Sep 9, 2021
@kip-hart
Copy link
Owner

kip-hart commented Sep 9, 2021

Sorry to hear that. Could you attach the 2D and 3D polymesh.txt files? You can also email them to support@microstructpy.org if that doesn't work. I'll debug with those; it's likely a floating point error for the 2D case.

For the .plot() function, this is the default behavior when there are no options. If you assign a color to each seed, then .plot(color=seed_colors, index_by='attribute') should produce the desired result. You can also pass other matplotlib arguments, like .plot(..., edgecolor='k', lw=0.5). All the non-MicroStructPy keyword args are passed through to matplotlib. More details at https://docs.microstructpy.org/en/latest/api/meshing/trimesh.html#microstructpy.meshing.TriMesh.plot

@Kresa96
Copy link
Author

Kresa96 commented Sep 9, 2021

Yes you are right the implemented 2D plot works for me:
image
But it seems like that the array gets turned and the anomalies arrise when as_array() gets applied to the grid.

I will send the files to the mail!
Thanks!

@kip-hart
Copy link
Owner

kip-hart commented Sep 9, 2021

Since (0,0) is in the lower-left corner, plotting with plt.imshow(..., origin='lower') should prevent the turning. You can also use the extent=(left, right, bottom, top) option to scale the axes into physical space, rather than matrix coordinates

I'll debug the 2D plaid effect and 3D issues now and update you with results.

@kip-hart
Copy link
Owner

kip-hart commented Sep 9, 2021

I've cleared up the plaid issue in 2D. It was a floating point error. Working on the 3D case now.

Screen Shot 2021-09-09 at 11 44 52 AM

@Kresa96
Copy link
Author

Kresa96 commented Sep 9, 2021

Perfect! Thank you very much!

@kip-hart
Copy link
Owner

kip-hart commented Sep 9, 2021

3D is now working. I'll update the package shortly to v1.5.1 with the fixes.

Screen Shot 2021-09-09 at 12 58 02 PM

@kip-hart kip-hart mentioned this issue Sep 9, 2021
3 tasks
@Kresa96
Copy link
Author

Kresa96 commented Sep 9, 2021

Thanks a lot I am looking forward to the new version!

kip-hart added a commit that referenced this issue Sep 9, 2021
@kip-hart
Copy link
Owner

kip-hart commented Sep 9, 2021

v1.5.1 has been released. Please let me know if you encounter any issues.

@Kresa96
Copy link
Author

Kresa96 commented Sep 9, 2021

Sorry but the coloring doesnt work in my case:
image
Error code:
ValueError Traceback (most recent call last)
in
30 colors_list.append(random_color)
31
---> 32 rmesh.plot(index_by='attribute', color=['blue', 'red'])

~\anaconda3\envs\micropy4\lib\site-packages\microstructpy\meshing\trimesh.py in plot(self, index_by, material, loc, **kwargs)
1344 for i, x in enumerate(pt):
1345 grids[i][tuple(pt_tup)] = x
-> 1346 ax.voxels(*grids, inds >= 0, **plt_kwargs)
1347
1348 # Add legend

~\anaconda3\envs\micropy4\lib\site-packages\mpl_toolkits\mplot3d\axes3d.py in voxels(self, facecolors, edgecolors, shade, lightsource, *args, **kwargs)
2907 )
2908
-> 2909 poly = art3d.Poly3DCollection(
2910 faces, facecolors=facecolor, edgecolors=edgecolor, **kwargs)
2911 self.add_collection3d(poly)

~\anaconda3\envs\micropy4\lib\site-packages\mpl_toolkits\mplot3d\art3d.py in init(self, verts, zsort, *args, **kwargs)
671 and _edgecolors properties.
672 """
--> 673 super().init(verts, *args, **kwargs)
674 self.set_zsort(zsort)
675 self._codes3d = None

~\anaconda3\envs\micropy4\lib\site-packages\matplotlib\collections.py in init(self, verts, sizes, closed, **kwargs)
1121 Forwarded to .Collection.
1122 """
-> 1123 Collection.init(self, **kwargs)
1124 self.set_sizes(sizes)
1125 self.set_verts(verts, closed)

~\anaconda3\envs\micropy4\lib\site-packages\matplotlib\cbook\deprecation.py in wrapper(*inner_args, **inner_kwargs)
409 else deprecation_addendum,
410 **kwargs)
--> 411 return func(*inner_args, **inner_kwargs)
412
413 return wrapper

~\anaconda3\envs\micropy4\lib\site-packages\matplotlib\collections.py in init(self, edgecolors, facecolors, linewidths, linestyles, capstyle, joinstyle, antialiaseds, offsets, transOffset, norm, cmap, pickradius, hatch, urls, offset_position, zorder, **kwargs)
211
212 self._path_effects = None
--> 213 self.update(kwargs)
214 self._paths = None
215

~\anaconda3\envs\micropy4\lib\site-packages\matplotlib\artist.py in update(self, props)
996 raise AttributeError(f"{type(self).name!r} object "
997 f"has no property {k!r}")
--> 998 ret.append(func(v))
999 if ret:
1000 self.pchanged()

~\anaconda3\envs\micropy4\lib\site-packages\matplotlib\collections.py in set_color(self, c)
749 For setting the edge or face color individually.
750 """
--> 751 self.set_facecolor(c)
752 self.set_edgecolor(c)
753

~\anaconda3\envs\micropy4\lib\site-packages\mpl_toolkits\mplot3d\art3d.py in set_facecolor(self, colors)
792
793 def set_facecolor(self, colors):
--> 794 PolyCollection.set_facecolor(self, colors)
795 self._facecolor3d = PolyCollection.get_facecolor(self)
796

~\anaconda3\envs\micropy4\lib\site-packages\matplotlib\collections.py in set_facecolor(self, c)
778 """
779 self._original_facecolor = c
--> 780 self._set_facecolor(c)
781
782 def get_facecolor(self):

~\anaconda3\envs\micropy4\lib\site-packages\matplotlib\collections.py in _set_facecolor(self, c)
762 except AttributeError:
763 pass
--> 764 self._facecolors = mcolors.to_rgba_array(c, self._alpha)
765 self.stale = True
766

~\anaconda3\envs\micropy4\lib\site-packages\matplotlib\colors.py in to_rgba_array(c, alpha)
339 return np.zeros((0, 4), float)
340 else:
--> 341 return np.array([to_rgba(cc, alpha) for cc in c])
342
343

~\anaconda3\envs\micropy4\lib\site-packages\matplotlib\colors.py in (.0)
339 return np.zeros((0, 4), float)
340 else:
--> 341 return np.array([to_rgba(cc, alpha) for cc in c])
342
343

~\anaconda3\envs\micropy4\lib\site-packages\matplotlib\colors.py in to_rgba(c, alpha)
187 rgba = None
188 if rgba is None: # Suppress exception chaining of cache lookup failure.
--> 189 rgba = _to_rgba_no_colorcycle(c, alpha)
190 try:
191 _colors_full_map.cache[c, alpha] = rgba

~\anaconda3\envs\micropy4\lib\site-packages\matplotlib\colors.py in _to_rgba_no_colorcycle(c, alpha)
263 raise ValueError(f"Invalid RGBA argument: {orig_c!r}")
264 if len(c) not in [3, 4]:
--> 265 raise ValueError("RGBA sequence should have length 3 or 4")
266 if not all(isinstance(x, Number) for x in c):
267 # Checks that don't work: map(float, ...), np.array(..., float) and

ValueError: RGBA sequence should have length 3 or 4

@kip-hart kip-hart reopened this Sep 9, 2021
@kip-hart kip-hart mentioned this issue Sep 9, 2021
3 tasks
@kip-hart
Copy link
Owner

kip-hart commented Sep 9, 2021

This will be fixed in v1.5.2. For 3D, you'll need to use the keyword facecolors= instead of color=.

Specifying facecolors=['blue', 'red'] only gives colors to seeds 0 and 1. In this case, other seeds will be rendered invisible:

Screen Shot 2021-09-09 at 4 27 13 PM

What I've done in the past is create a color list like:

n = max(rmesh.element_attributes) + 1
colors = ['C' + str(i % 10) for i in range(n)]
rmesh.plot(index_by='attribute', facecolors=colors)
plt.show()

This will produce a plot like:

Screen Shot 2021-09-09 at 4 31 19 PM

You can also create a color list by sampling colormaps. Additionally, you can set grains of the same material phase to have the same color.

Hopefully this helps. Version 1.5.2 will be available shortly.

kip-hart added a commit that referenced this issue Sep 9, 2021
@Kresa96
Copy link
Author

Kresa96 commented Sep 10, 2021

It works! Great work and very fast implementation!
Thank you very much!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
2 participants