pnts

Author: CR10-dv

geospatial_learn/learning.py

@@ -51,6 +51,9 @@
 
 import pandas as pd
 import simpledbf
+from plyfile import PlyData
+from pyntcloud import PyntCloud
+import open3d as o3d
 
 
 gdal.UseExceptions()
@@ -1439,140 +1442,121 @@ def get_training(inShape, inRas, bands, field, outFile = None):
 
     return outData, rejects
 
+def ply_features(incld, outcld=None):
 
-def get_training_point(inShape, inRas, bands, field):
-    """ Collect training as a np array for use with create model function using 
-          point data
-          
-    Parameters
-    --------------
-        
-    inShape : string
-              the input shapefile - must be esri .shp at present
-        
-    inRas : string
-            the input raster from which the training is extracted
-        
-    bands : int
-            no of bands
-        
-    field : string
-            the attribute field containing the training labels
+    """ 
+    Calculate point cloud features and write to file
     
-    outFile : string (optional)
-              path to the training file saved as joblib format (eg - 'training.gz')
     
-    Returns
-    ---------------------
+    Parameters 
+    ----------- 
     
-    A tuple containing:
-    -np array of training data
-    -list of polygons with invalid geometry that were not collected 
+    incld: string
+              the input point cloud
+        
+    outcld: string
+               the output point cloud
+         
 
+    """  
 
-    UNFINISHED DO NOT USE
-    
+    pcd = PyntCloud.from_file(incld)
+    #o3d.io.read_point_cloud(incld)
+
+    #pcd.estimate_normals()
 
-    """
-    #t0 = time()
-    outData = list()
-    print('Loading & prepping data')    
-    raster = gdal.Open(inRas)
-    shp = ogr.Open(inShape)
-    lyr = shp.GetLayer()
-    labels = np.arange(lyr.GetFeatureCount())
-    rb = raster.GetRasterBand(1)
-    rgt = raster.GetGeoTransform()
-    mem_drv = ogr.GetDriverByName('Memory')
-    driver = gdal.GetDriverByName('MEM')  
-    rejects = []     
+    #cloud = PyntCloud.from_instance("open3d", 
+
+
+    pProps =['anisotropy', "curvature", "eigenentropy", "eigen_sum", "linearity",
+             "omnivariance", "planarity", "sphericity"]#, "inclination_deg",
+            # "inclination_rad", "orientation_deg", "orientation_rad"]
+    #, "HueSaturationValue",#"RelativeLuminance"," RGBIntensity"]
+
+    k_neighbors = pcd.get_neighbors(k=30)
+    eigenvalues = pcd.add_scalar_field("eigen_values", k_neighbors=k_neighbors)
+
+    [pcd.add_scalar_field(p, ev=eigenvalues) for p in pProps]
 
-    print('getting points')
-    for label in tqdm(labels):
-        #print(label)
-        feat = lyr.GetFeature(label)
-        if feat == None:
-            print('no geometry for feature '+str(label))
-            continue
-        iD = feat.GetField(field)
-        geom = feat.GetGeometryRef()
-        mx,my=geom.GetX(), geom.GetY()  #coord in map units
+    if outcld == None:
+        pcd.to_file(incld)
+    else:
+        pcd.to_file(outcld)
 
-        #Convert from map to pixel coordinates.
-        #Only works for geotransforms with no rotation.
-        px = int((mx - gt[0]) / gt[1]) #x pixel
-        py = int((my - gt[3]) / gt[5]) #y pixel
+def get_training_ply(incld, label_field="scalar_label", outFile=None):
 
-        structval=rb.ReadRaster(px, py,1,1,buf_type=gdal.GDT_UInt16) #Assumes 16 bit int aka 'short'
-        intval = struct.unpack('h' , structval) #use the 'short' format code (2 bytes) not int (4 bytes)
+    """ 
+    Get training from a point cloud
     
-        print(intval[0])
-        # Get raster georeference info
-            
-#        src_offset = bbox_to_pixel_offsets(rgt, geom)
-#        
-#        
-#        # calculate new geotransform of the feature subset
-#        new_gt = (
-#        (rgt[0] + (src_offset[0] * rgt[1])),
-#        rgt[1],
-#        0.0,
-#        (rgt[3] + (src_offset[1] * rgt[5])),
-#        0.0,
-#        rgt[5])
-            
-            
-        # Create a temporary vector layer in memory
-        mem_ds = mem_drv.CreateDataSource('out')
-        mem_layer = mem_ds.CreateLayer('poly', None, ogr.wkbPolygon)
-        mem_layer.CreateFeature(feat.Clone())
-
-        # Rasterize it
-        rvds = driver.Create('', src_offset[2], src_offset[3], 1, gdal.GDT_Byte)
-        rvds.SetGeoTransform(new_gt)
-        gdal.RasterizeLayer(rvds, [1], mem_layer, burn_values=[1])
-        rv_array = rvds.ReadAsArray()
-        
-        # Mask the source data array with our current feature
-        # we take the logical_not to flip 0<->1 to get the correct mask effect
-        # we also mask out nodata values explictly
-            
+    
+    Parameters 
+    ----------- 
+    
+    incld: string
+              the input point cloud
+    
+    label_field: string
+              the name of the field representing the training points which must
+              be positive integers
+                
+    outFile: string
+               path to training array to be saved as .gz via joblib
+         
 
+    """  
+    # TODO Clean up lack of loops funcs to do stuff
+    # classify ply TODO also
+    # convert Open3D.o3d.geometry to necessary vars
+    pcd = o3d.io.read_point_cloud(incld)
+    xyz = np.asarray(pcd.points)
+    cols = np.asarray(pcd.colors)
+    norms = np.asarray(pcd.normals)
+    pf = PlyData.read(incld)
+    
+#    pProps =['anisotropy', 'curvature', "eigenentropy", "eigen_sum",
+#             "linearity", "omnivariance", "planarity", "sphericity"]
+    
+    # doesn't matter if this is a float
+    label = np.array(pf.elements[0].data['scalar_label'], dtype='float64')
+    a = np.array(pf.elements[0].data['anisotropy(26)'], dtype='float64')
+    c = np.array(pf.elements[0].data["curvature(26)"], dtype='float64')
+    et = np.array(pf.elements[0].data["eigenentropy(26)"], dtype='float64')
+    es = np.array(pf.elements[0].data["eigen_sum(26)"], dtype='float64')
+    l = np.array(pf.elements[0].data["linearity(26)"], dtype='float64')
+    om = np.array(pf.elements[0].data["linearity(26)"], dtype='float64') 
+    pl = np.array(pf.elements[0].data["omnivariance(26)"], dtype='float64')
+    sp = np.array(pf.elements[0].data["sphericity(26)"], dtype='float64')
+
+    label.shape = (label.shape[0], 1)
+    a.shape=(label.shape[0], 1)
+    c.shape=(label.shape[0], 1)
+    et.shape=(label.shape[0], 1)
+    es.shape=(label.shape[0], 1)
+    l.shape=(label.shape[0], 1)
+    om.shape=(label.shape[0], 1)
+    pl.shape=(label.shape[0], 1)
+    sp.shape=(label.shape[0], 1)
+    
+    final = np.hstack((label, xyz, cols, norms, a, c, et, es, l, om, pl, sp))
+    
+    # all these could be dumped now
+    del xyz, cols, norms, pf,  a, c, et, es, l, om, pl, sp
+    
+    # prep for sklearn
+    X_train = final[final[:,0] >= 0]
 
-        
-        rb = raster.GetRasterBand(1)
-        src_array = rb.ReadAsArray(src_offset[0], src_offset[1], src_offset[2],
-                           src_offset[3])
-        if np.shape(src_array) is ():
-            rejects.append(label)
-            continue
-        # Read raster as arrays
-        for band in range(1,bands+1): 
-            
-            rb = raster.GetRasterBand(band)
-            src_array = rb.ReadAsArray(src_offset[0], src_offset[1], src_offset[2],
-                               src_offset[3])
-            if src_array is None:
-                src_array = rb.ReadAsArray(src_offset[0]-1, src_offset[1], src_offset[2],
-                                           src_offset[3])
-
-            masked = np.ma.MaskedArray(src_array, 
-                                       mask=np.logical_or(src_array == 0,
-                                                          np.logical_not(rv_array)))
-            
-
-            datafinal = masked.flatten()
+    # Remove non-finite values
+    X_train = X_train[np.isfinite(X_train).all(axis=1)]
+    
+    # y labels
+    #y_train = X_train[:,0] 
+    
+    if outFile != None:
+        jb.dump(X_train, outFile, compress=2)
+    
+    return X_train
+    
 
-            if band == 1:
-                X = np.zeros(shape = (datafinal.shape[0], bands+1))
-            X[:,0] = iD
-            
-            X[:,band] = datafinal
-        #print(label,fieldval,xcount, ycount)   
-        outData.append(X)
-    outData = np.asarray(outData)
-    outData = np.concatenate(outData).astype(None)
-    return outData, rejects
 
 
 def rmse_vector_lyr(inShape, attributes):

geospatial_learn/utilities.py

@@ -1893,133 +1893,6 @@ def rotate_im(image, angle):
     return image
 
 
-def get_corners(bboxes):
-    
-    """Get corners of bounding boxes
-    
-    Parameters
-    ----------
-    
-    bboxes: numpy.ndarray
-        Numpy array containing bounding boxes of shape `N X 4` where N is the 
-        number of bounding boxes and the bounding boxes are represented in the
-        format `x1 y1 x2 y2`
-    
-    returns
-    -------
-    
-    numpy.ndarray
-        Numpy array of shape `N x 8` containing N bounding boxes each described by their 
-        corner co-ordinates `x1 y1 x2 y2 x3 y3 x4 y4`      
-        
-    """
-    width = (bboxes[:,2] - bboxes[:,0]).reshape(-1,1)
-    height = (bboxes[:,3] - bboxes[:,1]).reshape(-1,1)
-    
-    x1 = bboxes[:,0].reshape(-1,1)
-    y1 = bboxes[:,1].reshape(-1,1)
-    
-    x2 = x1 + width
-    y2 = y1 
-    
-    x3 = x1
-    y3 = y1 + height
-    
-    x4 = bboxes[:,2].reshape(-1,1)
-    y4 = bboxes[:,3].reshape(-1,1)
-    
-    corners = np.hstack((x1,y1,x2,y2,x3,y3,x4,y4))
-    
-    return corners
-
-
-def rotate_box(corners,angle,  cx, cy, h, w):
-    
-    """Rotate the bounding box.
-    
-    
-    Parameters
-    ----------
-    
-    corners : numpy.ndarray
-        Numpy array of shape `N x 8` containing N bounding boxes each described by their 
-        corner co-ordinates `x1 y1 x2 y2 x3 y3 x4 y4`
-    
-    angle : float
-        angle by which the image is to be rotated
-        
-    cx : int
-        x coordinate of the center of image (about which the box will be rotated)
-        
-    cy : int
-        y coordinate of the center of image (about which the box will be rotated)
-        
-    h : int 
-        height of the image
-        
-    w : int 
-        width of the image
-    
-    Returns
-    -------
-    
-    numpy.ndarray
-        Numpy array of shape `N x 8` containing N rotated bounding boxes each described by their 
-        corner co-ordinates `x1 y1 x2 y2 x3 y3 x4 y4`
-    """
-
-    corners = corners.reshape(-1,2)
-    corners = np.hstack((corners, np.ones((corners.shape[0],1), dtype = type(corners[0][0]))))
-    
-    M = cv2.getRotationMatrix2D((cx, cy), angle, 1.0)
-    
-    
-    cos = np.abs(M[0, 0])
-    sin = np.abs(M[0, 1])
-    
-    nW = int((h * sin) + (w * cos))
-    nH = int((h * cos) + (w * sin))
-    # adjust the rotation matrix to take into account translation
-    M[0, 2] += (nW / 2) - cx
-    M[1, 2] += (nH / 2) - cy
-    # Prepare the vector to be transformed
-    calculated = np.dot(M,corners.T).T
-    
-    calculated = calculated.reshape(-1,8)
-    
-    return calculated
-
-def get_enclosing_box(corners):
-    """Get an enclosing box for ratated corners of a bounding box
-    
-    Parameters
-    ----------
-    
-    corners : numpy.ndarray
-        Numpy array of shape `N x 8` containing N bounding boxes each described by their 
-        corner co-ordinates `x1 y1 x2 y2 x3 y3 x4 y4`  
-    
-    Returns 
-    -------
-    
-    numpy.ndarray
-        Numpy array containing enclosing bounding boxes of shape `N X 4` where N is the 
-        number of bounding boxes and the bounding boxes are represented in the
-        format `x1 y1 x2 y2`
-        
-    """
-    x_ = corners[:,[0,2,4,6]]
-    y_ = corners[:,[1,3,5,7]]
-    
-    xmin = np.min(x_,1).reshape(-1,1)
-    ymin = np.min(y_,1).reshape(-1,1)
-    xmax = np.max(x_,1).reshape(-1,1)
-    ymax = np.max(y_,1).reshape(-1,1)
-    
-    final = np.hstack((xmin, ymin, xmax, ymax,corners[:,8:]))
-    
-    return final
-
 def spinim(self, img, bboxes):
 
     angle = random.uniform(*self.angle)