From 7c798c88484b16f999e41aae03e86cb6d0f8f6ad Mon Sep 17 00:00:00 2001 From: Weiwei Duan Date: Fri, 4 Dec 2020 10:18:28 -0800 Subject: [PATCH] Add files via upload --- pipeline/evaluation/buffer.py | 38 ++++++++ pipeline/evaluation/skeletonize.py | 33 +++++++ pipeline/evaluation/utils_shp.py | 147 +++++++++++++++++++++++++++++ pipeline/evaluation/workspace.py | 63 +++++++++++++ 4 files changed, 281 insertions(+) create mode 100644 pipeline/evaluation/buffer.py create mode 100644 pipeline/evaluation/skeletonize.py create mode 100644 pipeline/evaluation/utils_shp.py create mode 100644 pipeline/evaluation/workspace.py diff --git a/pipeline/evaluation/buffer.py b/pipeline/evaluation/buffer.py new file mode 100644 index 0000000..78da17e --- /dev/null +++ b/pipeline/evaluation/buffer.py @@ -0,0 +1,38 @@ +import numpy as np +import cv2 + +def buffer(pred_ske, buffer_size=1): + if np.max(pred_ske)==0: + print('blank img') + return + if np.max(pred_ske)!=1: + pred_ske= pred_ske/np.max(pred_ske) + + for i in range(buffer_size): + nonzeros_idx = np.where(pred_ske!=0) + xs, ys = nonzeros_idx + # west + d1 = (xs-1, ys) + # east + d2 = (xs+1, ys) + # north + d3 = (xs, ys-1) + # south + d4 = (xs, ys+1) + # northwest + d5 = (xs-1, ys-1) + # northeast + d6 = (xs+1, ys-1) + # southwest + d7 = (xs-1, ys+1) + # southeast + d8 = (xs+1, ys+1) + pred_ske[d1]=1 + pred_ske[d2]=1 + pred_ske[d3]=1 + pred_ske[d4]=1 + pred_ske[d5]=1 + pred_ske[d6]=1 + pred_ske[d7]=1 + pred_ske[d8]=1 + return pred_ske diff --git a/pipeline/evaluation/skeletonize.py b/pipeline/evaluation/skeletonize.py new file mode 100644 index 0000000..2d0831a --- /dev/null +++ b/pipeline/evaluation/skeletonize.py @@ -0,0 +1,33 @@ +from skimage.morphology import skeletonize +import numpy as np +import cv2 +import os + +def preprocess(img): + img = (img == 255).astype(np.bool) + return img + +def make_skeleton(root, fix_borders=True, debug=False): + replicate = 5 + clip = 2 + rec = replicate + clip + # open and skeletonize + img = cv2.imread(root, cv2.IMREAD_GRAYSCALE) + print(img.shape) + + if fix_borders: + img = cv2.copyMakeBorder(img, replicate, replicate, replicate, replicate, cv2.BORDER_REPLICATE) + img_copy = None + if debug: + if fix_borders: + img_copy = np.copy(img[replicate:-replicate,replicate:-replicate]) + else: + img_copy = np.copy(img) + img = preprocess(img) + if not np.any(img): + return None, None + ske = skeletonize(img).astype(np.uint16) + if fix_borders: + ske = ske[rec:-rec, rec:-rec] + ske = cv2.copyMakeBorder(ske, clip, clip, clip, clip, cv2.BORDER_CONSTANT, value=0) + return img_copy, ske diff --git a/pipeline/evaluation/utils_shp.py b/pipeline/evaluation/utils_shp.py new file mode 100644 index 0000000..f3072fb --- /dev/null +++ b/pipeline/evaluation/utils_shp.py @@ -0,0 +1,147 @@ +import subprocess +import os, cv2 +import xml.etree.ElementTree as ET +from osgeo import gdal, gdalconst, ogr, osr +import zipfile +import numpy as np +from gdalconst import * + +#full path to gdal executables> +gdalsrsinfo = r'C:\OSGeo4W64\bin\gdalsrsinfo.exe' +ogr2ogr = r'C:\OSGeo4W64\bin\ogr2ogr.exe' +gdal_translate = r'C:\OSGeo4W64\bin\gdal_translate.exe' +#------------------------------------------------------------------------------ + +def unzip(zip_path, unzip_path): + with zipfile.ZipFile(zip_path ,"r") as zip_ref: + zip_ref.extractall(unzip_path) + temp = unzip_path.split('\\') + map_name = temp[-1].replace('_tif', '.tif') + return os.path.join(unzip_path,map_name) + + +def reproject(input_vec, in_crs, out_crs, ogr2ogr): + vector_proj = input_vec.replace('.shp','_proj.shp') + call = ogr2ogr+' -t_srs "'+out_crs+'" -s_srs "'+in_crs+'" "'+vector_proj+'" "'+input_vec+'"' + print (call) + response=subprocess.check_output(call, shell=False) + print(response) + return vector_proj + +def get_bbox(xml_file): + bbox = [] #west, east, north, south + tree = ET.parse(xml_file) + root = tree.getroot() + for child in root: + if child.tag == 'idinfo': + for cchild in child: + if cchild.tag =='spdom': + for ccchild in cchild: # ccchild is bounding + for coor in ccchild: + # print coor.text + bbox.append(float(coor.text)) + # xmin,xmax,ymin,ymax = bbox + return bbox + +# convert bbox points to goecoordinate of maps +def point_convertor(x, y, input_epsg, out_srs): + point = ogr.Geometry(ogr.wkbPoint) + point.AddPoint(x, y) + # print(point) + # create coordinate transformation + inSpatialRef = osr.SpatialReference() + inSpatialRef.ImportFromEPSG(input_epsg) + coordTransform = osr.CoordinateTransformation(inSpatialRef, out_srs) + # transform point + point.Transform(coordTransform) + return [point.GetX(), point.GetY()] + +def clip(input_vec,xmin,ymin,xmax,ymax,ogr2ogr): + vector_clip = input_vec.replace('.shp','_clip.shp') + # call = '%s -dim 2 -clipsrc %f %f %f %f -nlt POLYGON %s %s ' % (ogr2ogr, xmin,ymin,xmax,ymax, vector_clip, input_vec) + call = '%s -dim 2 -clipsrc %f %f %f %f %s %s ' % (ogr2ogr, xmin,ymin,xmax,ymax, vector_clip, input_vec) + print (call) + response=subprocess.check_output(call, shell=False) + print(response) + return vector_clip + + +def create_filtered_shapefile(value, filter_field, in_shapefile, out_shapefile): + driver = ogr.GetDriverByName("ESRI Shapefile") + dataSource = driver.Open(in_shapefile, 0) + input_layer = dataSource.GetLayer() + query_str = '{} = {}'.format(filter_field, value) + print(query_str) + err = input_layer.SetAttributeFilter(query_str) + print(err) + # Copy Filtered Layer and Output File + out_ds = driver.CreateDataSource(out_shapefile) + out_layer = out_ds.CopyLayer(input_layer, str(value)) + del input_layer, out_layer, out_ds + return out_shapefile + +def vector2raster(input_vec, map_tiff_geo): +# open data + raster_fn = input_vec.replace('.shp', '.tif') + raster = gdal.Open(map_tiff_geo) + shp = ogr.Open(input_vec) + lyr = shp.GetLayer() +# # Get raster georeference info + transform = raster.GetGeoTransform() + # print(transform) +# # Create memory target raster + x_res, y_res = raster.RasterXSize, raster.RasterYSize + target_ds = gdal.GetDriverByName('GTiff').Create(raster_fn, x_res, y_res, 1, gdal.GDT_Byte) + target_ds.SetGeoTransform(transform) + raster_srs = osr.SpatialReference() + raster_srs.ImportFromWkt(raster.GetProjectionRef()) + target_ds.SetProjection(raster_srs.ExportToWkt()) + # Rasterize + err = gdal.RasterizeLayer(target_ds, [1], lyr,burn_values=[255]) + if err != 0: + print(err) + del target_ds + return raster_fn + +def tif2png(in_tif): + png_fn = in_tif.replace('.tif', '.png') + call = 'gdal_translate -of PNG -ot Byte "'+in_tif+'" "'+png_fn+'"' + print(call) + response=subprocess.check_output(call, shell=True) + print(response) + xml_fn = png_fn + '.aux.xml' + os.remove(xml_fn) + return png_fn + +def dilated(png_fn): + dilated_fn = png_fn.replace('.png', '_dilated.png') + img = cv2.imread(png_fn,0) + kernel = np.ones((2,2),np.uint8) + dilation = cv2.dilate(img,kernel,iterations = 1) + cv2.imwrite(dilated_fn, dilation) + return dilated_fn + +def geo2img_coor(x, y, path): + dataset = gdal.Open( path, GA_ReadOnly ) + adfGeoTransform = dataset.GetGeoTransform() + dfGeoX=float(x) + dfGeoY =float(y) + det = adfGeoTransform[1] * adfGeoTransform[5] - adfGeoTransform[2] *adfGeoTransform[4] + # X = (int)(((dfGeoX - adfGeoTransform[0]) / adfGeoTransform[1])) + # Y = (int)(((dfGeoY - adfGeoTransform[3]) / adfGeoTransform[5])) + X = ((dfGeoX - adfGeoTransform[0]) * adfGeoTransform[5] - (dfGeoY - + adfGeoTransform[3]) * adfGeoTransform[2]) / det + Y = ((dfGeoY - adfGeoTransform[3]) * adfGeoTransform[1] - (dfGeoX - + adfGeoTransform[0]) * adfGeoTransform[4]) / det + return [int(Y),int(X)] + +def points_generator(raster, save_path, n_points=400): + img = cv2.imread(raster, 0) + index_x, index_y = np.where(img==255) + print(index_x.shape) + rand = np.random.randint(0, index_x.shape[0], n_points) + points = np.array([index_x[rand], index_y[rand]]) + points = np.swapaxes(points,0,1) + print(points.shape) + np.savetxt(save_path, points, fmt='%d', delimiter=',') + return save_path diff --git a/pipeline/evaluation/workspace.py b/pipeline/evaluation/workspace.py new file mode 100644 index 0000000..548aa11 --- /dev/null +++ b/pipeline/evaluation/workspace.py @@ -0,0 +1,63 @@ +import numpy as np +import os, cv2 +import copy +from skeletonize import make_skeleton +from buffer import buffer +import utils_shp + +''' +inputs for the evaluation are +ground truth shape file +segmentation results +''' + +data_dir = 'C:\Users\weiweiduan\Documents\Map_proj_data\CA\CA_Bray_100414_2001_24000_geo_tif' +location_name = 'CA_Bray_100414_2001_24000_geo' +gt_name = 'CA_Bray_railroads_gt_buffer0.shp' +pred_name = 'CA_Bray_100414_2001_24000_geo_pointrend_unet_pred.png' +buffer_size = 2 + +pred_path = os.path.join(data_dir, 'res', pred_name) +gt_shp_path = os.path.join(data_dir, 'ground_truth', gt_name) +map_tif_path = os.path.join(data_dir, location_name+'.tif') + +# Step 1: rasterize the ground truth shapefile +gt_tif_path = utils_shp.vector2raster(gt_path, map_tif_path) +gt_png_path = utils_shp.tif2png(gt_tif_path) + +# Step 2: load the bounding box coordinates to remove the white borders of the map +bbox_file = os.path.join(data_dir, 'bbox.txt') +points = np.loadtxt(bbox_file, dtype='int32', delimiter=',') +start_point, end_point = points[0], points[1] +print(start_point, end_point) + +# Step 3: skeletonize the segmentation results +img_copy, pred_ske = make_skeleton(pred_path) +pred_ske = pred_ske.astype('uint') + +# Step 4: buffer the segmentation results and ground truth +pred_buffer = buffer(pred_ske, buffer_size=buffer_size) +gt_map = cv2.imread(gt_png_path,0) / 255 +gt_buffer = buffer(gt_map, buffer_size=buffer_size) + +# Step 5 calculate correctness and completeness +overlap_map = gt_buffer * pred_ske +fp_map = pred_ske - overlap_map + + +tp = np.count_nonzero(overlap_map) +fp = np.count_nonzero(fp_map) + +correctness = tp / (tp+fp) +print('correctness = ', correctness) + +overlap_comp_map = gt_map * pred_buffer + +fn_map = gt_map - overlap_comp_map + +tp_comp = np.count_nonzero(overlap_comp_map) +fn = np.count_nonzero(fn_map) + +completeness = tp_comp / (tp_comp+fn) +print('completeness = ', completeness) +