Merge branch 'revert-c29534a9' into 'master'

Revert "Merge branch 'master' of"

See merge request !1

Common.py

+'''
+Created on 15.03.2017
+
+@author: luroth
+'''
+
+import os
+import math
+import numpy as np
+from matplotlib import path
+from osgeo import gdal
+
+DELTA_x = 2693000
+DELTA_y = 1256000
+
+def rotation_matrix(axis, theta):
+    axis = np.asarray(axis)
+    axis = axis/math.sqrt(np.dot(axis, axis))
+    a = math.cos(theta/2.0)
+    b, c, d = -axis*math.sin(theta/2.0)
+    aa, bb, cc, dd = a*a, b*b, c*c, d*d
+    bc, ad, ac, ab, bd, cd = b*c, a*d, a*c, a*b, b*d, c*d
+    return np.array([[aa+bb-cc-dd, 2*(bc+ad), 2*(bd-ac)],
+                     [2*(bc-ad), aa+cc-bb-dd, 2*(cd+ab)],
+                     [2*(bd+ac), 2*(cd-ab), aa+dd-bb-cc]])
+    
+def find_mins_maxs(mesh):
+    minx = maxx = miny = maxy = minz = maxz = None
+    for p in mesh.points:
+        # p contains (x, y, z)
+        if minx is None:
+            minx = p[stl.Dimension.X]
+            maxx = p[stl.Dimension.X]
+            miny = p[stl.Dimension.Y]
+            maxy = p[stl.Dimension.Y]
+            minz = p[stl.Dimension.Z]
+            maxz = p[stl.Dimension.Z]
+        else:
+            maxx = max(p[stl.Dimension.X], maxx)
+            minx = min(p[stl.Dimension.X], minx)
+            maxy = max(p[stl.Dimension.Y], maxy)
+            miny = min(p[stl.Dimension.Y], miny)
+            maxz = max(p[stl.Dimension.Z], maxz)
+            minz = min(p[stl.Dimension.Z], minz)
+    return minx, maxx, miny, maxy, minz, maxz
+    
+def get_polygon_normale(points):
+    a = np.array(points[0])
+    b = np.array(points[1])
+    c = np.array(points[2])
+
+    v1 = a-b
+    v2 = a-c
+    normale = np.cross(v1, v2)
+    
+    return normale
+
+def getWKT_PRJ (epsg_code):
+    try:
+        file_dir = os.path.dirname(os.path.realpath(__file__))
+    except:
+        file_dir = '/home/luroth/git_repos/MSc_Lukas_Roth/'
+    
+    with open(os.path.join(file_dir, 'EPSG_files', 'EPSG' + str(epsg_code) + '.txt'), 'r') as file:
+        remove_spaces = file.read()
+        remove_spaces = remove_spaces.replace(' ', '')
+        output = remove_spaces.replace("\n", "")
+    return output
+
+def point_in_polygon(x,y,poly,precision):        
+    path_ = path.Path(poly)
+    return path_.contains_point((x, y), radius=precision)
+
+from shapely.geometry import Point
+from shapely.geometry.polygon import Polygon
+
+def point_in_polygon2(x,y,poly):
+    point = Point(x, y)
+    polygon = Polygon(poly)
+    return polygon.contains(point)
+
+def point_in_poly(x,y,poly):
+    n = len(poly)
+    inside = False
+
+    p1x,p1y = poly[0]
+    for i in range(n+1):
+        p2x,p2y = poly[i % n]
+        if y > min(p1y,p2y):
+            if y <= max(p1y,p2y):
+                if x <= max(p1x,p2x):
+                    if p1y != p2y:
+                        xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
+                    if p1x == p2x or x <= xints:
+                        inside = not inside
+        p1x,p1y = p2x,p2y
+
+    return inside
+
+def create_GeoTiff(file_path, array, data_type):
+    
+    driver = gdal.GetDriverByName('GTiff')
+    file = file_path + '.tif'
+    
+    ds = driver.Create(file, array.shape[1], array.shape[0], 1, data_type)
+
+    band_1 = ds.GetRasterBand(1)
+    band_1.WriteArray(array)
+    band_1.SetNoDataValue(0)
+    
+    ds = None
+    
+    return file_path

GenerateSamplingMasks.py

@@ -25,13 +25,13 @@ from progressbar import ProgressBar
 
 from Common import *
 
-def mask_worker(work_queue, result_queue, metadata_folder, output_directory, lot, sub_plots):
+def mask_worker(work_queue, result_queue, metadata_folder, output_directory, sub_plots):
     for job in iter(work_queue.get, 'STOP'):
         PhotoID = job['PhotoID']
-        create_mask(metadata_folder, output_directory, lot, sub_plots, PhotoID)
+        create_mask(metadata_folder, output_directory, sub_plots, PhotoID)
         result_queue.put(PhotoID)
 
-def create_mask(metadata_folder, output_directory, lot, sub_plots, PhotoID):
+def create_mask(metadata_folder, output_directory, sub_plots, PhotoID):
 
     photo_meta_x_ds = gdal.Open(os.path.join(metadata_folder, PhotoID + '_metadata_coord_x.tif'))
     photo_meta_x = np.array(photo_meta_x_ds.GetRasterBand(1).ReadAsArray(), dtype=np.float64)
@@ -47,15 +47,6 @@ def create_mask(metadata_folder, output_directory, lot, sub_plots, PhotoID):
 
     photo_meta_xy = np.stack((photo_meta_x, photo_meta_y), axis=2)
     photo_meta_xy = photo_meta_xy.reshape(-1,2) 
-    
-    lot_path = path.Path(lot)
-    mask_lot = lot_path.contains_points(photo_meta_xy, radius=0.02)
-    if not np.any(mask_lot):
-        return
-    mask_lot = mask_lot.reshape(photo_meta_y.shape)
-    #gdalnumeric.SaveArray(1*mask_lot, os.path.join(output_directory, PhotoID + '_mask_lot.tif'), format='GTiff')
-    create_1band_GeoTiff(os.path.join(output_directory, PhotoID + '_mask_lot.tif'), 1*mask_lot, gdal.GDT_Byte)
-            
         
     mask_subplots = np.empty(photo_meta_xy.shape[0:1], dtype=int)
     
@@ -76,8 +67,8 @@ def create_mask(metadata_folder, output_directory, lot, sub_plots, PhotoID):
     
     if found:
         mask_subplots = mask_subplots.reshape(photo_meta_y.shape)
-        create_1band_GeoTiff(os.path.join(output_directory, PhotoID + '_mask_subsamples.tif'), mask_subplots, gdal.GDT_Int16)
-        #gdalnumeric.SaveArray(mask_subplots, os.path.join(output_directory, PhotoID + '_mask_subsamples.tif'), format='GTiff')
+
+        gdalnumeric.SaveArray(mask_subplots, os.path.join(output_directory, PhotoID + '_mask_subsamples.tif'), format='GTiff')
 
 def main(argv):
     camera_position_file = ''
@@ -105,10 +96,6 @@ def main(argv):
             sub_plots_file = arg                
     
     print('Metdadata folder is', metadata_folder)
-    
-    with open(lot_dim, 'r') as f:
-        lot_polygons = geojson.load(f)
-        lot = lot_polygons['features'][0]['geometry']['coordinates'][0]
             
     with open(sub_plots_file, 'r') as f:
         sub_plots_tmp = geojson.load(f)
@@ -139,7 +126,7 @@ def main(argv):
         progress = ProgressBar(min_value=0, max_value=max_jobs)
             
         for w in range(multiprocessing.cpu_count()-2):
-            p = Process(target=mask_worker, args=(jobs, results, metadata_folder, output_directory, lot, sub_plots))
+            p = Process(target=mask_worker, args=(jobs, results, metadata_folder, output_directory, sub_plots))
             p.daemon = True
             p.start()
             processes.append(p)

Segmentation.py

+'''
+Created on 16.03.2017
+
+@author: luroth
+'''
+
+#from pylab import plot,show
+from numpy import vstack,array
+from numpy.random import rand
+import numpy.ma as ma
+import numpy as np
+from scipy.cluster.vq import kmeans,vq
+from scipy import misc
+#import sklearn.cluster
+import matplotlib.pyplot as plt
+from osgeo import gdalnumeric
+import os
+import sys
+import getopt
+import csv
+from scipy import ndimage
+from Common import create_GeoTiff
+import gdal
+import colorsys
+
+def main(argv):
+    camera_position_file = ''
+    metadata_folder = ''
+    output_directory = ''
+    
+    
+    try:
+        opts, args = getopt.getopt(argv,"hc:m:p:s:o:")
+    except getopt.GetoptError:
+        print('test.py -c <camera position file> -m <metadata folder> -o <Output directory>')
+        sys.exit(2)
+    for opt, arg in opts:
+        if opt == '-h':
+            print('test.py -c <camera position file> -m <metadata folder> -o <Output directory>')
+            sys.exit()
+        elif opt == "-c":
+            camera_position_file = arg
+        elif opt == "-m":
+            metadata_folder = arg    
+        elif opt == "-o":
+            output_directory = arg                   
+    
+    print('Metdadata folder is', metadata_folder)
+
+
+    with open(camera_position_file, 'r') as csvfile:
+        spamreader = csv.reader(csvfile, delimiter='\t', quotechar='|')
+        next(spamreader, None)
+        next(spamreader, None)
+ 
+        for row in spamreader:  
+            PhotoID = os.path.splitext(row[0])[0]
+            
+            image = misc.imread(os.path.join(metadata_folder, '../../RGB', PhotoID + '.JPG'))
+            
+            if os.path.isfile(os.path.join(metadata_folder, PhotoID + '_mask_lot.tif')):
+                image_mask_all =  misc.imread(os.path.join(metadata_folder, PhotoID + '_mask_lot.tif'))
+                image_mask_soil =  misc.imread(os.path.join(metadata_folder, PhotoID + '_mask_soil.tif'))
+            
+                origin_shape = image.shape
+                
+                image = image.reshape(-1,3)
+                image_mask_all = image_mask_all.reshape(-1)
+                image_mask_soil = image_mask_soil.reshape(-1)
+                
+                image[image_mask_all==0,:] = 0
+                
+                # Colors
+                R = image[:,0]/255.0
+                G = image[:,1]/255.0
+                B = image[:,2]/255.0
+
+                # ExG
+                r= np.divide(R, (np.mean(R)+np.mean(G)+np.mean(B)), dtype=float)
+                g= np.divide(G, (np.mean(R)+np.mean(G)+np.mean(B)), dtype=float)
+                b= np.divide(B, (np.mean(R)+np.mean(G)+np.mean(B)), dtype=float)
+                
+                ExG = 2.0* g - r - b
+                ExR = 1.4 * r - b
+                
+                misc.imsave(os.path.join(output_directory, PhotoID + '_ExG.tif'), ExG.reshape(origin_shape[0], origin_shape[1]))
+                misc.imsave(os.path.join(output_directory, PhotoID + '_ExR.tif'), ExR.reshape(origin_shape[0], origin_shape[1]))
+                
+                ExG_ExR = ExG - ExR
+                
+                plants = np.array(ExG_ExR, dtype=float)
+                plants = plants.reshape(origin_shape[0], origin_shape[1])
+                
+                #gdalnumeric.SaveArray(1*(plants.reshape(origin_shape[0], origin_shape[1])), os.path.join(output_directory, PhotoID + '_segmented_ExG_ExR.tif'), format='GTiff')
+                create_GeoTiff(os.path.join(output_directory, PhotoID + '_segmented_ExG_ExR'), plants, gdal.GDT_Float32)
+                
+                #ExG_soil = ExG[image_mask_soil==1]
+                #print(ExG_soil.size)
+                #if ExG_soil.size ==0:
+                #    continue
+                
+                #ExG_threshold = np.percentile(ExG_soil, 99.95)
+                
+                #soil = np.array(ExG > ExG_threshold)
+                
+                #gdalnumeric.SaveArray(1*(soil.reshape(origin_shape[0], origin_shape[1])), os.path.join(output_directory, PhotoID + '_segmented_ExG.tif'), format='GTiff')
+                #misc.imsave(os.path.join(output_directory, PhotoID + '_segmented_ExG.tif'), soil.reshape(origin_shape[0], origin_shape[1]))
+                
+                
+                print(PhotoID, "done")
+
+if __name__ == "__main__":
+    main(sys.argv[1:])
+    
+    print("\n============\nAll images segmented, END\n")

ZonalStat.py

+'''
+Created on 16.03.2017
+
+@author: luroth
+'''
+
+from pylab import plot,show
+from numpy import vstack,array
+from numpy.random import rand
+import numpy.ma as ma
+import numpy as np
+from scipy.cluster.vq import kmeans,vq
+from scipy import misc
+import sklearn.cluster
+import multiprocessing
+import matplotlib.pyplot as plt
+from osgeo import gdalnumeric
+import os
+import sys
+import getopt
+import csv
+from scipy import ndimage
+import geojson
+import json
+from Common import *
+from osgeo import ogr, gdal
+import pandas as pd
+from multiprocessing import Process, Queue, Manager
+import time
+from progressbar import ProgressBar
+
+def zonal_worker(work_queue, result_queue, metadata_folder, sub_plots):
+    for job in iter(work_queue.get, 'STOP'):
+        PhotoID = job['PhotoID']
+        result = zonal_stat(metadata_folder, sub_plots, PhotoID)
+        result_queue.put(result)
+
+def zonal_stat(metadata_folder, sub_plots, PhotoID):
+    samples = []
+    
+    subplots_ds = gdal.Open(os.path.join(metadata_folder, PhotoID + '_mask_subsamples.tif'))
+    subplots = np.array(subplots_ds.GetRasterBand(1).ReadAsArray(), dtype=int)
+    origin_shape = subplots.shape
+    subplots = subplots.reshape(-1)
+    
+    photo_meta_x_ds = gdal.Open(os.path.join(metadata_folder, PhotoID + '_metadata_coord_x.tif'))
+    photo_meta_x = np.array(photo_meta_x_ds.GetRasterBand(1).ReadAsArray(), dtype=np.float64)
+    photo_meta_x += DELTA_x
+    
+    photo_meta_y_ds = gdal.Open(os.path.join(metadata_folder, PhotoID + '_metadata_coord_y.tif'))
+    photo_meta_y = np.array(photo_meta_y_ds.GetRasterBand(1).ReadAsArray(), dtype=np.float64)
+    photo_meta_y += DELTA_y
+
+    photo_meta_xy = np.stack((photo_meta_x, photo_meta_y), axis=2)
+    photo_meta_xy = photo_meta_xy.reshape(-1,2)
+    
+    photo_meta_x = photo_meta_x.reshape(-1)
+    photo_meta_y = photo_meta_y.reshape(-1)
+    
+    segmented  =  misc.imread(os.path.join(metadata_folder, PhotoID + '_segmented_ExG.tif'))
+    segmented = segmented.reshape(-1)
+    
+    photo_meta_z  =  misc.imread(os.path.join(metadata_folder, PhotoID + '_metadata_coord_z.tif'))
+    photo_meta_z = photo_meta_z.reshape(-1)
+    
+    photo_meta_angle_z  =  misc.imread(os.path.join(metadata_folder, PhotoID + '_metadata_angle_z.tif'))
+    photo_meta_angle_z = photo_meta_angle_z.reshape(-1)
+    photo_meta_angle_xy  =  misc.imread(os.path.join(metadata_folder, PhotoID + '_metadata_angle_xy.tif'))
+    photo_meta_angle_xy = photo_meta_angle_xy.reshape(-1)
+
+    for sub_plot in sub_plots:
+        plot_id = sub_plot['properties']['plot_number']
+        subsample = sub_plot['properties']['subsample']
+        
+        #print("calculating plot id", plot_id, ', subsample', subsample)
+        
+        mask_subplot = subplots == (plot_id*100+subsample)
+        
+        subplot_segmented = segmented[mask_subplot]
+        
+        if subplot_segmented.size > 0:
+            subplot_mean_cc = np.mean(subplot_segmented)
+            subplot_mean_x = np.mean(photo_meta_x[mask_subplot])
+            subplot_mean_y = np.mean(photo_meta_y[mask_subplot])
+            subplot_mean_z = np.mean(photo_meta_z[mask_subplot])
+            subplot_mean_angle_z = np.mean(photo_meta_angle_z[mask_subplot])
+            subplot_mean_angle_xy = np.mean(photo_meta_angle_xy[mask_subplot])
+            
+            image_x, image_y = np.where(mask_subplot.reshape(origin_shape))
+            image_x = np.mean(image_x)
+            image_y = np.mean(image_y)
+        
+            #print(subplot_mean_cc, end="", flush=True)
+            samples.append({"plot_id":plot_id, 'subsample':subsample, 'CC':subplot_mean_cc,
+                           'glob_X':subplot_mean_x,
+                           'glob_Y':subplot_mean_y,
+                           'glob_Z':subplot_mean_z,
+                           'angle_z':subplot_mean_angle_z,
+                           'angle_xy':subplot_mean_angle_xy,
+                           'image':PhotoID,
+                           'image_x':image_x,
+                           'image_y':image_y})
+        else:
+            pass
+    
+    return samples
+
+
+def main(argv):
+    camera_position_file = ''
+    metadata_folder = ''
+    output_directory = ''
+    sub_plots_file = ''
+    
+    try:
+        opts, args = getopt.getopt(argv,"hc:m:p:s:o:")
+    except getopt.GetoptError:
+        print('test.py -c <camera position file> -m <metadata folder> -o <Output directory>')
+        sys.exit(2)
+    for opt, arg in opts:
+        if opt == '-h':
+            print('test.py -c <camera position file> -m <metadata folder> -o <Output directory>')
+            sys.exit()
+        elif opt == "-c":
+            camera_position_file = arg
+        elif opt == "-m":
+            metadata_folder = arg    
+        elif opt == "-o":
+            output_directory = arg      
+        elif opt == "-s":
+            sub_plots_file = arg              
+    
+    print('Metdadata folder is', metadata_folder)
+    
+    all_cc = []
+    
+    sub_plots_file = r'E:\UAV\_Workspace_\2016_ethz_eschikon_FIP_50m_general\design_subsamples.geojson'
+    with open(sub_plots_file, 'r') as f:
+        sub_plots_tmp = geojson.load(f)
+        sub_plots = sub_plots_tmp['features']
+        
+    with open(camera_position_file, 'r') as csvfile:
+        spamreader = csv.reader(csvfile, delimiter='\t', quotechar='|')
+        next(spamreader, None)
+        next(spamreader, None)
+   
+        m = Manager()
+        jobs = m.Queue()
+        results = m.Queue()
+        processes = []
+        
+        max_jobs = 0
+        count = 0
+        for row in spamreader:  
+            PhotoID = os.path.splitext(row[0])[0]
+            
+            if (os.path.isfile(os.path.join(metadata_folder, PhotoID + '_mask_subsamples.tif')) and 
+                os.path.isfile(os.path.join(metadata_folder, PhotoID + '_segmented_ExG.tif'))):
+                
+                job = dict()
+                job['PhotoID'] = PhotoID
+                jobs.put(job)
+                 
+                max_jobs += 1
+            
+        progress = ProgressBar(min_value=0, max_value=max_jobs)
+            
+        for w in range(multiprocessing.cpu_count()-2):
+            p = Process(target=zonal_worker, args=(jobs, results, metadata_folder, sub_plots))
+            p.daemon = True
+            p.start()
+            processes.append(p)
+            jobs.put('STOP')
+        
+        print("jobs all started")
+        
+        samples = []
+        progress.update(0)
+        while count < (max_jobs):
+            samples_ = results.get()
+            samples.extend(samples_)
+            progress.update(count)
+            count+=1
+            
+        progress.finish()
+        
+        for p in processes:
+            p.join()
+        
+        pd_data = pd.DataFrame(samples)
+        pd_data.to_csv(os.path.join(metadata_folder, "CC.csv"))
+                
+if __name__ == "__main__":
+    main(sys.argv[1:])
+    
+    print("\n============\nAll images segmented, END\n")