Grains having disconnected regions.

[aditya-shukl]

CdTe.zip
Using topology. py in the attached file for disconnected regions leads to errors.

1. ~/cdte/axelscripts/topology.py in get_xray_endpoints(sample_polygon)
   183 Ymin = []
   184 Ymax = []
   --> 185 for poly in sample_polygon:
   186 xc, yc = poly.exterior.xy
   187 Xmin.append(np.min(xc))

TypeError: 'MultiPolygon' object is not iterable

and

in get_xray_endpoints(sample_polygon)
185 poly=sample_polygon
186 #for poly in sample_polygon:
--> 187 xc, yc = poly.exterior.xy
188 Xmin.append(np.min(xc))
189 Xmax.append(np.max(xc))

AttributeError: 'MultiPolygon' object has no attribute 'exterior'

[AxelHenningsson]

Hi! Could you please post a minimal failing example.

it looks like perhaps a version problem with shapely. The following code ray-traces as expected on my computer with shapely 1.7.0 :

import numpy as np
import matplotlib.pyplot as plt
import topology

if __name__ == "__main__":
    image = np.zeros((32,32))
    image[4:10,4:10]=1
    image[15:18,17:19]=1

    poly = topology.pixels_to_polygon(image, pixel_size=0.1, center=(0.5, 0.5))
    polygons = topology.voxels_to_polygon(image_stack=[image,image], pixel_size=0.1, center=(0.5, 0.5))

    angles = [45]
    ytrans = [0.0]
    zpos = [0]
    entry, exit, nhat, L, nsegs, bad_lines = topology.get_integral_paths(angles, ytrans, zpos, poly, nprocs=1, show_geom=False)

    print(entry)
    print(exit)
    print(L)
    print(poly)

    for p in poly:
        plt.plot(p.exterior.xy[0], p.exterior.xy[1])
    plt.plot([entry[0], exit[0]], [entry[1], exit[1]])
    plt.plot([entry[0+3], exit[0+3]], [entry[1+3], exit[1+3]])
    plt.show()

    plt.figure()
    plt.imshow(image)
    plt.show()

this represents two disconnected regions as:

these gets raytraced as:

[aditya-shukl]

Yes it was a version issues. I was running shapely 2.0.1

[aditya-shukl]

s3dxrd.utils.mesher.mesh_from_polygon also uses this poly.exterior attribute leading to error 'MultiPolygon' object has no attribute 'exterior'. Is the fix similar to what you did in topology.py

[AxelHenningsson]

Yeah, something like that. This implementation looks about right to me:

import numpy as np
import matplotlib.pyplot as plt
from shapely.geometry import Point as shapelyPoint
import topology

def mesh_from_polygon( w, multipolygon ):
    """Generate a quadratile 2d mesh over a Shapely polygon object.

    Args:
        w (:obj:`float`): Element side length.
        multipolygon (:obj:`Shapely MultiPolygon`): Polygon object to be meshed. The output mesh
            follows the :obj:`poly` coordinate system.

    Returns:
        (:obj:`numpy array`) The nodal coordinates of each element in the mesh.

    """

    xc = np.concatenate([[x for x in p.exterior.xy[0]] for p in multipolygon]).flatten()
    yc = np.concatenate([[y for y in p.exterior.xy[1]] for p in multipolygon]).flatten()

    mesh = []
    for x in np.arange(np.min(xc), np.max(xc) + w, w):
        for y in np.arange(np.min(yc), np.max(yc) + w, w):
            for poly in multipolygon:
                if poly.contains( shapelyPoint([(x+w/2., y+w/2.)]) ):
                    mesh.append([x, y, x+w, y, x+w, y+w, x, y+w])
                    break
    return np.array(mesh)


if __name__ == "__main__":
    image = np.zeros((50,50))
    image[4:10,4:10]=1
    image[15:18,17:19]=1
    image[17, 18]=0

    w = 0.01
    multipolygon = topology.pixels_to_polygon(image, pixel_size=w, center=(0.5, 0.5))

    mesh = mesh_from_polygon( w, multipolygon )
    fig,(ax1,ax2) = plt.subplots(1,2)
    for e in mesh:
        for i in range(3):
            ax2.plot([e[0+2*i], e[2+2*i]], [e[1+2*i], e[3+2*i]], 'k')
        ax2.plot([e[-2], e[0]], [e[-1], e[1]], 'k')

    ax1.imshow(np.fliplr(np.rot90(image)))
    plt.show()

It gives me somethign like this:

[AxelHenningsson]

Did this last feature solve the problem? Or are there further functions in need of patching?

Cheers
Axel

[aditya-shukl]

Yes, this solves the issue. Although, when we have many grains like in my case , the code takes a lot of time to go through all the layers. But I will close this issue.

Best
Adi

[aditya-shukl]

Hi again

I find some problems when I use the operations binary closing and binary opening for the grains. I can see that some of the grain shapes drastically change . Is there a work around ? Not using these operations leads to measurement_converter function failing. Here are the npy files of the grains that do not get converted to polygons nicely.
fail_case2.zip

Something like this .

Do you have any suggestions?

[aditya-shukl]

To clarify the left image is after using the binary operations and the right image is before using the binary operations.