diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..52b2c0a --- /dev/null +++ b/.gitignore @@ -0,0 +1,11 @@ +tohme_crops +*.pyc +old +*.swp +venv +crops/* +*.log +crop.log.bak.txt +.idea +*.csv +cropsbak diff --git a/CropRunner.py b/CropRunner.py new file mode 100644 index 0000000..362ae2d --- /dev/null +++ b/CropRunner.py @@ -0,0 +1,390 @@ +""" +** Crop Extractor for Project Sidewalk ** + +Given label metadata from the Project Sidewalk database, this script will +extract JPEG crops of the features that have been labeled. The required metadata +may be obtained by running the SQL query in "samples/getFullLabelList.sql" on the +Sidewalk database, and exporting the results in CSV format. You must supply the +path to the CSV file containing this data below. You can find an example of what +this file should look like in "samples/labeldata.csv". + +Additionally, you should have downloaded original panorama +images from Street View using DownloadRunner.py. You will need to supply the +path to the folder containing these files. + +""" + +# ***************************************** +# Update paths below * +# ***************************************** + +# Path to CSV data from database +csv_export_path = "labeldata.csv" +# Path to panoramas downloaded using DownloadRunner.py +gsv_pano_path = "/mnt/umiacs/Panoramas/scrapes4" +# Path to location for saving the crops +destination_path = "/home/anthony/Downloads/crops" + +# Mark the center of the crop? +mark_center = True + +import csv +import GSVImage +import fnmatch + +import logging + +logging.basicConfig(filename='crop.log', level=logging.DEBUG) +from utilities import * + +from PIL import Image, ImageDraw + +try: + from xml.etree import cElementTree as ET +except ImportError, e: + from xml.etree import ElementTree as ET + + +def extract_panoyawdeg(path_to_metadata_xml): + pano = {} + pano_xml = open(path_to_metadata_xml, 'rb') + tree = ET.parse(pano_xml) + root = tree.getroot() + for child in root: + if child.tag == 'projection_properties': + pano[child.tag] = child.attrib + print(pano['projection_properties']['pano_yaw_deg']) + + return pano['projection_properties']['pano_yaw_deg'] + + +def crop_box_helper(path_to_scrapes, path_to_labeldata_csv, target_label_type=2): + target = open("cropboxes.log", 'w') + target.truncate(); + + for root, dirnames, filenames in os.walk(path_to_scrapes): + for filename in fnmatch.filter(filenames, '*.depth.txt'): + depth_location = os.path.join(root, filename) + + pano_id = filename[:-10] + # Search through CSV for labels that exist in this panorama + csv_file = open(path_to_labeldata_csv) + csv_f = csv.reader(csv_file) + num_matching_labels = 0 + + path_to_xml = os.path.join(root, pano_id + ".xml") + image_name = os.path.join(root, pano_id + ".jpg") + pano_im = Image.open(image_name) + + for row in csv_f: + # Skip the header row + if row[0] == "gsv_panorama_id": + continue + csv_pano_id = row[0] + sv_image_x = float(row[1]) + + sv_image_y = float(row[2]) + label_type = int(row[3]) + photographer_heading = float(row[4]) + heading = float(row[5]) + label_id = int(row[7]) + + pano_yaw_deg = 180 - photographer_heading + if csv_pano_id == pano_id and label_type == target_label_type: + num_matching_labels += 1 + im_width = GSVImage.GSVImage.im_width + im_height = GSVImage.GSVImage.im_height + + draw = ImageDraw.Draw(pano_im) + # sv_image_x = sv_image_x - 100 + x = ((float(pano_yaw_deg) / 360) * im_width + sv_image_x) % im_width + y = im_height / 2 - sv_image_y + r = 50 + draw.ellipse((x - r, y - r, x + r, y + r), fill=128) + if num_matching_labels == 0: + print("Skipping image because no labels of interest were found.") + continue + else: + print("Found " + str(num_matching_labels) + " labels of interest in this image.") + + figure() + im = imshow(pano_im) + fig = gcf() + ax = gca() + + class EventHandler: + def __init__(self): + self.top_left_x = None + self.top_left_y = None + self.bottom_right_x = None + self.bottom_right_y = None + fig.canvas.mpl_connect('button_press_event', self.onpress) + + def onpress(self, event): + + if event.inaxes != ax: + return + xi, yi = (int(round(n)) for n in (event.xdata, event.ydata)) + if self.top_left_x is None and self.bottom_right_x is None: + print("Pressed 1") + self.top_left_x = xi + self.top_left_y = yi + elif self.top_left_x is not None and self.bottom_right_x is None: + print("Pressed 2") + self.bottom_right_x = xi + self.bottom_right_y = yi + # Find the center + center_x = (self.top_left_x + self.bottom_right_x) / 2 + center_y = (self.top_left_y + self.bottom_right_y) / 2 + # Get the depth at the center + depth = get_depth_at_location(depth_location, center_x, center_y) + crop_width = abs(self.top_left_x - self.bottom_right_x) + crop_height = abs(self.top_left_y - self.bottom_right_y) + target = open("cropboxes.log", 'a') + target.write( + str(center_x) + "," + str(center_y) + "," + str(depth[0]) + "," + str(depth[1]) + "," + str( + depth[2]) + "," + str(crop_width) + "," + str(crop_height) + "\n") + target.close() + self.top_left_x = None + self.top_left_y = None + self.bottom_right_x = None + self.bottom_right_y = None + + handler = EventHandler() + show() + + +def extract_tiltyawdeg(path_to_metadata_xml): + pano = {} + pano_xml = open(path_to_metadata_xml, 'rb') + tree = ET.parse(pano_xml) + root = tree.getroot() + for child in root: + if child.tag == 'projection_properties': + pano[child.tag] = child.attrib + print(pano['projection_properties']['tilt_yaw_deg']) + + return pano['projection_properties']['tilt_yaw_deg'] + + +def get_depth_at_location(path_to_depth_txt, xi, yi): + depth_location = path_to_depth_txt + + filename = depth_location + + print(filename) + + with open(filename, 'rb') as f: + depth = loadtxt(f) + + depth_x = depth[:, 0::3] + depth_y = depth[:, 1::3] + depth_z = depth[:, 2::3] + + val_x, val_y, val_z = interpolated_3d_point(xi, yi, depth_x, depth_y, depth_z) + print 'depth_x, depth_y, depth_z', val_x, val_y, val_z + return val_x, val_y, val_z + + +def predict_crop_size_by_position(x, y, im_width, im_height): + print("Predicting crop size by panorama position") + dist_to_center = math.sqrt((x - im_width / 2) ** 2 + (y - im_height / 2) ** 2) + # Calculate distance from point to center of left edge + dist_to_left_edge = math.sqrt((x - 0) ** 2 + (y - im_height / 2) ** 2) + # Calculate distance from point to center of right edge + dist_to_right_edge = math.sqrt((x - im_width) ** 2 + (y - im_height / 2) ** 2) + + min_dist = min([dist_to_center, dist_to_left_edge, dist_to_right_edge]) + + crop_size = (4.0 / 15.0) * min_dist + 200 + + logging.info("Depth data unavailable; using crop size " + str(crop_size)) + + return crop_size + + +def predict_crop_size(x, y, im_width, im_height, path_to_depth_file): + """ + # Calculate distance from point to image center + dist_to_center = math.sqrt((x-im_width/2)**2 + (y-im_height/2)**2) + # Calculate distance from point to center of left edge + dist_to_left_edge = math.sqrt((x-0)**2 + (y-im_height/2)**2) + # Calculate distance from point to center of right edge + dist_to_right_edge = math.sqrt((x - im_width) ** 2 + (y - im_height/2) ** 2) + + min_dist = min([dist_to_center, dist_to_left_edge, dist_to_right_edge]) + + crop_size = (4.0/15.0)*min_dist + 200 + + print("Min dist was "+str(min_dist)) + """ + crop_size = 0 + try: + depth = get_depth_at_location(path_to_depth_file, x, y) + depth_x = depth[0] + depth_y = depth[1] + depth_z = depth[2] + + distance = math.sqrt(depth_x ** 2 + depth_y ** 2 + depth_z ** 2) + print("Distance is " + str(distance)) + if distance == "nan": + print("Distance is not a number.") + # If no depth data is available, use position in panorama as fallback + # Calculate distance from point to image center + dist_to_center = math.sqrt((x - im_width / 2) ** 2 + (y - im_height / 2) ** 2) + # Calculate distance from point to center of left edge + dist_to_left_edge = math.sqrt((x - 0) ** 2 + (y - im_height / 2) ** 2) + # Calculate distance from point to center of right edge + dist_to_right_edge = math.sqrt((x - im_width) ** 2 + (y - im_height / 2) ** 2) + + min_dist = min([dist_to_center, dist_to_left_edge, dist_to_right_edge]) + + crop_size = (4.0 / 15.0) * min_dist + 200 + + logging.info("Depth data unavailable; using crop size " + str(crop_size)) + else: + # crop_size = (30700.0/37.0)-(300.0/37.0)*distance + # crop_size = 2600 - 220*distance + # crop_size = (5875.0/3.0)-(275.0/3.0)*distance + crop_size = 2050 - 110 * distance + crop_size = 8725.6 * (distance ** -1.192) + if crop_size < 50: + crop_size = 50 + elif crop_size > 1500: + crop_size = 1500 + + logging.info("Distance " + str(distance) + "Crop size " + str(crop_size)) + except IOError: + # If no depth data is available, use position in panorama as fallback + # Calculate distance from point to image center + dist_to_center = math.sqrt((x - im_width / 2) ** 2 + (y - im_height / 2) ** 2) + # Calculate distance from point to center of left edge + dist_to_left_edge = math.sqrt((x - 0) ** 2 + (y - im_height / 2) ** 2) + # Calculate distance from point to center of right edge + dist_to_right_edge = math.sqrt((x - im_width) ** 2 + (y - im_height / 2) ** 2) + + min_dist = min([dist_to_center, dist_to_left_edge, dist_to_right_edge]) + + crop_size = (4.0 / 15.0) * min_dist + 200 + + logging.info("Depth data unavailable; using crop size " + str(crop_size)) + + return crop_size + + +def make_single_crop(path_to_image, sv_image_x, sv_image_y, PanoYawDeg, output_filename, path_to_depth, draw_mark=False): + im_width = GSVImage.GSVImage.im_width + im_height = GSVImage.GSVImage.im_height + im = Image.open(path_to_image) + draw = ImageDraw.Draw(im) + # sv_image_x = sv_image_x - 100 + x = ((float(PanoYawDeg) / 360) * im_width + sv_image_x) % im_width + y = im_height / 2 - sv_image_y + + r = 10 + if draw_mark: + draw.ellipse((x - r, y - r, x + r, y + r), fill=128) + + print("Plotting at " + str(x) + "," + str(y) + " using yaw " + str(PanoYawDeg)) + + # Crop rectangle around label + cropped_square = None + try: + predicted_crop_size = predict_crop_size(x, y, im_width, im_height, path_to_depth) + crop_width = predicted_crop_size + crop_height = predicted_crop_size + print(x, y) + top_left_x = x - crop_width / 2 + top_left_y = y - crop_height / 2 + cropped_square = im.crop((top_left_x, top_left_y, top_left_x + crop_width, top_left_y + crop_height)) + except (ValueError, IndexError) as e: + + predicted_crop_size = predict_crop_size_by_position(x, y, im_width, im_height) + crop_width = predicted_crop_size + crop_height = predicted_crop_size + print(x, y) + top_left_x = x - crop_width / 2 + top_left_y = y - crop_height / 2 + cropped_square = im.crop((top_left_x, top_left_y, top_left_x + crop_width, top_left_y + crop_height)) + cropped_square.save(output_filename) + + return + + +def bulk_extract_crops(path_to_db_export, path_to_gsv_scrapes, destination_dir, mark_label=False): + csv_file = open(path_to_db_export) + csv_f = csv.reader(csv_file) + counter = 0 + no_metadata_fail = 0 + + no_pano_fail = 0 + for row in csv_f: + + if counter == 0: + counter += 1 + continue + + pano_id = row[0] + print(pano_id) + + sv_image_x = float(row[1]) + sv_image_y = float(row[2]) + label_type = int(row[3]) + photographer_heading = float(row[4]) + heading = float(row[5]) + label_id = int(row[7]) + + # Extract Yaw from metadata xml file + pano_xml_path = os.path.join(path_to_gsv_scrapes, pano_id[:2], pano_id + ".xml") + pano_img_path = os.path.join(path_to_gsv_scrapes, pano_id[:2], pano_id + ".jpg") + pano_depth_path = os.path.join(path_to_gsv_scrapes, pano_id[:2], pano_id + ".depth.txt") + print(pano_xml_path) # pano_yaw_deg = float(extract_panoyawdeg(pano_xml_path)) + + # Check that metadata exists for this image; if not skip it + try: + if (os.path.exists(pano_xml_path)): + pano_yaw_deg = float(extract_panoyawdeg(pano_xml_path)) + else: + print("Skipping label due to missing XML data") + logging.warn("Skipped label id " + str(label_id) + " due to missing XML.") + no_metadata_fail += 1 + continue + except (KeyError, ET.ParseError) as e: + print("Skipping label due to invalid XML data") + logging.warn("Skipped label id " + str(label_id) + " due to invalid XML.") + no_metadata_fail += 1 + continue + + print("Photographer heading is " + str(photographer_heading)) + print("Viewer heading is " + str(heading)) + + pano_yaw_deg = 180 - photographer_heading + + print("Yaw:" + str(pano_yaw_deg)) + + # Extract the crop + if os.path.exists(pano_img_path): + counter += 1 + destination_folder = os.path.join(destination_dir, str(label_type)) + if not os.path.isdir(destination_folder): + os.makedirs(destination_folder) + + crop_destination = os.path.join(destination_dir, str(label_type), str(counter) + ".jpg") + if not os.path.exists(crop_destination): + make_single_crop(pano_img_path, sv_image_x, sv_image_y, pano_yaw_deg, crop_destination, pano_depth_path, draw_mark=mark_label) + print("Successfully extracted crop to " + str(counter) + ".jpg") + logging.info(str(counter) + ".jpg" + " " + pano_id + " " + str(sv_image_x) + + " " + str(sv_image_y) + " " + str(pano_yaw_deg) + " " + str(label_id)) + logging.info("---------------------------------------------------") + else: + no_pano_fail += 1 + print("Panorama image not found.") + logging.warn("Skipped label id " + str(label_id) + " due to missing image.") + + print("Finished.") + print(str(no_pano_fail) + " extractions failed because panorama image was not found.") + print(str(no_metadata_fail) + " extractions failed because metadata was not found.") + + +bulk_extract_crops(csv_export_path, gsv_pano_path, destination_path, mark_label=mark_center) +# crop_box_helper("/mnt/umiacs/Panoramas/scrapes4", "labeldata.csv") diff --git a/DownloadRunner.py b/DownloadRunner.py new file mode 100644 index 0000000..87c7d14 --- /dev/null +++ b/DownloadRunner.py @@ -0,0 +1,195 @@ +# !/usr/bin/python + +from SidewalkDB import * + +import os +import httplib +import json +import logging +import cStringIO + +import urllib +import urllib2 + +from PIL import Image +from random import shuffle +import fnmatch + +# **************************************** +# Specify the file storage location below +# **************************************** +storage_location = "/home/anthony/Downloads/test/" + + +def check_download_failed_previously(panoId): + if panoId in open('scrape.log').read(): + return True + else: + return False + + +def fetch_pano_ids_from_webserver(): + unique_ids = [] + conn = httplib.HTTPConnection("sidewalk.umiacs.umd.edu") + conn.request("GET", "/adminapi/labels/panoid") + r1 = conn.getresponse() + data = r1.read() + # print(data) + jsondata = json.loads(data) + + for value in jsondata["features"]: + if value["properties"]["gsv_panorama_id"] not in unique_ids: + unique_ids.append(value["properties"]["gsv_panorama_id"]) + return unique_ids + + +def download_panorama_images(storage_path, pano_list=None): + logging.basicConfig(filename='scrape.log', level=logging.DEBUG) + + if pano_list is None: + unique_ids = fetch_pano_ids_from_webserver() + else: + unique_ids = pano_list + + counter = 0 + failed = 0 + im_dimension = (512 * 26, 512 * 13) + blank_image = Image.new('RGBA', im_dimension, (0, 0, 0, 0)) + base_url = 'http://maps.google.com/cbk?' + shuffle(unique_ids) + for pano_id in unique_ids: + print '-- Extracting images for', pano_id, + + destination_dir = os.path.join(storage_path, pano_id[:2]) + if not os.path.isdir(destination_dir): + os.makedirs(destination_dir) + + filename = pano_id + ".jpg" + + out_image_name = os.path.join(destination_dir, filename) + if os.path.isfile(out_image_name): + print 'File already exists.' + + counter = counter + 1 + print 'Completed ' + str(counter) + ' of ' + str(len(unique_ids)) + continue + + for y in range(13): + for x in range(26): + url_param = 'output=tile&zoom=5&x=' + str(x) + '&y=' + str( + y) + '&cb_client=maps_sv&fover=2&onerr=3&renderer=spherical&v=4&panoid=' + pano_id + url = base_url + url_param + + # Open an image, resize it to 512x512, and paste it into a canvas + req = urllib.urlopen(url) + file = cStringIO.StringIO(req.read()) + im = Image.open(file) + im = im.resize((512, 512)) + + blank_image.paste(im, (512 * x, 512 * y)) + + # Wait a little bit so you don't get blocked by Google + # sleep_in_milliseconds = float(delay) / 1000 + # sleep(sleep_in_milliseconds) + print '.', + print + + # In some cases (e.g., old GSV images), we don't have zoom level 5, + # so Google returns a tranparent image. This means we need to set the + # zoom level to 3. + + # Check if the image is transparent + # http://stackoverflow.com/questions/14041562/python-pil-detect-if-an-image-is-completely-black-or-white + extrema = blank_image.convert("L").getextrema() + if extrema == (0, 0): + print("Panorama %s is an old image and does not have the tiles for zoom level") + temp_im_dimension = (int(512 * 6.5), int(512 * 3.25)) + temp_blank_image = Image.new('RGBA', temp_im_dimension, (0, 0, 0, 0)) + for y in range(3): + for x in range(7): + url_param = 'output=tile&zoom=3&x=' + str(x) + '&y=' + str( + y) + '&cb_client=maps_sv&fover=2&onerr=3&renderer=spherical&v=4&panoid=' + pano_id + url = base_url + url_param + # Open an image, resize it to 512x512, and paste it into a canvas + req = urllib.urlopen(url) + file = cStringIO.StringIO(req.read()) + im = Image.open(file) + im = im.resize((512, 512)) + + temp_blank_image.paste(im, (512 * x, 512 * y)) + + # Wait a little bit so you don't get blocked by Google + # sleep_in_milliseconds = float(delay) / 1000 + # sleep(sleep_in_milliseconds) + print '.', + print + temp_blank_image = temp_blank_image.resize(im_dimension, Image.ANTIALIAS) # resize + temp_blank_image.save(out_image_name, 'jpeg') + else: + blank_image.save(out_image_name, 'jpeg') + print 'Done.' + counter += 1 + print 'Completed ' + str(counter) + ' of ' + str(len(unique_ids)) + return + + +def download_panorama_depthdata(storage_path, decode=True, pano_list=None): + ''' + This method downloads a xml file that contains depth information from GSV. It first + checks if we have a folder for each pano_id, and checks if we already have the corresponding + depth file or not. + ''' + + base_url = "http://maps.google.com/cbk?output=xml&cb_client=maps_sv&hl=en&dm=1&pm=1&ph=1&renderer=cubic,spherical&v=4&panoid=" + + if pano_list is None: + pano_ids = fetch_pano_ids_from_webserver() + else: + pano_ids = pano_list + + for pano_id in pano_ids: + print '-- Extracting depth data for', pano_id, '...', + # Check if the directory exists. Then check if the file already exists and skip if it does. + destination_folder = os.path.join(storage_path, pano_id[:2]) + if not os.path.isdir(destination_folder): + os.makedirs(destination_folder) + + filename = pano_id + ".xml" + destination_file = os.path.join(destination_folder, filename) + if os.path.isfile(destination_file): + print 'File already exists.' + continue + + url = base_url + pano_id + with open(destination_file, 'wb') as f: + req = urllib2.urlopen(url) + for line in req: + f.write(line) + + print 'Done.' + + return + + +def generate_depthmapfiles(path_to_scrapes): + # Iterate through all .xml files in specified path, recursively + for root, dirnames, filenames in os.walk(path_to_scrapes): + for filename in fnmatch.filter(filenames, '*.xml'): + xml_location = os.path.join(root, filename) + print(xml_location) + print(filename) + + # Pano id is XML filename minus the extension + pano_id = filename[:-4] + + # Generate a .depth.txt file for the .xml file + output_file = os.path.join(root, pano_id + ".depth.txt") + call(["./decode_depthmap", xml_location, output_file]) + + +pano_list = fetch_pano_ids_from_webserver() +pano_list = pano_list[:2] + +download_panorama_images(storage_location, pano_list=pano_list) # Trailing slash required +download_panorama_depthdata(storage_location, pano_list=pano_list) +generate_depthmapfiles(storage_location) diff --git a/GSV3DPointImage.py b/GSV3DPointImage.py new file mode 100644 index 0000000..d4dfd68 --- /dev/null +++ b/GSV3DPointImage.py @@ -0,0 +1,653 @@ +from PIL import Image +import math +import numpy as np +import os.path + +import GSVImage +import GSVTopDownImage +from pylab import * +from scipy import ndimage +from utilities import * + +import SidewalkDB as SDB + + +class GSV3DPointImage(object): + def __init__(self, path): + """ + A constructor. This method takes a path to depth.txt. For example: + '../data/GSV/5umV8SPGE1jidFGstzcQDA/' . It then reads the depth.txt and parse + the 3D point data. Note depth.txt contains 512x256x3 points, representing (x, y, z) + in a Cartesian coordinate of 512x256 points (which corresponds to the actual GSV image). + Also note that the origin of the space is (supposedly) the position of a camera on a SV car, + which means the z-value of the origin is 1 to 2 meters higher than the ground. + """ + if type(path) != str and len(str) > 0: + raise ValueError('path should be a string') + if path[-1] != '/': + path += '/' + + ensure_dir(path) + self.path = path + self.pano_id = self.path.split('/')[-1] + self.depth_filename = path + 'depth.txt' + self.image_filename = path + 'images/pano.jpg' + + with open(self.depth_filename, 'rb') as f: + depth = loadtxt(f) + + self.depth = depth + self.px = depth[:, 0::3] + self.py = depth[:, 1::3] + self.pz = depth[:, 2::3] + + # + # Height and width of the 3D point image. + self.height, self.width = self.px.shape + self.gsv_depth_height = self.height + self.gsv_depth_width = self.width + self.gsv_image_width = GSVImage.GSVImage.gsv_image_width + self.gsv_image_height = GSVImage.GSVImage.gsv_image_height + + + # + # Compute normal vectors at each point + # At each point, compute the vectors going outwards from the point by + # taking 8 neighborign points (at the edge of the image you won't get all 8 neighbors) + # Then calculate the normal vector by generating and solving a homogeneous equation (Ax = 0) by svd. + if not os.path.isfile(path + 'normal.txt'): + #if True: + normal_matrix = zeros((self.height, self.width * 3)) + for row_idx in range(self.height): + for col_idx in range(self.width): + vectors = [] + p = array([self.px[row_idx, col_idx], self.py[row_idx, col_idx], self.pz[row_idx, col_idx]]) + if math.isnan(p[0]) or math.isnan(p[1]) or math.isnan(p[2]): + # If p is nan, normal is also nan + normal = array([float('nan'), float('nan'), float('nan')]) + normal_matrix[row_idx, 3 * col_idx] = normal[0] + normal_matrix[row_idx, 3 * col_idx + 1] = normal[1] + normal_matrix[row_idx, 3 * col_idx + 2] = normal[2] + else: + vec_indices = [(row_idx - 1, col_idx - 1), (row_idx - 1, col_idx), (row_idx - 1, col_idx + 1), + (row_idx, col_idx - 1), (row_idx, col_idx + 1), + (row_idx + 1, col_idx - 1), (row_idx + 1, col_idx), (row_idx + 1, col_idx + 1)] + + for idx in vec_indices: + ri = idx[0] + ci = idx[1] + # Check for the corner cases and the case where one of the vector is nan. + if ri > 0 and ri < self.height and ci > 0 and ci < self.width: + if not math.isnan(self.px[idx]) and not math.isnan(self.py[idx]) and not math.isnan(self.pz[idx]): + vec = array([self.px[idx], self.py[idx], self.pz[idx]]) - p + vectors.append(vec) + + # You need at least 3 vectors to calculate normal vectors + if len(vectors) > 2: + vectors = array(vectors) + U, S, Vh = svd(vectors) + normal = array(Vh[-1]) + normal = normal / sqrt(sum(normal ** 2)) + else: + normal = array([float('nan'), float('nan'), float('nan')]) + + + normal_matrix[row_idx, 3 * col_idx] = normal[0] + normal_matrix[row_idx, 3 * col_idx + 1] = normal[1] + normal_matrix[row_idx, 3 * col_idx + 2] = normal[2] + + savetxt(path + 'normal.txt', normal_matrix) + else: + normal_matrix = loadtxt(path + 'normal.txt') + + self.normal = normal_matrix + self.nx = normal_matrix[:, 0::3] + self.ny = normal_matrix[:, 1::3] + self.nz = normal_matrix[:, 2::3] + + return + + def depth_to_png(self, destination, mode='gray'): + """ + Convert the depth map to grayscale png image + """ + destination_dir = '/'.join(destination.split('/')[:-1]) + dir = os.path.dirname(destination_dir) + dir = dir.replace('~', os.path.expanduser('~')) + if not os.path.exists(dir): + raise ValueError('The directory ' + str(dir) + ' does not exist.') + destination = destination.replace('~', os.path.expanduser('~')) + if mode != 'gray': + pass + else: + # + # Grayscale image + depth = sqrt(self.px ** 2 + self.py ** 2 + self.pz ** 2) + + # + # Morphology operation. Get rid of black peppers. + from scipy import ndimage + kernel = array([[-1, -1, -1], [-1, 8, -1], [-1, -1, -1]]) + dr = ndimage.convolve(depth, kernel, mode='reflect', cval=0.0) + # Image.fromarray(array(depth, dtype=float)).show() + # Image.fromarray(array(dr, dtype=float)).show() + # Image.open(self.image_filename).resize((512,256)).convert('L').show() + # + imarray = asarray(Image.open(self.image_filename).resize((512,256)).convert('L')) + + import matplotlib.pyplot as plt + from skimage.data import lena + + # depth = depth / 100 + depth_image = Image.fromarray(array(depth, dtype=float)) + # normal_image = normal_image.convert('1') + depth_image.show() + #normal_image.save(destination, 'PNG') + + def get_overlay_value(self, x, y, overlay="depth", verbose=False): + """ + Given the image coordinate (x, y), return the interpolated overlay value + """ + + if verbose: + print + print 'x, y: ', x, y + + gsv_im_width = self.gsv_image_width + gsv_im_height = self.gsv_image_height + x = (float(x) / gsv_im_width) * self.width + y = (float(y) / gsv_im_height) * self.height + x_floor = math.floor(x) + x_ceil = math.ceil(x) + y_floor = math.floor(y) + y_ceil = math.ceil(y) + + # + # Modify incase points do not form a rectangle + if x_floor == x_ceil: + if x_floor < self.width: + x_ceil += 1 + else: + x_floor -= 1 + if y_floor == y_ceil: + if y_floor < self.height: + y_ceil += 1 + else: + y_floor -= 1 + + if overlay == 'depth': + depth = sqrt(self.px ** 2 + self.py ** 2 + self.pz ** 2) + points = [(x_floor, y_floor, depth[y_floor, x_floor]), + (x_ceil, y_floor, depth[y_floor, x_ceil]), + (x_ceil, y_ceil, depth[y_ceil, x_ceil]), + (x_floor, y_ceil, depth[y_ceil, x_floor]) + ] + value = bilinear_interpolation(x, y, points) + elif overlay == 'normal_z_component': + # Get interpolated z-component of normal vectors + normal_z = self.nz * 255. + points = [(x_floor, y_floor, normal_z[y_floor, x_floor]), + (x_ceil, y_floor, normal_z[y_floor, x_ceil]), + (x_ceil, y_ceil, normal_z[y_ceil, x_ceil]), + (x_floor, y_ceil, normal_z[y_ceil, x_floor]) + ] + value = bilinear_interpolation(x, y, points) + else: + return False + + return value + + def normal_to_png(self, destination, morphology=False): + """ + Convert the normal vector map to monochrome png image. + + """ + + # + # Make sure there is a proper destination + destination_dir = '/'.join(destination.split('/')[:-1]) + dir = os.path.dirname(destination_dir) + dir = dir.replace('~', os.path.expanduser('~')) + if not os.path.exists(dir): + raise ValueError('The directory ' + str(dir) + ' does not exist.') + destination = destination.replace('~', os.path.expanduser('~')) + + # + # Convert a normal bitmap to a png image + normal_array = array(self.nz * 255, dtype=uint8) + if morphology: + binary_array = normal_array > 128 + binary_array = ndimage.binary_opening(binary_array, iterations=2) + binary_array = ndimage.binary_closing(binary_array, iterations=2) + binary_array = binary_array.astype(np.uint8) + normal_array = binary_array * 255 + normal_image = Image.fromarray(normal_array) + normal_image = normal_image.convert('1') + normal_image.save(destination, 'PNG') + return + + def point_to_3d_distance(self, x, y, verbose=True): + """ + This method converts an image point (x, y) into a 3D distance. + """ + # + # First get 4 points available on the depth image that are closest to the passed point + depth_x = self.gsv_depth_width * (float(x) / self.gsv_image_width) + depth_y = self.gsv_depth_height * (float(y) / self.gsv_image_height) + depth_x_floor = floor(depth_x) + depth_y_floor = floor(depth_y) + depth_x_ceil = ceil(depth_x) + depth_y_ceil = ceil(depth_y) + + if depth_x_ceil == self.gsv_depth_width: + depth_x_ceil_ = 0 + else: + depth_x_ceil_ = depth_x_ceil + + points = [ + (depth_x_floor, depth_y_floor), + (depth_x_floor, depth_y_ceil), + (depth_x_ceil, depth_y_floor), + (depth_x_ceil, depth_y_ceil) + ] + + depth_3d_x_values = [(point[0], point[1], self.px[point[1], depth_x_ceil_]) for point in points] + depth_3d_y_values = [(point[0], point[1], self.py[point[1], depth_x_ceil_]) for point in points] + depth_3d_z_values = [(point[0], point[1], self.pz[point[1], depth_x_ceil_]) for point in points] + + depth_3d_x_value = bilinear_interpolation(depth_x, depth_y, depth_3d_x_values) + depth_3d_y_value = bilinear_interpolation(depth_x, depth_y, depth_3d_y_values) + depth_3d_z_value = bilinear_interpolation(depth_x, depth_y, depth_3d_z_values) + distance = math.sqrt(depth_3d_x_value ** 2 + depth_3d_y_value ** 2 + depth_3d_z_value ** 2) + + return distance + + def point_to_latlng(self, x, y, verbose=True): + """ + This method converts an image point (x, y) into a latlng coordinate + + Todos. + - You need to take care of corner cases + """ + # + # First get 4 points available on the depth image that are closest to the passed point + depth_x = self.gsv_depth_width * (float(x) / self.gsv_image_width) + depth_y = self.gsv_depth_height * (float(y) / self.gsv_image_height) + depth_x_floor = floor(depth_x) + depth_y_floor = floor(depth_y) + depth_x_ceil = ceil(depth_x) + depth_y_ceil = ceil(depth_y) + + if depth_x_ceil == self.gsv_depth_width: + depth_x_ceil_ = 0 + else: + depth_x_ceil_ = depth_x_ceil + + points = [ + (depth_x_floor, depth_y_floor), + (depth_x_floor, depth_y_ceil), + (depth_x_ceil, depth_y_floor), + (depth_x_ceil, depth_y_ceil) + ] + + + + + depth_3d_x_values = [(point[0], point[1], self.px[point[1], depth_x_ceil_]) for point in points] + depth_3d_y_values = [(point[0], point[1], self.py[point[1], depth_x_ceil_]) for point in points] + + # + # Corner cases + # if depth_x_ceil == depth_x_floor: + # if depth_y_ceil == depth_y_floor: + # if x is at the left most corner + # print depth_3d_x_values + depth_3d_x_value = bilinear_interpolation(depth_x, depth_y, depth_3d_x_values) + depth_3d_y_value = bilinear_interpolation(depth_x, depth_y, depth_3d_y_values) + distance = math.sqrt(depth_3d_x_value ** 2 + depth_3d_y_value ** 2) + + with open(self.path + 'meta.xml', 'rb') as xml: + tree = ET.parse(xml) + root = tree.getroot() + yaw_deg = float(root.find('projection_properties').get('pano_yaw_deg')) + yaw_deg = (yaw_deg + 180) % 360 + heading = (360 * (float(x) / self.gsv_image_width) + yaw_deg) % 360 + latlng = distance_to_latlng(self.path, distance, heading) + + #latlng = point_to_latlng(self.path, (depth_3d_x_value, depth_3d_y_value)) # This function is is utilities + if verbose: latlng + return latlng + + def save_overlay(self, outfile, overlay_type="depth", mask=None): + """ + Save the image with overlay layer + """ + overlay_layer_intensity = 186 + panorama_image_intensity = 70 + if overlay_type == 'depth': + # Show depth map on top of a GSV image + pass + elif overlay_type == 'vertical_coordinate': + # Highlight points with vertical coordinate less than a threshold + pass + elif overlay_type == 'normal_z_component': + # Show z-component of normal vectors at each point on top of a GSV image + pass + elif overlay_type == 'mask': + # Caution! + # I am assuming mask is 256x512 numpy array with dtype=uint8 + mask_image = Image.fromarray(mask) + mask_image = mask_image.convert('RGB') + overlay = asarray(mask_image) + overlay = overlay.astype('float') + overlay = (overlay / 255) * overlay_layer_intensity + else: + pass + + # Format the panorama image + panorama_image = Image.open(self.image_filename) + panorama_image = panorama_image.resize((512, 256)) + panorama_image = panorama_image.convert('RGB') + panorama_image = asarray(panorama_image) + panorama_image = panorama_image.astype('float') + panorama_image = (panorama_image / 256) * panorama_image_intensity + + # Blend panorama and overlay + im_array = overlay + panorama_image + im_array = im_array.astype('uint8') + im = Image.fromarray(im_array) + im.save(outfile, 'PNG') + return + + def show_overlay(self, im_shape=(1024, 512), overlay_type="depth", transparency=0.7, morphology=False, mask=None): + """ + Show the GSV image with another layer on top of it. + overlay_type: depth, normal_z_component, + """ + overlay_layer_intensity = int(255. * transparency) + panorama_image_intensity = int(255. * (1 - transparency)) + if overlay_type == 'depth': + # Show depth map on top of a GSV image + depth_threshold = 50 # 50 + + depth = sqrt(self.px ** 2 + self.py ** 2 + self.pz ** 2) + temp_mask = isnan(depth) + depth[isnan(depth)] = depth_threshold + depth[depth > depth_threshold] = depth_threshold + depth = depth - depth.min() + depth = depth / depth.max() + depth = 255 - 255 * depth + + # Treating an array with NaN + # http://stackoverflow.com/questions/5480694/numpy-calculate-averages-with-nans-removed + #minimum_depth = ma.masked_array(depth, isnan(depth)).min() + #depth = depth - minimum_depth + #maximum_depth = ma.masked_array(depth, isnan(depth)).max() + #depth = 256 - 256 * (depth / maximum_depth) + + overlay = depth.astype('uint8') + overlay = Image.fromarray(overlay).convert('RGB') + overlay = asarray(overlay) + overlay = overlay.astype('float') + overlay[:,:,2] = 0. + overlay[:,:,0] = 255. + overlay[:,:,1] = 255. - overlay[:,:,1] + overlay[temp_mask, :] = 0 + overlay = (overlay / 256) * overlay_layer_intensity + + elif overlay_type == 'vertical_coordinate': + # Highlight points with vertical coordinate less than a threshold + vertical_axis_threshold = -2.5 + + overlay = self.pz + overlay[isnan(overlay)] = 0 + overlay[overlay > vertical_axis_threshold] = 0 + overlay[overlay < vertical_axis_threshold] = 255 + overlay = overlay.astype('uint8') + + overlay = Image.fromarray(overlay).convert('RGB') + overlay = asarray(overlay) + overlay = overlay.astype('float') + overlay = (overlay / 256) * overlay_layer_intensity + elif overlay_type == 'normal_z_component': + # Show z-component of normal vectors at each point on top of a GSV image + + normal_array = array(self.nz * 255, dtype=uint8) + if morphology: + binary_array = normal_array > 128 + binary_array = ndimage.binary_opening(binary_array, iterations=2) + binary_array = ndimage.binary_closing(binary_array, iterations=2) + binary_array = binary_array.astype(np.uint8) + normal_array = binary_array * 255 + + normal_image = Image.fromarray(normal_array) + + normal_image = normal_image.convert('RGB') + + # Blending. Superimpose one image on top of another. + # http://stackoverflow.com/questions/5605174/python-pil-function-to-divide-blend-two-images + overlay = asarray(normal_image) + #overlay.flags.writeable = True + #thresh = median(overlay[overlay>0]) - 0.5 + #overlay[overlay > thresh] = 255 + #overlay[overlay < thresh] = 0 + overlay = overlay.astype('float') + overlay = (overlay / 255) * overlay_layer_intensity + elif overlay_type == 'mask': + # Caution! + # I am assuming mask is 512x1024 numpy array with dtype=uint8 + mask_image = Image.fromarray(mask.astype(np.uint8)) + mask_image = mask_image.convert('RGB') + overlay = asarray(mask_image) + overlay = overlay.astype(np.float) + overlay = (overlay / 255) * overlay_layer_intensity + else: + panorama_image = Image.open(self.image_filename) + panorama_image = panorama_image.resize((512, 256)) + panorama_image.show() + return + + # Format the panorama image + # Resize Image.fromarray(overlay.astype(np.uint8), 'RGB').resize((13312, 6656)).show() + panorama_image = Image.open(self.image_filename) + panorama_image = panorama_image.resize(im_shape) + panorama_image = panorama_image.convert('RGB') + panorama_image = asarray(panorama_image) + panorama_image = panorama_image.astype(np.float) + panorama_image = (panorama_image / 256) * panorama_image_intensity + + # Blend panorama and overlay + if overlay.shape != panorama_image.shape: + overlay = Image.fromarray(overlay.astype(np.uint8), 'RGB').resize((panorama_image.shape[1], panorama_image.shape[0])) + overlay = np.asarray(overlay) + im_array = overlay + panorama_image + im_array = im_array.astype(np.uint8) + im = Image.fromarray(im_array) + im.show() + return + + def top_down_to_pano(self, src_image, size=(1024, 512)): + """ + This method maps the top down view back to the original panorama view. + The method can also map the mask image (currently only gray scale images) to orignal + panorama view. You can pass the image parameter to do so. + + Fist, the method creates the mapping from top down image to panorama image + Based on the mapping that you've created, you map the pixel values. + """ + # + # Get the output image size and scale relative to the 3d point data + pano_im_width = size[0] + pano_im_height = size[1] + scale = pano_im_width / 512 + + # + # Get th size of 3d point image + height_3d, width_3d = self.px.shape + x_grid, y_grid = mgrid[0:height_3d, 0:width_3d] + + # + # Interpolate the depth data to match the output image size + grid_points = np.array([x_grid.flatten(), y_grid.flatten()]).T + grid_points = grid_points * scale + x_image_grid, y_image_grid = mgrid[0:pano_im_height, 0:pano_im_width] + + px = self.px.flatten() + py = self.py.flatten() + pz = self.pz.flatten() + + from scipy.interpolate import griddata + px_new = griddata(grid_points, px, (x_image_grid, y_image_grid), method='linear') + py_new = griddata(grid_points, py, (x_image_grid, y_image_grid), method='linear') + pz_new = griddata(grid_points, pz, (x_image_grid, y_image_grid), method='linear') + + # + # Vectorize + x = arange(0, pano_im_width) + y = arange(0, pano_im_height) + xx, yy = np.meshgrid(x, y) + + def pixel_mapping(x, y): + """ numpy function. This takes numpy arrays x and y. This function then maps pixels in + a Street View image to a 2-D top down image. + """ + depth_3d_x_values = px_new[y, x] + depth_3d_y_values = py_new[y, x] + + # + # Create a masked array + # http://blog.remotesensing.io/2013/05/Dealing-with-no-data-using-NumPy-masked-array-operations/ + # invert + # http://docs.scipy.org/doc/numpy/reference/generated/numpy.invert.html + masked_depth_3d_x_values = np.ma.masked_invalid(depth_3d_x_values) + masked_depth_3d_y_values = np.ma.masked_invalid(depth_3d_y_values) + + # + # Intersect x-mask and y-mask and create new masked_depth values just in + # case we have something like (x, y, z) = (nan, 1, nan), which we want to neglect + mask = masked_depth_3d_x_values.mask + masked_depth_3d_y_values.mask + masked_depth_3d_x_values = np.ma.array(depth_3d_x_values, mask=mask) + masked_depth_3d_y_values = np.ma.array(depth_3d_y_values, mask=mask) + + # + # Get the top down image pixel coordinate + top_down_y_pixel_coordinates = (masked_depth_3d_y_values * 10 + GSVTopDownImage.GSVTopDownImage.height / 2).astype(int) + top_down_x_pixel_coordinates = (masked_depth_3d_x_values * 10 + GSVTopDownImage.GSVTopDownImage.width / 2).astype(int) + + # + # Eliminate pixels that are outside of the top down image + mask = (top_down_x_pixel_coordinates >= 1000) + (top_down_x_pixel_coordinates < 0) + (top_down_y_pixel_coordinates >= 1000) + (top_down_y_pixel_coordinates < 0) + top_down_x_pixel_coordinates = np.ma.array(top_down_x_pixel_coordinates, mask=mask) + top_down_y_pixel_coordinates = np.ma.array(top_down_y_pixel_coordinates, mask=mask) + + # + # Make arrays of top_down_image pixel coordinates + mask = np.invert(top_down_x_pixel_coordinates.mask) + xx = top_down_x_pixel_coordinates.data[mask] + yy = top_down_y_pixel_coordinates.data[mask] + + return y_image_grid[mask], x_image_grid[mask], src_image[yy, xx] + + y_coordinates, x_coordinates, pixel_values = pixel_mapping(xx, yy) + + blank_image = np.zeros((size[1], size[0], pixel_values[0].size)).astype(np.uint8) + blank_image[x_coordinates, y_coordinates, :] = array([array(pixel_values)]).T + + return blank_image.astype(np.uint8) + +def batch_decode_normal_image(): + # + # Retrive task panoramas and store them into TaskPanoramaTable + sql = "SELECT * FROM TaskPanoramas WHERE TaskDescription=%s" + with SDB.SidewalkDB() as db: + records = db.query(sql, ('PilotTask_v2_MountPleasant')) + + # + # The constructor of GSV3DPointImage creates a normal image + for record in records: + pano_id = record[1] + gsv = GSV3DPointImage('../data/GSV/' + pano_id + '/') + +def batch_get_normal_mask(): + # This function retrieves panorama ids and + # + # Retrive task panoramas and store them into TaskPanoramaTable + sql = "SELECT * FROM TaskPanoramas WHERE TaskDescription=%s" + with SDB.SidewalkDB() as db: + records = db.query(sql, ('PilotTask_v2_MountPleasant')) + + # + # The constructor of GSV3DPointImage creates a normal image + for record in records: + pano_id = record[1] + gsv_3d_point_image = GSV3DPointImage('../data/GSV/' + pano_id + '/') + filename = '../data/temp/masked_images_2/' + pano_id + '.png' + gsv_3d_point_image.normal_to_png(filename, morphology=True) + return + + +def depth_features_for_UIST_scenes(): + filename = "../data/temp/StreetPixels.csv" + out_filename = "../data/temp/SceneDepth.csv" + import csv + with open(out_filename, "w") as of: + writer = csv.writer(of) + writer.writerow(["PanoramaId", "DepthMean", "DepthMedian", "DepthStdev", "DepthMax", "DepthMin"]) + with open(filename, 'r') as f: + reader = csv.reader(f) + next(f) + for line in reader: + #print line + path = '../data/GSV/%s/' % (line[0]) + gsv_3d = GSV3DPointImage(path) + arr = np.ma.masked_invalid(np.sqrt(gsv_3d.px ** 2 + gsv_3d.py ** 2 + gsv_3d.pz ** 2)) + depth_mean = arr.mean() + depth_median = np.ma.median(arr) + depth_stdev = arr.std() + depth_max = arr.max() + depth_min = arr.min() + # sprint "%s,%0.2f,%0.2f,%0.2f,%0.2f,%0.2f" % (line[0], depth_mean, depth_median, depth_stdev, depth_max, depth_min) + writer.writerow([line[0], depth_mean, depth_median, depth_stdev, depth_max, depth_min]) + + +def script(): + # panorama_ids = [ + # "_AUz5cV_ofocoDbesxY3Kw", + # "0C6PG3Zpuwz11kZKfG_vUg", + # "D-2VNbhqOqYAKTU0hFneIw", + # "-dlUzxwCI_-k5RbGw6IlEg", + # "inU6vka3Fzz8xFkD50OkcA", + # "MoooLIfIUXg4Id8UfCobmw", + # "pc68VExSXDKWHaXqmYyNew", + # "q3hvBgZVagEv2wgrPgJ6og", + # "QaklWtS6F4qXTdmXzynhxQ", + # "qbPS050BhVdsW9Jh7bbLRA", + # "QF4m3RsaH8qNayiRq7GeSQ", + # "SSkaybYviuU_u0MHRoZMJw", + # "tId8wkF-MITThzEUOlIWXA", + # "U0koc8H_RE_E2eM-DFsoYQ", + # "uCxyTYfPDd7efWvxnYQSSg", + # "VsP0gcbV2Yv-WX_NlbVHTQ", + # "vVlss1eLpiYU9AsLi-Didg", + # "ZtTnE0fh5firrU4yDAhCCw", + # "ZWl97D059PURIRRbvb5AEA" + # ] + panorama_ids = [ + # "h7ZW0_VasRt3vhevz1mjeg", + "Aw67wmndIEG7DT3jLFXH6g" + ] + for panorama_id in panorama_ids: + im_3d = GSV3DPointImage('../data/GSV/%s/' % panorama_id) + im_3d.show_overlay((4096, 2048)) + + return + +if __name__ == "__main__": + print "GSV3DPointImage.py" + script() + + + + + + diff --git a/GSVImage.py b/GSVImage.py new file mode 100644 index 0000000..5f33f6e --- /dev/null +++ b/GSVImage.py @@ -0,0 +1,1106 @@ +try: + import cv2 + import cv2.cv as cv +except ImportError, e: + cv2 = None + cv = None +# http://stackoverflow.com/questions/9226258/why-does-python-cv2-modules-depend-on-old-cv + +import os +import unittest + +import GSV3DPointImage as gsv3d + +from copy import deepcopy +from PIL import Image, ImageDraw +from pylab import * +from random import randint +from utilities import * + +try: + from xml.etree import cElementTree as ET +except ImportError, e: + from xml.etree import ElementTree as ET + +gsv_image_width = 13312 +gsv_image_height = 6656 +im_width = gsv_image_width +im_height = gsv_image_height + +class GSVImage(object): + gsv_image_width = 13312 + gsv_image_height = 6656 + im_width = gsv_image_width + im_height = gsv_image_height + + + def __init__(self, path): + """ + A constructor. This method takes a path to GSV data files. + For example: ../data/GSV/5umV8SPGE1jidFGstzcQDA/ + """ + self.im_width = 13312 # deprecated + self.im_height = 6656 # deprecated + self.gsv_image_width = 13312 + self.gsv_image_height = 6656 + self.pano_id = path.split('/')[-1] + + ensure_dir(path) + self.path = path + + def adjust_aspect_ratio(self, boundingbox, aspect_ratio): + patch_height = boundingbox['y_max'] - boundingbox['y_min'] + patch_width = patch_height * aspect_ratio + + if boundingbox['boundary']: + x_min = boundingbox['x_min'] - GSVImage.im_width + else: + x_min = boundingbox['x_min'] + boundingbox_width = boundingbox['x_max'] - x_min + + # If patch_width is larger than the boundingbox_width, use the patch_width. + # Otherwise use the bounding_box width and modify the height of the bounding box + if boundingbox_width < patch_width: + x_center = (boundingbox['x_max'] + x_min) / 2 + + # Set x_max and x_min + x_min = x_center - patch_width / 2 + + if x_min < 0: + boundingbox['boundary'] = True + boundingbox['x_min'] = GSVImage.im_width + x_min + else: + boundingbox['boundary'] = False + boundingbox['x_min'] = x_min + + boundingbox['x_max'] = x_center + patch_width / 2 + + boundingbox['x_max'] = int(boundingbox['x_max']) + boundingbox['x_min'] = int(boundingbox['x_min']) + else: + patch_width = boundingbox_width + patch_height = boundingbox_width / aspect_ratio + y_center = (boundingbox['y_max'] + boundingbox['y_min']) / 2 + boundingbox['y_max'] = int(y_center + patch_height / 2) + boundingbox['y_min'] = int(y_center - patch_height / 2) + return boundingbox + + def crop_user_bounding_boxes(self, boundingboxes, output_filename, aspect_ratio=None, show_image=False, + overwrite=False, verbose=False): + filename = self.path + 'images/pano.jpg' + # ensure_dir(directory) + im = Image.open(filename) + filenames = [] + for idx, boundingbox in enumerate(boundingboxes): + # + # Crop + filename = output_filename + '_' + str(idx) + '.jpg' + filenames.append(filename) + if verbose: + print filename + if os.path.isfile(filename) and not overwrite: + print filename, 'File alread exists' + return + + if aspect_ratio is not None: + boundingbox = self.adjust_aspect_ratio(boundingbox, aspect_ratio) + + if boundingbox['boundary']: + patch_height = boundingbox['y_max'] - boundingbox['y_min'] + patch_width = (GSVImage.im_width - boundingbox['x_min']) + boundingbox['x_max'] + im_dimension = (patch_width, patch_height) + patch = Image.new('RGBA', im_dimension, (0, 0, 0, 0)) + + # Crop and paste the first half of the bounding box + box = (int(boundingbox['x_min']), int(boundingbox['y_min']), GSVImage.im_width, + int(boundingbox['y_max'])) + partial = im.crop(box) + patch.paste(partial, (0, 0)) + + # Crop and paste the last half of the bounding box + box = (0, int(boundingbox['y_min']), int(boundingbox['x_max']), int(boundingbox['y_max'])) + partial = im.crop(box) + patch.paste(partial, (GSVImage.im_width - boundingbox['x_min'], 0)) + else: + # http://stackoverflow.com/questions/1076638/trouble-using-python-pil-library-to-crop-and-save-image + box = (int(boundingbox['x_min']), int(boundingbox['y_min']), int(boundingbox['x_max']), + int(boundingbox['y_max'])) + patch = im.crop(box) + + if show_image: + patch.show() + + patch.save(filename, 'JPEG') + + return filenames + + def crop_negative_bounding_boxes(self, boundingboxes, negative_filenames, overlay='normal_z_component', + overlap_ratio=0.5, aspect_ratio=None, show_image=True, verbose=False): + """ + This method takes a list of ground truth bounding boxes (*) and a filename header (**). + (*) E.g., [{ + 'boundary': False, + 'x_min': 300, + 'x_max': 700, + 'y_min': 300, + 'y_max': 600 + }, ...]), and + (**) "../data/temp/negative/ + """ + if len(boundingboxes) == 0: + return + num_ground = 0 + num_patches = 1 + + # + # Find the smallest boudning box and use it as the cropping size + filename = self.path + 'images/pano.jpg' + # ensure_dir(directory) + im = Image.open(filename) + gsv_3d_point_image = gsv3d.GSV3DPointImage(self.path) + positive_boundingboxes = deepcopy(boundingboxes) + + area = float('inf') + for boundingbox in boundingboxes: + if boundingbox['boundary']: + bb_width = boundingbox['x_max'] + (self.im_width - boundingbox['x_min']) + else: + bb_width = boundingbox['x_max'] - boundingbox['x_min'] + + bb_area = bb_width * (boundingbox['y_max'] - boundingbox['y_min']) + if area > bb_area: + area = bb_area + bb = boundingbox + + if bb['boundary']: + crop_width = bb['x_max'] - (bb['x_min'] - self.im_width) + else: + crop_width = bb['x_max'] - bb['x_min'] + crop_height = bb['y_max'] - bb['y_min'] + + # + # For each label, crop 9 negative image patches + half_crop_height = int(crop_height / 2) + half_crop_width = int(crop_width / 2) + + padding_top = int((crop_height + 1) / 2) + padding_bottom = padding_top + padding_left = int((crop_width + 1) / 2) + padding_right = padding_left + negative_bounding_boxes = [] + + for current_bb in boundingboxes: + # + # Randomly choose a point on the image. + # Make sure the negative example does not overlap with other bounding boxes. + # If overlay is specified, crop image patches only from the masked area + bb_w = current_bb['x_max'] - current_bb['x_min'] # Bounding box width + bb_h = current_bb['y_max'] - current_bb['y_min'] # Bounding box height + num_around_ramp = num_patches - num_ground + for idx in range(num_ground): + print 'neg', len(negative_bounding_boxes) # debug + cropped = False + while not cropped: + y = randint(padding_top, self.im_height - padding_bottom) + x = randint(padding_left, self.im_width - padding_right) + crop_bb = { + 'boundary': False, + 'x_min': x - half_crop_width, + 'x_max': x + half_crop_width, + 'y_min': y - half_crop_height, + 'y_max': y + half_crop_height + } + # + # Debug + #print 'neg', len(negative_bounding_boxes), negative_bounding_boxes + #box = (crop_bb['x_min'], crop_bb['y_min'], crop_bb['x_max'], crop_bb['y_max']) + #patch = im.crop(box) + #patch.show() + + do_continue = False + for boundingbox in positive_boundingboxes: + areaoverlap = bounding_box_area_overlap(crop_bb, boundingbox) + if areaoverlap > overlap_ratio: + # Area overlap too big + do_continue = True + break + if do_continue: + continue + for boundingbox in negative_bounding_boxes: + areaoverlap = bounding_box_area_overlap(crop_bb, boundingbox) + if areaoverlap > overlap_ratio: + do_continue = True + break + if do_continue: + continue + # if overlay is specified, check if (x, y) is on the masked area + if overlay: + overlay_value = gsv_3d_point_image.get_overlay_value(x, y, overlay='normal_z_component') + overlay_threshold = 10 + if isnan(overlay_value) or overlay_value < overlay_threshold: + # print 'Weak overlay value: ', overlay_value + do_continue = True + if do_continue: + continue + negative_bounding_boxes.append(crop_bb) + cropped = True + # + # Crop negative image patches around the current boundning box (current_bb) + for idx in range(num_ground, num_patches): + cropped = False + num_loop = 0 + + + while not cropped: + # + # Crop an image patch. + # Make sure the area overlap between a negative patch and a positive patch is below a threshold + # Try to collect patches around positive patches, but if it is not possible, collect patches from + # the masked area (ground) + num_loop += 1 + #y = randint(current_bb['y_min'] - 2 * padding_top, current_bb['y_max'] + 2 * padding_bottom) + y = randint(current_bb['y_min'] - crop_height, current_bb['y_max']) # + crop_height) + + if idx % 2 == 0: + # x = randint(current_bb['x_min'] - padding_left, current_bb['x_min']) + x = randint(current_bb['x_min'] - 1.5 * crop_width, current_bb['x_min'] - 0.5 * crop_width) + if x < padding_left or x > self.im_width - padding_right: + x = randint(current_bb['x_max'], current_bb['x_max'] + 0.5 * crop_width) + else: + # x = randint(current_bb['x_max'], current_bb['x_max'] + padding_left) + x = randint(current_bb['x_max'], current_bb['x_max'] + 0.5 * crop_width) + if x < padding_left or x > self.im_width - padding_right: + x = randint(current_bb['x_min'] - 1.5 * crop_width, current_bb['x_min'] - 0.5 * crop_width) + + if num_loop > 200: + print "random" + y = randint(padding_top, self.im_height - padding_bottom) + x = randint(padding_left, self.im_width - padding_right) + + if x < padding_left or x > self.im_width - padding_right or y < padding_top or y > self.im_height - padding_bottom: + print 'bounding box out of range' + continue + crop_bb = { + 'boundary': False, + 'x_min': x - half_crop_width, + 'x_max': x + half_crop_width, + 'y_min': y - half_crop_height, + 'y_max': y + half_crop_height + } + # + # Debug + #print 'neg', len(negative_bounding_boxes), negative_bounding_boxes + #box = (crop_bb['x_min'], crop_bb['y_min'], crop_bb['x_max'], crop_bb['y_max']) + #patch = im.crop(box) + #patch.show() + + do_continue = False + for pos_bb in positive_boundingboxes: + areaoverlap = bounding_box_area_overlap(crop_bb, pos_bb) + # box = (pos_bb['x_min'], pos_bb['y_min'], pos_bb['x_max'], pos_bb['y_max']) + # patch = im.crop(box) + # patch.show() + + if areaoverlap > overlap_ratio: + # Area overlap too big + do_continue = True + break + + if do_continue: + continue + + for neg_bb in negative_bounding_boxes: + areaoverlap = bounding_box_area_overlap(crop_bb, neg_bb) + if areaoverlap > 0.5: + do_continue = True + break + + if do_continue: + continue + # if overlay is specified, check if (x, y) is on the masked area + if overlay: + overlay_value = gsv_3d_point_image.get_overlay_value(x, y, overlay='normal_z_component') + overlay_threshold = 10 + if isnan(overlay_value) or overlay_value < overlay_threshold: + # print 'Weak overlay value: ', overlay_value + do_continue = True + + if do_continue: + continue + + negative_bounding_boxes.append(crop_bb) + cropped = True + + # + # Crop bounding boxes for negative image patches + filenames = [] + for i, boundingbox in enumerate(negative_bounding_boxes): + outline = str(boundingbox['x_min']) + ' ' + str(boundingbox['y_min']) + ' ' + outline += str(boundingbox['x_max']) + ' ' + str(boundingbox['y_min']) + ' ' + outline += str(boundingbox['x_max']) + ' ' + str(boundingbox['y_max']) + ' ' + outline += str(boundingbox['x_min']) + ' ' + str(boundingbox['y_max']) + + box = (boundingbox['x_min'], boundingbox['y_min'], boundingbox['x_max'], boundingbox['y_max']) + patch = im.crop(box) + #if show_image: + #patch.show() + + + filename = negative_filenames + '_' + str(i) + '.jpg' + filenames.append(filename) + + patch.save(filename, 'JPEG') + + return filenames, negative_bounding_boxes + + def crop_negative_image_patches(self, outlines, output_filename, overlap_ratio=0.5, aspect_ratio=None, overlay=None, show_image=False): + """ + This method takes outlines, and crops image patches that do not overlap with outlines. + If overlay (mask) is specified, crop negative patches only from masked region. + """ + filename = self.path + '/images/pano.png' + # ensure_dir(directory) + im = Image.open(filename) + gsv_3d_point_image = gsv3d.GSV3DPointImage(self.path) + + boundingboxes = [] + for outline in outlines: + print '---' + xys = outline.strip().split(' ') + xs = [] + ys = [] + points = [] + for x, y in zip(xys[0::2], xys[1::2]): + p = user_point_to_sv_image_point(self.path, {'x': x, 'y': y}) + p = (int(x), int(y)) + points.append(p) + + # Compute the bounding box of the passed label points + # If aspect_ratio is passed, format the bounding box so the + # aspect ratio of the cropped image will be consistent across + # all the image patches + boundingbox = sv_image_points_to_bounding_box(points) + if aspect_ratio and (type(aspect_ratio) == int or type(aspect_ratio) == float): + patch_height = boundingbox['y_max'] - boundingbox['y_min'] + patch_width = patch_height * aspect_ratio + + if boundingbox['boundary']: + x_min = boundingbox['x_min'] - GSVImage.im_width + else: + x_min = boundingbox['x_min'] + boundingbox_width = boundingbox['x_max'] - x_min + + # If patch_width is larger than the boundingbox_width, use the patch_width. + # Otherwise use the bounding_box width and modify the height of the bounding box + if boundingbox_width < patch_width: + x_center = (boundingbox['x_max'] + x_min) / 2 + + # Set x_max and x_min + x_min = x_center - patch_width / 2 + + if x_min < 0: + boundingbox['boundary'] = True + boundingbox['x_min'] = GSVImage.im_width + x_min + else: + boundingbox['boundary'] = False + boundingbox['x_min'] = x_min + + boundingbox['x_max'] = x_center + patch_width / 2 + + boundingbox['x_max'] = int(boundingbox['x_max']) + boundingbox['x_min'] = int(boundingbox['x_min']) + else: + patch_width = boundingbox_width + patch_height = boundingbox_width / aspect_ratio + y_center = (boundingbox['y_max'] + boundingbox['y_min']) / 2 + boundingbox['y_max'] = int(y_center + patch_height / 2) + boundingbox['y_min'] = int(y_center - patch_height / 2) + boundingboxes.append(boundingbox) + + if len(boundingboxes) == 0: + return + # + # Find the smallest boudning box and use it as the cropping size + area = float('inf') + for boundingbox in boundingboxes: + + if boundingbox['boundary']: + bb_width = boundingbox['x_max'] + (self.im_width - boundingbox['x_min']) + else: + bb_width = boundingbox['x_max'] - boundingbox['x_min'] + + print 'bb list: ', boundingbox + print 'bb width: ', bb_width + + bb_area = bb_width * (boundingbox['y_max'] - boundingbox['y_min']) + if area > bb_area: + area = bb_area + bb = boundingbox + + print 'smallest bb: ', bb + print 'smallest bb width: ', bb_width + if bb['boundary']: + crop_width = bb['x_max'] - (bb['x_min'] - self.im_width) + else: + crop_width = bb['x_max'] - bb['x_min'] + crop_height = bb['y_max'] - bb['y_min'] + + # + # For each label, crop 9 negative image patches + half_crop_height = int(crop_height / 2) + half_crop_width = int(crop_width / 2) + + padding_top = int((crop_height + 1) / 2) + padding_bottom = padding_top + padding_left = int((crop_width + 1) / 2) + padding_right = padding_left + print 'crop width, crop height', crop_width, crop_height + negative_bounding_boxes = [] + for outline_idx, outline in enumerate(outlines): + # + # Randomly choose a point on the image. + # Make sure the negative example does not overlap with other bounding boxes. + # If overlay is specified, crop image patches only from the masked area + for idx in range(9): + cropped = False + while not cropped: + y = randint(padding_top, self.im_height - padding_bottom) + x = randint(padding_left, self.im_width - padding_right) + print "image size: width, height", self.im_width, self.im_height + print "crop size: width, height = ", crop_width, crop_height + print "paddings: ", padding_left, padding_top, padding_right, padding_bottom + #print '-- Random: x, y = ', x, y + crop_bb = { + 'boundary': False, + 'x_min': x - half_crop_width, + 'x_max': x + half_crop_width, + 'y_min': y - half_crop_height, + 'y_max': y + half_crop_height + } + print crop_bb + + crop_bb_does_not_overlap_with_other_bb = True + for boundingbox in boundingboxes: + areaoverlap = bounding_box_area_overlap(crop_bb, boundingbox) + if areaoverlap > overlap_ratio: + # Area overlap too big + continue + for boundingbox in negative_bounding_boxes: + areaoverlap = bounding_box_area_overlap(crop_bb, boundingbox) + if areaoverlap > overlap_ratio: + continue + + # if overlay is specified, check if (x, y) is on the masked area + if overlay: + overlay_value = gsv_3d_point_image.get_overlay_value(x, y, overlay='normal_z_component') + overlay_threshold = 10 + if isnan(overlay_value) or overlay_value < overlay_threshold: + # print 'Weak overlay value: ', overlay_value + continue + + negative_bounding_boxes.append(crop_bb) + cropped = True + + # Crop bounding boxes + for i, boundingbox in enumerate(negative_bounding_boxes): + outline = str(boundingbox['x_min']) + ' ' + str(boundingbox['y_min']) + ' ' + outline += str(boundingbox['x_max']) + ' ' + str(boundingbox['y_min']) + ' ' + outline += str(boundingbox['x_max']) + ' ' + str(boundingbox['y_max']) + ' ' + outline += str(boundingbox['x_min']) + ' ' + str(boundingbox['y_max']) + + box = (boundingbox['x_min'], boundingbox['y_min'], boundingbox['x_max'], boundingbox['y_max']) + patch = im.crop(box) + if show_image: + patch.show() + + filename = output_filename + '_' + str(i) + '.jpg' + patch.save(filename, 'JPEG') + + return + + def crop_user_outline(self, outline, output_filename, aspect_ratio=None, show_image=False): + """ + This method takes an outline, output_filename, + """ + + if type(outline) != str: + raise ValueError('First parameter should be str.') + if type(output_filename) != str: + raise ValueError('Second parameter should be str.') + + outline_length = len(outline.strip().split(' ')) + if outline_length % 2 != 0 or outline_length < 6: + raise ValueError('Illegal number of outline point coordinates') + + return crop_user_outline(self.path, outline, output_filename, aspect_ratio, show_image) + + def crop_user_outlines(self, outlines, output_filename, aspect_ratio=None, show_image=False): + """ + This method takes outlines, output_filename, + """ + path = self.path + + if type(outlines) != list: + raise ValueError('First parameter should be list.') + if type(output_filename) != str: + raise ValueError('Second parameter should be str.') + + for i, outline in enumerate(outlines): + outline_length = len(outline.strip().split(' ')) + if outline_length % 2 != 0 or outline_length < 6: + raise ValueError('Illegal number of outline point coordinates') + + filename = output_filename + '_' + str(i) + '.jpg' + crop_user_outline(self.path, outline, filename, aspect_ratio, show_image) + return + + def get_image_latlng(self): + """ + This method returns a latlng position of this image. + """ + xml = open(self.path + 'meta.xml', 'rb') + tree = ET.parse(xml) + data = tree.find('data_properties').attrib + lat = float(data['lat']) + lng = float(data['lng']) + return lat, lng + + def get_pano_id(self): + """ + Return the panorama id of this image + """ + return self.pano_id + + def plot_bounding_boxes(self, bounding_boxes, image_size=None, width=5, outline='red', output_file=None): + """ + This method renders detected bounding boxes on an image. Each bounding box should have + the format of [(x1, y1), (x2, y2)]. If image size is give, then shrink the image accordingly + + Color Examples: "red", "#92d050", "#00b050 + Open image file + http://scikit-image.org/docs/dev/auto_examples/applications/plot_morphology.html + """ + filename = self.path + 'images/pano.jpg' + + """ + import matplotlib.pyplot as plt + from skimage import io + im = io.imread(filename) + plt.imshow(im) + """ + + im = Image.open(filename) + + # + # Draw bounding boxes + # http://effbot.org/imagingbook/imagedraw.htm + draw = ImageDraw.Draw(im) + + for box in bounding_boxes: + x1 = box[0][0] + y1 = box[0][1] + x2 = box[1][0] + y2 = box[1][1] + for i in range(0, width): + if ((x1 + i) >= (x2 - i)) or ((y1 + i) >= (y2 - i)): + break + new_box = [(x1 + i, y1 + i), (x2 - i, y2 - i)] + draw.rectangle(new_box, outline=outline) + + del draw + if image_size: + im = im.resize(image_size) + im.show() + + if output_file: + im.save(output_file, 'PNG') + return + + def plot_user_outline(self, outline, user_point=True): + """ + This function plots a user provided GSV outline (a set of points) on an actual GSV image. + + :param path: + A path to a directory where GSV data are stored + :type path: + str. E.g., '12082 -490 12118 -411 12017 -388 11764 -374 11764 -420 11852 -462 11934 -490' + Note that + :param outline: + A set of GSV image points provided through the Street View Labeler Interface. Set user_point to True if + you are passing image points. + + :type outline: + list. + """ + # Go through every 2 items + # http://stackoverflow.com/questions/5389507/iterating-over-every-two-elements-in-a-list + filename = self.path + 'images/pano.jpg' + im = Image.open(filename) + draw = ImageDraw.Draw(im) + + xys = outline.strip().split(' ') + xs = [] + ys = [] + im_xys = [] + + points = [] + for x, y in zip(xys[0::2], xys[1::2]): + if user_point: + p = self.user_point_to_sv_image_point({'x': x, 'y': y}) + else: + p = (int(x), int(y)) + points.append(p) + r = 10 + draw.ellipse((p[0]-r, p[1]-r, p[0]+r, p[1]+r), fill='rgb(200,0,0)') + r = 5 + draw.ellipse((p[0]-r, p[1]-r, p[0]+r, p[1]+r), fill='white') + + for i, p in enumerate(points): + draw.line((points[i-1][0], points[i-1][1], points[i][0], points[i][1]), fill='red') + + figure() + imshow(im) + title('PanoId: ' + self.path[:-1].split('/')[-1]) + + show() + return + + def plot_user_outlines(self, outlines): + """ + This function plots a set of user provided GSV outlines + + :param outline: + A set of user provided GSV image points + :type outline: + list. + """ + # Go through every 2 items + # http://stackoverflow.com/questions/5389507/iterating-over-every-two-elements-in-a-list + filename = self.path + '/images/pano.png' + im = Image.open(filename) + draw = ImageDraw.Draw(im) + + for outline in outlines: + xys = outline.strip().split(' ') + xs = [] + ys = [] + im_xys = [] + + points = [] + for x, y in zip(xys[0::2], xys[1::2]): + p = self.user_point_to_sv_image_point({'x': x, 'y': y}) + points.append(p) + r = 20 + draw.ellipse((p[0]-r, p[1]-r, p[0]+r, p[1]+r), fill='rgb(200,0,0)') + r = 15 + draw.ellipse((p[0]-r, p[1]-r, p[0]+r, p[1]+r), fill='white') + + for i, p in enumerate(points): + draw.line((points[i-1][0], points[i-1][1], points[i][0], points[i][1]), fill='red', width=10) + + figure() + imshow(im) + title('PanoId: ' + self.path[:-1].split('/')[-1]) + + show() + + def user_point_to_sv_image_point(self, point): + """ + This function converts a GSV image point coordinate provided by user through CSI interface to + a true GSV image coordinate + """ + return user_point_to_sv_image_point(self.path, point) + + def show(self, size=False): + """ + This method shows the corresponding GSV panorama image + + options can take size + http://www.learnpython.org/Multiple_Function_Arguments + """ + if os.path.isfile(self.path + 'images/pano.jpg'): + im = Image.open(self.path + 'images/pano.jpg') + + if size and type(size) == tuple: + im = im.resize(size, Image.ANTIALIAS) + im.show() + else: + raise Exception(self.path + 'images/pano.jpg does not exist') + return + + def sv_image_point_to_user_point(self, point, image_size=None): + """ + This method converts sv_image point (x, y) on a street view image (e.g., points that constitutes + a curb ramp bounding box detected by a program) into user point (or point on SV image on SV API). + """ + return sv_image_point_to_user_point(self.path, point, image_size=image_size) + + def sv_image_point_to_pov(self, point, image_size=None): + """ + This method converts finds pov for the point + """ + return sv_image_point_to_pov(self.path, point, image_size=image_size) + + + def sv_image_points_to_bounding_box(self, points): + return sv_image_points_to_bounding_box(points) + + +def sv_image_point_to_pov(path, point, image_size=None): + """ + Pass an image point coordinate and convert it to pov + """ + x = point[0] + y = point[1] + + if image_size: + w = image_size[0] + h = image_size[1] + x = x * GSVImage.gsv_image_width / w + y = y * GSVImage.gsv_image_height / h + + xml = open(path + 'meta.xml', 'rb') + tree = ET.parse(xml) + root = tree.getroot() + pano_yaw_deg = float(root.find('projection_properties').get('pano_yaw_deg')) + + heading = (360. * (float(x) / GSVImage.gsv_image_width)) + (pano_yaw_deg - 180) + heading = (heading + 360) % 360 + + pitch = 90 - 180 * (float(y) / GSVImage.gsv_image_height) + + pov = {'heading': heading, + 'pitch': pitch, + 'zoom': 1 + } + return pov + +def crop_user_outline(input_path, outline, output_filename, aspect_ratio=None, show_image=False): + """ + This function crops an image patch from the GSV image passed in input_path by forming a bounding box from + the passed outline. It will save the output image patch in output_path + + :param input_path: A path to a directory where GSV files are stored. + :type input_path: str. + :param output_filename: A name of an output image path + :type output_filename: str. + :param outline: A set of user provided GSV image points + :type outline: list. + :param aspect_ratio: An aspect ratio (width/height) of image patches to crop. (1:aspect_ratio) + :param show_image: A flag to indicate whether the function should show the cropped image or not. + :type show_image: bool. + """ + filename = input_path + '/images/pano.png' + # ensure_dir(directory) + im = Image.open(filename) + + xys = outline.strip().split(' ') + xs = [] + ys = [] + points = [] + for x, y in zip(xys[0::2], xys[1::2]): + p = user_point_to_sv_image_point(input_path, {'x': x, 'y': y}) + points.append(p) + + # Compute the bounding box of the passed label points + # If aspect_ratio is passed, format the bounding box so the + # aspect ratio of the cropped image will be consistent across + # all the image patches + boundingbox = sv_image_points_to_bounding_box(points) + if aspect_ratio and (type(aspect_ratio) == int or type(aspect_ratio) == float): + patch_height = boundingbox['y_max'] - boundingbox['y_min'] + patch_width = patch_height * aspect_ratio + + if boundingbox['boundary']: + x_min = boundingbox['x_min'] - GSVImage.im_width + else: + x_min = boundingbox['x_min'] + boundingbox_width = boundingbox['x_max'] - x_min + + # If patch_width is larger than the boundingbox_width, use the patch_width. + # Otherwise use the bounding_box width and modify the height of the bounding box + if boundingbox_width < patch_width: + x_center = (boundingbox['x_max'] + x_min) / 2 + + # Set x_max and x_min + x_min = x_center - patch_width / 2 + + if x_min < 0: + boundingbox['boundary'] = True + boundingbox['x_min'] = GSVImage.im_width + x_min + else: + boundingbox['boundary'] = False + boundingbox['x_min'] = x_min + + boundingbox['x_max'] = x_center + patch_width / 2 + + boundingbox['x_max'] = int(boundingbox['x_max']) + boundingbox['x_min'] = int(boundingbox['x_min']) + else: + patch_width = boundingbox_width + patch_height = boundingbox_width / aspect_ratio + y_center = (boundingbox['y_max'] + boundingbox['y_min']) / 2 + boundingbox['y_max'] = int(y_center + patch_height / 2) + boundingbox['y_min'] = int(y_center - patch_height / 2) + + # boundingbox['boundary'] indicates whether the bounding box goes over the boundary of a SV image + # If it does (i.e., boundingbox['boundary'] == True), take care of it. Otherwise just crop it. + if boundingbox['boundary']: + patch_height = boundingbox['y_max'] - boundingbox['y_min'] + patch_width = (GSVImage.im_width - boundingbox['x_min']) + boundingbox['x_max'] + im_dimension=(patch_width, patch_height) + patch = Image.new('RGBA', im_dimension, (0, 0, 0, 0)) + + # Crop and paste the first half of the bounding box + box = (boundingbox['x_min'], boundingbox['y_min'], GSVImage.im_width, boundingbox['y_max']) + partial = im.crop(box) + patch.paste(partial, (0, 0)) + + # Crop and paste the last half of the bounding box + box = (0, boundingbox['y_min'], boundingbox['x_max'], boundingbox['y_max']) + partial = im.crop(box) + patch.paste(partial, (GSVImage.im_width - boundingbox['x_min'], 0)) + else: + # http://stackoverflow.com/questions/1076638/trouble-using-python-pil-library-to-crop-and-save-image + box = (boundingbox['x_min'], boundingbox['y_min'], boundingbox['x_max'], boundingbox['y_max']) + patch = im.crop(box) + + if show_image: + figure() + imshow(patch) + title('PanoId: ' + input_path[:-1].split('/')[-1]) + show() + + patch.save(output_filename, 'JPEG') + + return + + +def user_point_to_sv_image_point(path, point): + """ + This function converts a GSV image point coordinate provided by user through CSI interface to + a true GSV image coordinate + + :param path: + A path to a directory where GSV files for the target panorama are stored + E.g., path: '../data/GSV/MO1a01Frnzs4IgoGxo1XvQ/' + :type path: str. + :param point: + A GSV image point provided by a user through CSI interface. + If string, pass it like '3152 50', i.e., 'x-coordinate y-coordinate' + If dict, pass it like '{x: 3152, y:50}' + :type point: str. + """ + + if type(point) == str: + p = point.strip().split(' ') + if len(p) != 2: + raise ValueError("plot_user_point() expect the second input format to be 'x y'") + x = int(p[0].strip()) + y = int(p[1].strip()) + elif type(point) == dict: + x = int(point['x']) + y = int(point['y']) + elif type(point) == tuple: + x = int(point[0]) + y = int(point[1]) + else: + x = int(point[0]) + y = int(point[1]) + + # Extract the sv meta data. + xml = open(path + 'meta.xml', 'rb') + tree = ET.parse(xml) + root = tree.getroot() + pano_yaw_deg = float(root.find('projection_properties').get('pano_yaw_deg')) + # tilt_yaw_deg = float(root.find('projection_properties').get('tilt_yaw_deg')) + yaw_deg = pano_yaw_deg # - tilt_yaw_deg + + im_width = GSVImage.im_width + im_height = GSVImage.im_height + + # Translate a point to adjust its coordinate to the local image. + y = y + x = ((540 - yaw_deg) / 360) * im_width + x + x = x % im_width + y = im_height / 2 - y + + x = int(x) + y = int(y) + return x, y + + +def sv_image_points_to_bounding_box(points): + ''' + :param points: + A list of GSV image points provided by users through the CSI interface. + E.g., [(2049, 3818), (2085, 3739), (1984, 3716), (1731, 3702), (1731, 3748), (1819, 3790), (1901, 3818)] + :returns: + A bounding box. x_min/x_max/y_min/y_max and whether the boudning box is split by the vertical image boundary or not. + If the bounding box is split by the vertical image boundary, then x_max is the largest value between x = [0, gsv_im_width/2], and + x_min is the smallest value between x = [gsv_im_width/2, gsv_im_width] + ''' + boundary = False + x_min = 1000000 + x_max = -1 + y_min = 1000000 + y_max = -1000000 + + # + # Check if the outline points are split by the vertical image boundary. + for point in points: + if point[0] < x_min: + x_min = point[0] + if point[0] > x_max: + x_max = point[0] + if point[1] < y_min: + y_min = point[1] + if point[1] > y_max: + y_max = point[1] + + if x_max - x_min > 3500: + boundary = True + + # + # Split in two cases. + if boundary: + x_min = 1000000 + x_max = -1 + for point in points: + # x min and max + if point[0] < GSVImage.im_width / 2: + if point[0] > x_max: + x_max = point[0] + if point[0] > GSVImage.im_width / 2: + if point[0] < x_min: + x_min = point[0] + + return {'boundary' : boundary, + 'x_min': x_min, + 'x_max': x_max, + 'y_min': y_min, + 'y_max': y_max} + + +def sv_image_point_to_user_point(file_path, point, image_size=None): + """ + This method converts a Street View image point (x, y) on a street view image (e.g., one of bounding points that + forms a curb ramp bounding box that is detected by a detector) into user point (or point on SV API image). + """ + (x, y) = point + x = int(x) + y = int(y) + + im_width = GSVImage.im_width + im_height = GSVImage.im_height + + if image_size: + w = image_size[0] + h = image_size[1] + x = x * im_width / w + y = y * im_height / h + + # + # Extract the sv meta data. + xml = open(file_path + 'meta.xml', 'rb') + tree = ET.parse(xml) + root = tree.getroot() + pano_yaw_deg = float(root.find('projection_properties').get('pano_yaw_deg')) + yaw_deg = pano_yaw_deg + # + # Translate a point to adjust its coordinate to the local image. + y = im_height / 2 - y + x = x - ((540 - yaw_deg) / 360) * im_width + x = (x + im_width) % im_width + + + x = int(x) + y = int(y) + return x, y + + +def bounding_box_area_overlap(bb1, bb2): + """ + This function takes two bounding boxes and calculates the area overlap. + + Todo. There is a corner case that I'm not capturing. + There could be a boundary label and non-boundary label that could intersect. + """ + if bb1['boundary']: + bb1_x_min = bb1['x_min'] - GSVImage.im_width + else: + bb1_x_min = bb1['x_min'] + if bb2['boundary']: + bb2_x_min = bb2['x_min'] - GSVImage.im_width + else: + bb2_x_min = bb2['x_min'] + + # + # To make the function simple, manipulate data so the bb2 will always have smaller x_min + if bb1_x_min < bb2_x_min: + temp = bb1_x_min + bb1_x_min = bb2_x_min + bb2_x_min = temp + temp = bb1 + bb1 = bb2 + bb2 = temp + + bb1_x_max = bb1['x_max'] + bb2_x_max = bb2['x_max'] + + bb1_y_max = bb1['y_max'] + bb2_y_max = bb2['y_max'] + bb1_y_min = bb1['y_min'] + bb2_y_min = bb2['y_min'] + + if bb2_x_max < bb1_x_min: + overlap = 0 + elif bb2_x_max < bb1_x_max: + if bb2_y_min < bb1_y_min: + if bb2_y_max < bb1_y_min: + overlap = 0 + elif bb2_y_max < bb1_y_max: + overlap = (bb2_y_max - bb1_y_min) * (bb2_x_max - bb1_x_min) + else: + overlap = (bb1_y_max - bb1_y_min) * (bb2_x_max - bb1_x_min) + elif bb2_y_min < bb1_y_max: + # + # bb1_y_min < bb2_y_min < bb1_y_max: + if bb2_y_max < bb1_y_max: + overlap = (bb2_y_max - bb2_y_min) * (bb2_x_max - bb1_x_min) + else: + overlap = (bb1_y_max - bb2_y_min) * (bb2_x_max - bb1_x_min) + else: + # + # bb2_y_min > bb1_y_max + overlap = 0 + else: + # + # bb2_x_max > bb1_x_max + if bb2_y_min < bb1_y_min: + if bb2_y_max < bb1_y_min: + overlap = 0 + elif bb2_y_max < bb1_y_max: + overlap = (bb2_y_max - bb1_y_min) * (bb1_x_max - bb1_x_min) + else: + overlap = (bb1_y_max - bb1_y_min) * (bb1_x_max - bb1_x_min) + elif bb2_y_min < bb1_y_max: + # + # bb1_y_min < bb2_y_min < bb1_y_max: + if bb2_y_max < bb1_y_max: + overlap = (bb2_y_max - bb2_y_min) * (bb1_x_max - bb1_x_min) + else: + overlap = (bb1_y_max - bb2_y_min) * (bb1_x_max - bb1_x_min) + else: + # + # bb2_y_min > bb1_y_max + overlap = 0 + + area_bb1 = (bb1_y_max - bb1_y_min) * (bb1_x_max - bb1_x_min) + area_bb2 = (bb2_y_max - bb2_y_min) * (bb2_x_max - bb2_x_min) + area = area_bb1 + area_bb2 - overlap + return float(overlap) / area + +if __name__ == '__main__': + print "GSVImage.py" + diff --git a/GSVScraper.py b/GSVScraper.py new file mode 100644 index 0000000..597d8b5 --- /dev/null +++ b/GSVScraper.py @@ -0,0 +1,882 @@ +''' +Created on May 10, 2013 + +@author: kotarohara +''' +import cStringIO +import math +import os +import subprocess +import urllib +import urllib2 +import GSVImage + +from copy import deepcopy +# from GSVImage import user_point_to_sv_image_point, sv_image_points_to_bounding_box +from PIL import Image, ImageDraw +from pylab import * +from subprocess import call, check_output +from time import sleep +from utilities import * + +from SidewalkDB import * +from pymysql import OperationalError + +try: + from xml.etree import cElementTree as ET +except ImportError, e: + from xml.etree import ElementTree as ET + +class GSVScraper(object): + ''' + classdocs + ''' + def __init__(self, data_dir='../data/GSV/', database="sidewalk"): + ''' + Constructor + ''' + self.pano_ids = [] + self.coordinates = [] + self.data_dir = data_dir + try: + self.db = SidewalkDB(database=database) + except OperationalError, e: + self.db = None + return + + def decode_depthmap(self): + """ + Decode depth.xml + """ + for pano_id in self.pano_ids: + if os.path.isfile(self.data_dir + pano_id + '/depth.txt'): + print 'File already exists.' + else: + decode_depthmap('../data/GSV/' + pano_id + '/depth.xml', '../data/GSV/' + pano_id + '/depth.txt', verbose=True) + return + + def depth_first_search(self, depth=5, bounding_box=None, poly=None, verbose=False): + """ + This function + + :param depth: The depth of the map. + """ + if len(self.pano_ids) <= 0: + raise ValueError('No pano id provided') + all_panoramas = [] + + seed_panoramas = deepcopy(self.pano_ids) + + for pano in seed_panoramas: + if verbose: + print 'Seed pano: ', pano + print 'Extracting connected panoramas', + visited_panoramas = [] + passed_panorama = [] + panorama_stack = [pano] + + while len(panorama_stack) > 0: + if verbose: print '.', + if verbose: + print panorama_stack + curr_pano = panorama_stack[-1] + if curr_pano not in visited_panoramas: + visited_panoramas.append(curr_pano) + + # If the length of the stack is higher than the depth, + # pop the stack and continue + # Otherwise, see the top panorama in the stack + # Mark the top panorama as visited + # Find the next panorama that is not visited and push it on the stack + if len(panorama_stack) > depth: + leaf_pano = panorama_stack.pop() + else: + links = self.get_pano_links(curr_pano) + all_done = True + for link_pano in links: + if link_pano not in visited_panoramas and link_pano not in passed_panorama: + # visited_panoramas.append(link_pano) + + # + # Check if link_pano is a Google provided panoramas as opposed to user provied panoramas + if not self.pano_is_provided_by_users(link_pano): + # Check if link_pano is in the bounding box + if poly: + if self.pano_is_in_polygon(link_pano, poly): + panorama_stack.append(link_pano) + else: + passed_panorama.append(link_pano) + elif bounding_box: + if self.pano_is_in_bounding_box(link_pano, bounding_box): + panorama_stack.append(link_pano) + else: + passed_panorama.append(link_pano) + else: + panorama_stack.append(link_pano) + else: + continue + all_done = False + break + if all_done: + panorama_stack.pop() + all_panoramas += visited_panoramas + if verbose: + print + print + + all_panoramas = list(set(all_panoramas)) + return all_panoramas + + def get_intersections(self, panoramas, thresh=2): + """ + This method takes a list of panorama ids and returns the ones that has more than 2(thresh) links (intersectiosn) + """ + + intersections = [] + for pano in panoramas: + links = self.get_pano_links(pano) + links = filter(lambda x: self.pano_is_provided_by_google(x), links) + if len(links) > thresh: + intersections.append(pano) + return intersections + + def get_pano_coordinate(self, pano_id): + """ + This method takes a panorama id and returns the lat/lng coordinate + + :param pano_id: panorama id + """ + self.get_pano_metadata([pano_id]) + xml = open(self.data_dir + pano_id + '/meta.xml', 'rb') + tree = ET.parse(xml) + data = tree.find('data_properties').attrib + lat = float(data['lat']) + lng = float(data['lng']) + return (lat, lng) + + def get_projection_properties(self, pano_id): + """ + This method takes a panorama id and returns the projection properties including yaw_degree and + """ + self.get_pano_metadata([pano_id]) + xml = open(self.data_dir + pano_id + '/meta.xml', 'rb') + tree = ET.parse(xml) + data = tree.find('projection_properties').attrib + yaw = float(data['pano_yaw_deg']) + pitch = float(data['tilt_pitch_deg']) + return (yaw, pitch) + + def get_pano_depthdata(self, decode=True, delay=1000.): + ''' + This method downloads a xml file that contains depth information from GSV. It first + checks if we have a folder for each pano_id, and checks if we already have the corresponding + depth file or not. + ''' + + base_url = "http://maps.google.com/cbk?output=xml&cb_client=maps_sv&hl=en&dm=1&pm=1&ph=1&renderer=cubic,spherical&v=4&panoid=" + for pano_id in self.pano_ids: + print '-- Extracting depth data for', pano_id, '...', + # Check if the directory exists. Then check if the file already exists and skip if it does. + ensure_dir(self.data_dir + pano_id) + if os.path.isfile(self.data_dir + pano_id + '/depth.xml'): + print 'File already exists.' + continue + + url = base_url + pano_id + with open(self.data_dir + pano_id + '/depth.xml', 'wb') as f: + req = urllib2.urlopen(url) + for line in req: + f.write(line) + + # Wait a little bit so you don't get blocked by Google + sleep_in_seconds = float(delay) / 1000 + sleep(sleep_in_seconds) + + print 'Done.' + + if decode: + self.decode_depthmap() + + return + + def get_pano_id(self, lat, lng, verbose=False): + """ + This method gets the closest panorama id from the given latlng coordinate + """ + url_header = 'http://cbk0.google.com/cbk?output=xml&ll=' + url = url_header + str(lat) + ',' + str(lng) + pano_id = None + + try: + pano_xml = urllib.urlopen(url) + tree = ET.parse(pano_xml) + root = tree.getroot() + + pano_id = root.find('data_properties').get('pano_id') + except AttributeError: + pass + # Wait a little bit so you don't get blocked by Google + sleep_in_milliseconds = float(1000) / 1000 + sleep(sleep_in_milliseconds) + + return pano_id + + def get_pano_image(self, delay=100.): + ''' + This function collects panorama images and stitch them together + With zoom=5, there are 26x13 images. + http://stackoverflow.com/questions/7391945/how-do-i-read-image-data-from-a-url-in-python + ''' + 'http://maps.google.com/cbk?output=tile&zoom=5&x=1&y=12&cb_client=maps_sv&fover=2&onerr=3&renderer=spherical&v=4&panoid=rP_WcfFFp3V23ESWa59p4Q' + im_dimension = (512 * 26, 512 * 13) + blank_image = Image.new('RGBA', im_dimension, (0, 0, 0, 0)) + base_url = 'http://maps.google.com/cbk?' + + for pano_id in self.pano_ids: + print '-- Extracting images for', pano_id, + ensure_dir(self.data_dir + pano_id) + ensure_dir(self.data_dir + pano_id + '/images/') + out_image_name = self.data_dir + pano_id + '/images/pano.jpg' + if os.path.isfile(out_image_name): + print 'File already exists.' + continue + + for y in range(13): + for x in range(26): + url_param = 'output=tile&zoom=5&x=' + str(x) + '&y=' + str(y) + '&cb_client=maps_sv&fover=2&onerr=3&renderer=spherical&v=4&panoid=' + pano_id + url = base_url + url_param + + # Open an image, resize it to 512x512, and paste it into a canvas + req = urllib.urlopen(url) + file = cStringIO.StringIO(req.read()) + im = Image.open(file) + im = im.resize((512, 512)) + + blank_image.paste(im, (512 * x, 512 * y)) + + # Wait a little bit so you don't get blocked by Google + sleep_in_milliseconds = float(delay) / 1000 + sleep(sleep_in_milliseconds) + print '.', + print + + # In some cases (e.g., old GSV images), we don't have zoom level 5, so + # we need to set the zoom level to 3. + if array(blank_image)[:, :, :3].sum() == 0: + print "Panorama %s is an old image and does not have the tiles for zoom level" + temp_im_dimension = (int(512 * 6.5), int(512 * 3.25)) + temp_blank_image = Image.new('RGBA', temp_im_dimension, (0, 0, 0, 0)) + for y in range(3): + for x in range(7): + url_param = 'output=tile&zoom=3&x=' + str(x) + '&y=' + str(y) + '&cb_client=maps_sv&fover=2&onerr=3&renderer=spherical&v=4&panoid=' + pano_id + url = base_url + url_param + # Open an image, resize it to 512x512, and paste it into a canvas + req = urllib.urlopen(url) + file = cStringIO.StringIO(req.read()) + im = Image.open(file) + im = im.resize((512, 512)) + + temp_blank_image.paste(im, (512 * x, 512 * y)) + + # Wait a little bit so you don't get blocked by Google + sleep_in_milliseconds = float(delay) / 1000 + sleep(sleep_in_milliseconds) + print '.', + print + temp_blank_image = temp_blank_image.resize(im_dimension, Image.ANTIALIAS) # resize + temp_blank_image.save(out_image_name, 'jpeg') + else: + blank_image.save(out_image_name, 'jpeg') + print 'Done.' + return + + def get_pano_links(self, pano): + """ + This method takes a panorama id and returns a set of linked panorama ids + + :param pano: A GSV panorama id + """ + self.get_pano_metadata([pano]) + xml = open(self.data_dir + pano + '/meta.xml', 'rb') + tree = ET.parse(xml) + links = tree.findall('annotation_properties/link') + + linked_panos = [] + for link in links: + linked_panos.append(link.attrib['pano_id']) + #linked_panos.append(link.attrib) + """ + yaw_deg = float(root.find('projection_properties').get('pano_yaw_deg')) + lat = float(root.find('data_properties').get('lat')) + lng = float(root.find('data_properties').get('lng')) + yaw_radian = radians(yaw_deg) + rotation_matrix = array([[cos(yaw_radian), -sin(yaw_radian)], [sin(yaw_radian), cos(yaw_radian)]]) + """ + return linked_panos + + def get_pano_metadata(self, pano_ids=None, delay=1000., save_as_file=True, target_dir=None, verbose=False): + """ + This function collects Google Street View panorama metadata that corresponds to the nearest GSV panoramas. + E.g., + http://jamiethompson.co.uk/web/2010/05/15/google-streetview-static-api/ + http://cbk0.google.com/cbk?output=xml&ll=51.494966,-0.146674 + """ + + if not pano_ids: + pano_ids = self.pano_ids + elif type(pano_ids) != list: + raise ValueError('pano_ids must be a list of GSV panorama ids') + + + api_header = 'http://cbk0.google.com/cbk?output=xml' + for pano_id in pano_ids: + if verbose: + print '-- Extracting metadata for', pano_id, '...', + # Check if the directory exists. Then check if the file already exists and skip if it does. + # Check file: http://stackoverflow.com/questions/82831/how-do-i-check-if-a-file-exists-using-python + if target_dir is None: + target_dir = self.data_dir + ensure_dir(target_dir + pano_id + '/') + # ensure_dir(self.data_dir + pano_id + '/') + if os.path.isfile(target_dir + pano_id + '/meta.xml'): + if verbose: + print 'File already exists.' + continue + + url = api_header + '&panoid=' + pano_id + req = urllib2.urlopen(url) + + if save_as_file: + + with open(target_dir + pano_id + '/meta.xml', 'w+') as my_file: + for line in req: + my_file.write(line) + + + + # Wait a little bit so you don't get blocked by Google + sleep_in_milliseconds = float(delay) / 1000 + sleep(sleep_in_milliseconds) + if verbose: + print 'Done.' + + return + + def pano_is_provided_by_users(self, link_pano): + """ + This method checks if the panorama has level_id attribute (which exists only in user provided panorama images.) + """ + self.get_pano_metadata([link_pano]) + xml = open(self.data_dir + link_pano + '/meta.xml', 'rb') + try: + # print link_pano + tree = ET.parse(xml) + except ET.ParseError: + print link_pano + raise + + if tree.find('levels') != None: + return True + elif tree.find('data_properties/attribution_name') != None: + return True + else: + return False + + def pano_is_provided_by_google(self, link_pano): + """ + This method checks if the panorama is provided by Google + """ + if self.pano_is_provided_by_users(link_pano): + return False + else: + return True + + + def pano_is_in_bounding_box(self, link_pano, bounding_box): + """ + :param link_pano: A panorama id. + :param bounding_box: + A tuple of latitudes and longitudes that defines a bounding box of which part of the map you wnat to look at. + Format is (min_lat, max_lat, min_lng, max_lng) + E.g., (38.896231,38.897934,-77.029755,-77.025109) + """ + min_lat = bounding_box[0] + max_lat = bounding_box[1] + min_lng = bounding_box[2] + max_lng = bounding_box[3] + lat, lng = self.get_pano_coordinate(link_pano) + """self.get_pano_metadata([link_pano]) + xml = open(self.data_dir + link_pano + '/meta.xml', 'rb') + tree = ET.parse(xml) + data = tree.find('data_properties').attrib + lat = float(data['lat']) + lng = float(data['lng'])""" + # coffee + + if min_lat < lat and max_lat > lat and min_lng < lng and max_lng > lng: + return True + else: + return False + + def pano_is_in_polygon(self, pano_id, poly): + """ + This method takes a panorama id, retrieves the latlng coordinate, and checks if + the coordinate is in the polygon. + """ + lat, lng = self.get_pano_coordinate(pano_id) + return point_inside_polygon(lng, lat, poly) + + def set_pano_ids(self, pano_ids): + ''' + This method sets pano_ids of your interest. This method creates data folders + to store all bunch of crap you will download from Street View and name it with + pano_id. + ''' + self.pano_ids = pano_ids + + for pano_id in pano_ids: + ensure_dir(self.data_dir + pano_id + '/') + +""" + Helper functions +""" + + +def read_depth_file(path, show_image=True): + """ + This function reads a 3D point-cloud data from the file generated by ./decode_depthmap. (depth.txt) + The depth.txt contains (x, y, z) points in the Cartesian coordinate. The origin is set to the position + of a camera on a SV car. So the z of the origin is about 1 or 2 meters higher from the ground. + + Todo. Need to investigate if the 3D data takes into account of camera tilt at non-flat locations. + """ + filename = path + 'depth.txt' + image_name = path + 'images/pano.png' + pano_im = array(Image.open(image_name)) + + with open(filename, 'rb') as f: + depth = loadtxt(f) + + depth_x = depth[:, 0::3] + depth_y = depth[:, 1::3] + depth_z = depth[:, 2::3] + + figure() + im = imshow(pano_im) + fig = gcf() + ax = gca() + + class EventHandler: + def __init__(self): + self.prev_x = 0 + self.prev_y = 0 + self.prev_z = 0 + fig.canvas.mpl_connect('button_press_event', self.onpress) + + def onpress(self, event): + ''' + On press, do bilinear interpolation + http://en.wikipedia.org/wiki/Bilinear_interpolation + http://stackoverflow.com/questions/8661537/how-to-perform-bilinear-interpolation-in-python + ''' + if event.inaxes != ax: + return + xi, yi = (int(round(n)) for n in (event.xdata, event.ydata)) + # value = im.get_array()[xi,yi] + # color = im.cmap(im.norm(value)) + + val_x, val_y, val_z = interpolated_3d_point(xi, yi, depth_x, depth_y, depth_z) + print 'depth_x, depth_y, depth_z', val_x, val_y, val_z + + user_points = [(val_x, val_y)] + latlngs = points_to_latlng(path, user_points) + lat = latlngs[0][0] + lng = latlngs[0][1] + print 'lat, lng:', lat, lng + print 'Distance from previous point:', math.sqrt(math.pow((val_x - self.prev_x),2) + math.pow((val_y - self.prev_y), 2) + math.pow((val_z - self.prev_z), 2)) + self.prev_x = val_x + self.prev_y = val_y + self.prev_z = val_z + + handler = EventHandler() + show() + return + + +def decode_depthmap(file_in, file_out, verbose=True): + """ + This function executes ./decode_depthmap . The decode_depthmap retrieves 3D point-cloud data + from the file_in (depth.xml) and spits out the result. + + call function + http://stackoverflow.com/questions/89228/calling-an-external-command-in-python + """ + if verbose: print '-- Decoding depth data...', + if os.path.isfile(file_out): + print 'File already exists.' + return + + import platform + + operating_system = platform.system() + + if operating_system == 'Windows': + # Windows + # + # Caution!!! I have worked on this for a couple of hours, but I could not run the decode_depthmap_win.exe + # from PyLab using subprocess.call. Quick walk around is to run the python script from the cmd.exe + # Will investigate the solution in future. + # http://stackoverflow.com/questions/3022013/windows-cant-find-the-file-on-subprocess-call + # http://stackoverflow.com/questions/10236260/subprocess-pydev-console-vs-cmd-exe + + # pwd = os.path.dirname(os.path.abspath(__file__)) + # bin_dir = "\\".join(pwd.split("\\")[:-1]) + "\\bin" + # my_env = os.environ.copy() + # my_env["PATH"] += os.pathsep + bin_dir + + call(["../bin/decode_depthmap_win.exe", file_in, file_out]) + #popen = subprocess.Popen(["../bin/decode_depthmap_win.exe", file_in, file_out], creationflags=subprocess.CREATE_NEW_CONSOLE) + #popen.wait() + #out = check_output([bin_dir + "\decode_depthmap_win.exe", file_in, file_out], env=my_env) + #if verbose: print out + else: + # Mac + call(["../bin/decode_depthmap", file_in, file_out]) + print 'Done.' + return + + +def plot_user_points(pano_id): + # Image constant + records = [] + sql = """ +SELECT LabelTypeId, svImageX, svImageY, PanoYawDeg FROM Label +INNER JOIN LabelPosition +ON Label.LabelId = LabelPosition.LabelId +INNER JOIN Panorama +ON LabelGSVPanoramaId = Panorama.GSVPanoramaId +INNER JOIN PanoramaProjectionProperty +ON Panorama.GSVPanoramaId = PanoramaProjectionProperty.GSVPanoramaId +WHERE LabelGSVPanoramaId = %s +AND LabelTypeId=1 + """ + from BusStopDB import BusStopDB + with BusStopDB() as db: + records = db.query(sql, (pano_id)) + + im_width = GSVImage.GSVImage.im_width + im_height = GSVImage.GSVImage.im_height + PanoYawDeg = float(records[0][3]) + + filename = '../data/GSV/' + pano_id + '/images/pano.png' + im = Image.open(filename) + draw = ImageDraw.Draw(im) + for i, record in enumerate(records): + # PIL draw circle + # http://stackoverflow.com/questions/2980366/draw-circle-pil-python + # http://www.pythonware.com/library/pil/handbook/imagedraw.htm + # User input data + sv_image_x = int(record[1]) - 100 + sv_image_y = int(record[2]) + x = ((PanoYawDeg / 360) * im_width + sv_image_x) % im_width + y = im_height / 2 - sv_image_y + r = 30 + draw.ellipse((x-r, y-r, x+r, y+r), fill=128) + + figure() + imshow(im) + show() + return + +""" + Helper functions +""" +def batch_decode_depth_data(): + # + # Retrive task panoramas and store them into TaskPanoramaTable + sql = "SELECT * FROM TaskPanoramas WHERE TaskDescription=%s" + with SidewalkDB() as db: + records = db.query(sql, ('PilotTask_v2_MountPleasant')) + + pano_ids = [record[1] for record in records] + scraper = GSVScraper() + scraper.set_pano_ids(pano_ids) + scraper.get_pano_depthdata() + +def format_pano_metadata(pano_id, delay=1000.0, verbose=False): + """ + This function takes a pano_id (e.g., dWeBDzGMXwQv5fu1GoNy8Q) and + returns a Google Street View panorama metadata that corresponds to the nearest GSV panorama. + E.g., + + http://jamiethompson.co.uk/web/2010/05/15/google-streetview-static-api/ + http://cbk0.google.com/cbk?output=xml&ll=51.494966,-0.146674 + """ + + NOT_PROVIDED = 'Not provided' + api_header = 'http://cbk0.google.com/cbk?' + api_parameter = 'output=xml' + api_parameter += '&panoid=' + pano_id + api_path = api_header + api_parameter + try: + pano = {'pano_id': pano_id} + do_sleep = False + if os.path.isfile('../data/GSV/' + pano_id + '/meta.xml'): + pano_xml = '../data/GSV/' + pano_id + '/meta.xml' + else: + do_sleep = True + pano_xml = urllib.urlopen(api_path) + tree = ET.parse(pano_xml) + root = tree.getroot() + + for child in root: + if child.tag == 'data_properties': + pano[child.tag] = child.attrib + pano[child.tag]['copyright'] = child.find('copyright').text.strip().replace('\xa9', 'Copyright') + + if child.find('text') != None and child.find('text').text != None: + pano[child.tag]['text'] = child.find('text').text.strip() + else: + pano[child.tag]['text'] = NOT_PROVIDED + + if child.find('street_range') != None and child.find('street_range').text is not None: + pano[child.tag]['street_range'] = child.find('street_range').text.strip() + else: + pano[child.tag]['street_range'] = NOT_PROVIDED + + if child.find('region') is not None and child.find('region').text is not None: + pano[child.tag]['region'] = child.find('region').text.strip() + else: + pano[child.tag]['region'] = NOT_PROVIDED + + if child.find('country') is not None and child.find('country').text is not None: + pano[child.tag]['country'] = child.find('country').text.strip() + else: + pano[child.tag]['country'] = NOT_PROVIDED + elif child.tag == 'projection_properties': + pano[child.tag] = child.attrib + elif child.tag == 'annotation_properties': + pano['links'] = [] + for item in child: + if item.tag == 'link': + link_attrib = {} + link_attrib = item.attrib + + if item.find('link_text') != None and item.find('link_text').text != None: + link_attrib['link_text'] = item.find('link_text').text.strip() + else: + link_attrib['link_text'] = NOT_PROVIDED + + pano['links'].append(link_attrib) + pano['intersection'] = {'lat' : pano['data_properties']['lat'], 'lng' : pano['data_properties']['lng']} + if verbose: + print pano + except: + raise XMLAcquisitionError('Exception: Failed reading xml.') + + if do_sleep: + sleep_in_milliseconds = delay / 1000 + sleep(sleep_in_milliseconds) + return pano + + +def get_nearby_pano_ids(pano_id, max_step_size=2, delay=2000.0, verbose=False): + """ + This function performs breadth first search of GSV panoarama scenes + """ + queue = [{'step_size': 0, 'pano_id': pano_id, 'origin_pano_id': pano_id}] + visited = [] + ret = [] + + while queue: + pano_item = queue.pop(0) + if pano_item['step_size'] > max_step_size: + break + + if pano_item['pano_id'] not in visited: + visited.append(pano_item['pano_id']) + ret.append(pano_item) + pano_data = get_pano_metadata(pano_item['pano_id'], verbose) + linked_pano_ids = [link['pano_id'] for link in pano_data['links']] + for linked_pano_id in linked_pano_ids: + queue.append({'step_size' : pano_item['step_size'] + 1, 'pano_id': linked_pano_id, 'origin_pano_id': pano_id}) + + return ret + + +def get_pano_metadata(pano_id, verbose=False): + """ + This function takes a pano_id (e.g., dWeBDzGMXwQv5fu1GoNy8Q) and + returns a Google Street View panorama metadata that corresponds to the nearest GSV panorama. + E.g., + + http://jamiethompson.co.uk/web/2010/05/15/google-streetview-static-api/ + http://cbk0.google.com/cbk?output=xml&ll=51.494966,-0.146674 + """ + gsv = GSVScraper.GSVScraper() + NOT_PROVIDED = 'Not provided' + + api_header = 'http://cbk0.google.com/cbk?' + api_parameter = 'output=xml' + api_parameter += '&panoid=' + pano_id + api_path = api_header + api_parameter + try: + pano = {'pano_id' : pano_id} + + gsv.get_pano_metadata([pano_id]) + pano_xml = open('../data/GSV/' + pano_id + '/meta.xml', 'rb') + tree = ET.parse(pano_xml) + root = tree.getroot() + + for child in root: + if child.tag == 'data_properties': + pano[child.tag] = child.attrib + pano[child.tag]['copyright'] = child.find('copyright').text.strip().replace('\xa9', 'Copyright') + + if child.find('text') != None and child.find('text').text != None: + pano[child.tag]['text'] = child.find('text').text.strip() + else: + pano[child.tag]['text'] = NOT_PROVIDED + + if child.find('street_range') != None and child.find('street_range').text != None: + pano[child.tag]['street_range'] = child.find('street_range').text.strip() + else: + pano[child.tag]['street_range'] = NOT_PROVIDED + + if child.find('region') != None and child.find('region').text != None: + pano[child.tag]['region'] = child.find('region').text.strip() + else: + pano[child.tag]['region'] = NOT_PROVIDED + + if child.find('country') != None and child.find('country').text != None: + pano[child.tag]['country'] = child.find('country').text.strip() + else: + pano[child.tag]['country'] = NOT_PROVIDED + elif child.tag == 'projection_properties': + pano[child.tag] = child.attrib + elif child.tag == 'annotation_properties': + pano['links'] = [] + for item in child: + if item.tag == 'link': + link_attrib = {} + link_attrib = item.attrib + + if item.find('link_text') != None and item.find('link_text').text != None: + link_attrib['link_text'] = item.find('link_text').text.strip() + else: + link_attrib['link_text'] = NOT_PROVIDED + + pano['links'].append(link_attrib) + pano['bus_stop'] = {'lat' : pano['data_properties']['lat'], 'lng' : pano['data_properties']['lng']} + if verbose: + print pano + except: + raise + + return pano + +def get_nearest_pano_metadata(latlng, delay=1000.0, verbose='True'): + """ + This function takes a latlng object (e.g. {'lat': '38.9015110', 'lng': '-77.0188500'}) and + returns a Google Street View panorama metadata that corresponds to the nearest GSV panorama. + E.g., + + delay: delay in milliseconds + + http://jamiethompson.co.uk/web/2010/05/15/google-streetview-static-api/ + http://cbk0.google.com/cbk?output=xml&ll=51.494966,-0.146674 + """ + NOT_PROVIDED = 'Not provided' + + api_header = 'http://cbk0.google.com/cbk?' + api_parameter = 'output=xml' + api_parameter += '&ll=' + latlng['lat'] + ',' + latlng['lng'] + api_path = api_header + api_parameter + + if verbose: + print api_path + + try: + pano = {'bus_stop': latlng} + pano_xml = urllib.urlopen(api_path) + tree = ET.parse(pano_xml) + root = tree.getroot() + + for child in root: + if child.tag == 'data_properties': + pano[child.tag] = child.attrib + pano[child.tag]['copyright'] = child.find('copyright').text.strip().replace('\xa9', 'Copyright') + + if child.find('text') is not None and child.find('text').text is not None: + pano[child.tag]['text'] = child.find('text').text.strip() + else: + pano[child.tag]['text'] = NOT_PROVIDED + + if child.find('street_range') is not None and child.find('street_range').text is not None: + pano[child.tag]['street_range'] = child.find('street_range').text.strip() + else: + pano[child.tag]['street_range'] = NOT_PROVIDED + + if child.find('region') is not None and child.find('region').text is not None: + pano[child.tag]['region'] = child.find('region').text.strip() + else: + pano[child.tag]['region'] = NOT_PROVIDED + + if child.find('country') is None and child.find('country').text is not None: + pano[child.tag]['country'] = child.find('country').text.strip() + else: + pano[child.tag]['country'] = NOT_PROVIDED + elif child.tag == 'projection_properties': + pano[child.tag] = child.attrib + elif child.tag == 'annotation_properties': + pano['links'] = [] + for item in child: + if item.tag == 'link': + link_attrib = item.attrib + + if item.find('link_text') is not None and item.find('link_text').text is not None: + link_attrib['link_text'] = item.find('link_text').text.strip() + else: + link_attrib['link_text'] = NOT_PROVIDED + + pano['links'].append(link_attrib) + + if verbose: + print pano + except: + raise XMLAcquisitionError('Exception: Failed reading xml.') + + sleep_in_milliseconds = float(delay) / 1000 + sleep(sleep_in_milliseconds) + return pano + + +def collect_all_busstop_depth_data(): + with SidewalkDB(database="busstop") as db: + sql = "select SourceGSVPanoramaId, TargetGSVPanoramaId from PanoramaLink" + records = db.query(sql) + + panorama_ids = [] + for record in records: + panorama_ids.append(record[0]) + panorama_ids.append(record[1]) + + panorama_ids = list(set(panorama_ids)) + + scraper = GSVScraper() + scraper.set_pano_ids(['2V_YrQbwf45Mx9WTJ79ojg']) + scraper.get_pano_metadata() + scraper.get_pano_image() + scraper.get_pano_depthdata() + for panorama_id in panorama_ids: + scraper.set_pano_ids([panorama_id]) + scraper.get_pano_depthdata() + +if __name__ == '__main__': + print "GSVScraper" + collect_all_busstop_depth_data() + + + + + diff --git a/GSVTopDownImage.py b/GSVTopDownImage.py new file mode 100644 index 0000000..2903486 --- /dev/null +++ b/GSVTopDownImage.py @@ -0,0 +1,53 @@ +import GSVImage +import numpy as np +import os + +from PIL import Image + +class GSVTopDownImage(object): + width = 1000 + height = 1000 + def __init__(self, path): + """ + The constructor. This method takes a path to top + """ + if type(path) != str and len(str) > 0: + raise ValueError('path should be a string') + if path[-1] != '/': + path += '/' + + + self.path = path + self.pano_id = self.path.split('/')[-1] + self.filename = 'topdown.png' + self.file_path = self.path + 'images/topdown.png' + if not os.path.isfile(self.file_path): + raise ValueError(self.file_path + ' does not exist.') + + return + + def get_image(self): + """ + This method returns the image array + """ + return np.asarray(Image.open(self.file_path)) + + def show(self, size=False): + """ + This method shows the topdown.png image + """ + if os.path.isfile(self.file_path): + im = Image.open(self.file_path) + + if size and type(size) == tuple: + im = im.resize(size, Image.ANTIALIAS) + im.show() + else: + raise Exception(self.path + 'images/pano.jpg does not exist') + return + +if __name__ == '__main__': + print "GSVTopDownImage.py" + + pano_id = '-015tl-_IqAuhn4X2_km6Q' + gsv_topdown = GSVTopDownImage('../data/GSV/' + pano_id + '/') diff --git a/SidewalkDB.py b/SidewalkDB.py new file mode 100644 index 0000000..8a33118 --- /dev/null +++ b/SidewalkDB.py @@ -0,0 +1,999 @@ +try: + from xml.etree import cElementTree as ET +except ImportError, e: + from xml.etree import ElementTree as ET + +import redis +import pandas.io.sql as pdsql +import pprint +import pymysql +import cPickle as pickle + +from GSVScraper import * +from time import sleep +from pandas import Series, DataFrame + +pp = pprint.PrettyPrinter(indent=2) + + +_sidewalk_db_cache = {} + +class RecordError(Exception): + def __init__(self, value): + self.value = value + + def __str__(self): + return repr(self.value) + + +class XMLAcquisitionError(Exception): + def __init__(self, value): + self.value = value + + def __str__(self): + return repr(self.value) + + +class SidewalkDB(object): + def __init__(self, database='sidewalk', use_sscursor=False, **kwargs): + self._database = database + + if "host" in kwargs: + self._host = kwargs["host"] + else: + self._host = '127.0.0.1' + + if "port" in kwargs: + self._port = kwargs["port"] + else: + self._port = 3306 + + if "user" in kwargs: + self._user = kwargs["user"] + else: + self._user = 'root' + + if "password" in kwargs: + self._password = kwargs["password"] + else: + self._password = '' + if use_sscursor: + self._conn = pymysql.connect(host=self._host, port=self._port, user=self._user, passwd=self._password, + db=self._database, cursorclass=pymysql.SSCursor) + else: + self._conn = pymysql.connect(host=self._host, port=self._port, user=self._user, passwd=self._password, + db=self._database) + self._cur = self._conn.cursor() + + def __enter__(self): + return self + + def __exit__(self, type, value, traceback): + if type is None: + self._cur.connection.commit() + else: + self._cur.connection.rollback() + + self._cur.close() + self._conn.close() + + def fetch_city(self, panorama_id, return_dataframe=False): + sql = """SELECT * FROM Intersections +WHERE NearestGSVPanoramaId = %s""" + return self.query(sql, (panorama_id,), return_dataframe=return_dataframe) + + def fetch_label_points(self, task_description='VerificationExperiment_2', return_dataframe=True, **kwargs): + """ + This method returns label points + """ + sql = """SELECT Labels.LabelGSVPanoramaId, Labels.LabelTypeId, Labels.Deleted, + LabelPoints.LabelId, LabelPoints.svImageX, LabelPoints.svImageY + FROM LabelPoints +INNER JOIN Labels +ON Labels.LabelId = LabelPoints.LabelId +INNER JOIN LabelingTasks +ON LabelingTasks.LabelingTaskId = Labels.LabelingTaskId +INNER JOIN TaskPanoramas +ON TaskPanoramas.TaskPanoramaId = LabelingTasks.TaskPanoramaId +WHERE TaskPanoramas.TaskDescription = %s """ + + if 'only_original' in kwargs and kwargs['only_original']: + sql += " AND LabelingTasks.PreviousLabelingTaskId IS NULL" + + return self.query(sql, (task_description, ), return_dataframe=return_dataframe) + + def fetch_golden_label_points(self, label_type_id, return_dataframe=True): + """ + This method retrieves the golden label points + """ + sql = """SELECT LabelPoints.LabelId, svImageX, svImageY, Labels.LabelGSVPanoramaId, Assignments.*, + LabelingTasks.PreviousLabelingTaskId, LabelingTasks.LabelingTaskId +FROM Assignments +INNER JOIN LabelingTasks +ON LabelingTasks.AssignmentId = Assignments.AssignmentId +INNER JOIN TaskPanoramas +ON TaskPanoramas.TaskGSVPanoramaId = LabelingTasks.TaskGSVPanoramaId +INNER JOIN Labels +ON Labels.LabelingTaskId = LabelingTasks.LabelingTaskId +INNER JOIN LabelPoints +ON LabelPoints.LabelId = Labels.LabelId +WHERE +Assignments.AmazonTurkerId = 'Researcher_Kotaro' AND +TaskPanoramas.TaskDescription = 'GoldenInsertion' AND +Labels.Deleted = 0 AND +Labels.LabelTypeId = %s""" + df = self.query(sql, (str(label_type_id)), return_dataframe=return_dataframe) + previous_task_ids = df.PreviousLabelingTaskId[~df.PreviousLabelingTaskId.isnull()].unique() + df = df[~df.LabelingTaskId.isin(previous_task_ids)] + return df + + def fetch_labeling_tasks(self, task_description, panorama_id=None, assignment_description=None, workers=None, + return_dataframe=False): + sql = """ + SELECT Assignments.AmazonTurkerId, Assignments.TaskDescription AS AssignmentDescription, + LabelingTasks.LabelingTaskId, LabelingTasks.AssignmentId, LabelingTasks.TaskPanoramaId, + LabelingTasks.TaskGSVPanoramaId AS GSVPanoramaId, LabelingTasks.NoLabel, LabelingTasks.Description, + LabelingTasks.PreviousLabelingTaskId, TaskPanoramas.TaskDescription, Intersections.City + FROM Assignments + INNER JOIN LabelingTasks + ON LabelingTasks.AssignmentId = Assignments.AssignmentId + INNER JOIN TaskPanoramas + ON TaskPanoramas.TaskPanoramaId = LabelingTasks.TaskPanoramaId + INNER JOIN Intersections + ON Intersections.NearestGSVPanoramaId = LabelingTasks.TaskGSVPanoramaId + WHERE TaskPanoramas.TaskDescription = %s + AND Assignments.Completed = 1 + """ + + if panorama_id is not None: + sql += " AND TaskPanoramas.TaskGSVPanoramaId = '%s'" % panorama_id + + if assignment_description is not None: + sql += " AND Assignments.TaskDescription = '%s'" % assignment_description + + if workers is not None: + workers = ["'%s'" % worker for worker in workers] + sql += " AND (Assignments.AmazonTurkerId = " + \ + " OR Assignments.AmazonTurkerId = ".join(workers) + ")" + + return self.query(sql, (task_description,), return_dataframe=return_dataframe) + + def fetch_labeling_tasks_with_task_ids(self, task_ids, return_dataframe=True): + sql = """ + SELECT Assignments.AmazonTurkerId, Assignments.TaskDescription AS AssignmentDescription, + LabelingTasks.LabelingTaskId, LabelingTasks.AssignmentId, LabelingTasks.TaskPanoramaId, + LabelingTasks.TaskGSVPanoramaId AS GSVPanoramaId, LabelingTasks.NoLabel, LabelingTasks.Description, + LabelingTasks.PreviousLabelingTaskId, TaskPanoramas.TaskDescription, Intersections.City + FROM Assignments + INNER JOIN LabelingTasks + ON LabelingTasks.AssignmentId = Assignments.AssignmentId + INNER JOIN TaskPanoramas + ON TaskPanoramas.TaskPanoramaId = LabelingTasks.TaskPanoramaId + INNER JOIN Intersections + ON Intersections.NearestGSVPanoramaId = LabelingTasks.TaskGSVPanoramaId + WHERE Assignments.Completed = 1 + """ + task_ids = map(str, task_ids) + sql += " AND (LabelingTasks.LabelingTaskId = " + \ + " OR LabelingTasks.LabelingTaskId = ".join(task_ids) + ")" + + return self.query(sql, return_dataframe=return_dataframe) + + def fetch_labels(self, label_type_id=1, workers=["Researcher_Kotaro"], return_dataframe=True): + sql = """ + SELECT Labels.*, LabelingTasks.PreviousLabelingTaskId, Intersections.City from Labels + INNER JOIN LabelingTasks + ON LabelingTasks.LabelingTaskId = Labels.LabelingTaskId + INNER JOIN Assignments + ON Assignments.AssignmentId = LabelingTasks.AssignmentId + INNER JOIN Intersections + ON Intersections.NearestGSVPanoramaId = LabelingTasks.TaskGSVPanoramaId + """ + if workers is not None: + workers = ["'%s'" % worker for worker in workers] + sql += " WHERE (Assignments.AmazonTurkerId = " + \ + " OR Assignments.AmazonTurkerId = ".join(workers) + ")" + + df = self.query(sql, return_dataframe=return_dataframe) + df = df[df.Deleted != 1] + previous_task_ids = df.PreviousLabelingTaskId[~df.PreviousLabelingTaskId.isnull()].unique() + df = df[~df.LabelingTaskId.isin(previous_task_ids)] + + return df + + def fetch_task_label_points(self, labeling_task_id, return_dataframe=True, **kwargs): + if 'task_description' in kwargs: + sql = """SELECT LabelingTasks.LabelingTaskId, LabelingTasks.TaskGSVPanoramaId AS GSVPanoramaId, PanoYawDeg, LabelType, +Labels.LabelId, Labels.Deleted, LabelPoints.LabelPointId, svImageX, svImageY, heading, pitch, zoom, labelLat, labelLng +FROM LabelingTasks +INNER JOIN TaskPanoramas +ON TaskPanoramas.TaskPanoramaId = LabelingTasks.TaskPanoramaId +INNER JOIN PanoramaProjectionProperties +ON PanoramaProjectionProperties.GSVPanoramaId = LabelingTasks.TaskGSVPanoramaId +INNER JOIN Labels +ON Labels.LabelingTaskId = LabelingTasks.LabelingTaskId +INNER JOIN LabelTypes +ON Labels.LabelTypeId = LabelTypes.LabelTypeId +INNER JOIN LabelPoints +ON LabelPoints.LabelId = Labels.LabelId +WHERE TaskPanoramas.TaskDescription = %s""" + + global _sidewalk_db_cache + key = 'SidewalkDB.fetch_task_label_points.%s' % kwargs['task_description'] + + if key in _sidewalk_db_cache: + df = _sidewalk_db_cache[key] + else: + df = self.query(sql, (kwargs['task_description'],), return_dataframe=True) + _sidewalk_db_cache[key] = df + return df[df.LabelingTaskId == labeling_task_id] + else: + sql = """SELECT LabelingTasks.LabelingTaskId, LabelingTasks.TaskGSVPanoramaId AS GSVPanoramaId, PanoYawDeg, LabelType, +Labels.LabelId, Labels.Deleted, LabelPoints.LabelPointId, svImageX, svImageY, heading, pitch, zoom, labelLat, labelLng +FROM LabelingTasks +INNER JOIN PanoramaProjectionProperties +ON PanoramaProjectionProperties.GSVPanoramaId = LabelingTasks.TaskGSVPanoramaId +INNER JOIN Labels +ON Labels.LabelingTaskId = LabelingTasks.LabelingTaskId +INNER JOIN LabelTypes +ON Labels.LabelTypeId = LabelTypes.LabelTypeId +INNER JOIN LabelPoints +ON LabelPoints.LabelId = Labels.LabelId +WHERE LabelingTasks.LabelingTaskId = %s""" + return self.query(sql, (str(labeling_task_id),), return_dataframe=return_dataframe) + + def fetch_image_level_labels(self, task_description='UISTTurkTasks', assignment_group='TurkerTask', header=False, + return_dataframe=False, **kwargs): + """ + This method fetches + """ + sql = """ + SELECT Assignments.AmazonTurkerId, LabelingTasks.NoLabel, LabelingTasks.Description, + LabelingTasks.PreviousLabelingTaskId, Labels.LabelId, Labels.LabelingTaskId, LabelingTasks.TaskGSVPanoramaId, + Labels.LabelTypeId, Labels.Deleted FROM Assignments + INNER JOIN LabelingTasks + ON LabelingTasks.AssignmentId = Assignments.AssignmentId + INNER JOIN TaskPanoramas + ON TaskPanoramas.TaskPanoramaId = LabelingTasks.TaskPanoramaId + LEFT JOIN Labels + ON Labels.LabelingTaskId = LabelingTasks.LabelingTaskId + WHERE Assignments.Completed = 1 + AND Assignments.TaskDescription = %s + AND TaskPanoramas.TaskDescription = %s + AND (Labels.Deleted = 0 OR Labels.Deleted IS NULL) + """ + + if 'workers' in kwargs and (type(kwargs['workers']) == tuple or type(kwargs['workers']) == list) \ + and len(kwargs['workers']) > 0: + workers = ["'%s'" % worker for worker in kwargs['workers']] + sql += " AND (Assignments.AmazonTurkerId = " + \ + " OR Assignments.AmazonTurkerId = ".join(workers) + ")" + + if return_dataframe: + return self.query(sql, (assignment_group, task_description), return_dataframe=return_dataframe) + else: + if header: + header_row = ("AmazonTurkerId", "NoLabel", "Description", "PreviousLabelingTaskId", "LabelId", + "LabelingTaskId", "GSVPanoramaId", "LabelTypeId", "Deleted") + return [header_row] + list(self.query(sql, (assignment_group, task_description))) + + else: + return self.query(sql, (assignment_group, task_description)) + + def fetch_intersections(self, return_dataframe=False): + """ + This method fetches the intersections + """ + sql = "SELECT * FROM Intersections" + return self.query(sql, return_dataframe=return_dataframe) + + def fetch_metadata(self, pano_id): + """ + This method fetches the panorama metadata of the specified panorama + """ + sql = "SELECT * FROM Panoramas WHERE GSVPanoramaId = %s" + return self.query(sql, (pano_id,)) + + def fetch_outlines(self, task_description='PilotTask_v2_MountPleasant', header=False, + neglect_deleted=False, limit=False, return_dataframe=False, **kwargs): + """ + This method fetches the outlines provided by labelers + !!! Use the helper function helper_clean_up_outlines(). Because the records fetched with + this method include the outlines that are modified and outdated!!! + """ + sql = """ + SELECT LabelingTasks.TaskGSVPanoramaId, Assignments.AmazonTurkerId, PanoYawDeg, LabelType, LabelPointId, + Labels.LabelId, svImageX, svImageY, heading, pitch, zoom, labelLat, labelLng, Intersections.City, + LabelingTasks.LabelingTaskId, LabelingTasks.PreviousLabelingTaskId, Labels.Deleted + FROM LabelingTasks + INNER JOIN Assignments + ON Assignments.AssignmentId = LabelingTasks.AssignmentId + INNER JOIN PanoramaProjectionProperties + ON PanoramaProjectionProperties.GSVPanoramaId = LabelingTasks.TaskGSVPanoramaId + INNER JOIN Labels + ON Labels.LabelingTaskId = LabelingTasks.LabelingTaskId + INNER JOIN LabelTypes + ON Labels.LabelTypeId = LabelTypes.LabelTypeId + INNER JOIN LabelPoints + ON Labels.LabelId = LabelPoints.LabelId + INNER JOIN TaskPanoramas + ON LabelingTasks.TaskPanoramaId = TaskPanoramas.TaskPanoramaId + INNER JOIN Intersections + ON Intersections.NearestGSVPanoramaId = LabelingTasks.TaskGSVPanoramaId + WHERE TaskPanoramas.TaskDescription = %s + AND Assignments.Completed = 1 + """ + + if 'workers' in kwargs and (type(kwargs['workers']) == tuple or type(kwargs['workers']) == list) \ + and len(kwargs['workers']) > 0: + workers = ["'%s'" % worker for worker in kwargs['workers']] + sql += " AND (Assignments.AmazonTurkerId = " + \ + " OR Assignments.AmazonTurkerId = ".join(workers) + ")" + + if 'pano_ids' in kwargs and (type(kwargs['pano_ids']) == set or type(kwargs['pano_ids']) == list or + type(kwargs['pano_ids']) == tuple) and len(kwargs['pano_ids']) > 0: + pano_ids = ["'%s'" % pano_id for pano_id in kwargs['pano_ids']] + sql += " AND (LabelingTasks.TaskGSVPanoramaId = " + \ + " OR LabelingTasks.TaskGSVPanoramaId = ".join(pano_ids) + ")" + + if neglect_deleted: + sql += " AND Labels.Deleted = 0" + + if 'only_original' in kwargs and kwargs['only_original']: + sql += " AND LabelingTasks.PreviousLabelingTaskId IS NULL" + + if limit: + sql += " LIMIT %s" % (limit) + + if return_dataframe: + return self.query(sql, (task_description,), return_dataframe=return_dataframe) + else: + if header: + header_row = ("TaskGSVPanoramaId", "AmazonTurkerId", "PanoYawDeg", "LabelType", "LabelPointId", + "LabelId", "svImageX", "svImageY", "heading", "pitch", "zoom", "labelLat", "labelLng", + "City", "LabelingTaskId", "PreviousLabelingTaskId", "Deleted") + records = [header_row] + list(self.query(sql, (task_description,))) + return records + else: + #return self.query(sql, (task_description,), lazy=True) + return self.yield_query(sql, (task_description,)) + + def fetch_quick_checks(self, task_description="UISTTurkTasks_2", assignment_group='TurkerTask', + return_dataframe=False): + sql = """ + SELECT LabelingTaskInteractions.*, Assignments.AmazonTurkerId FROM LabelingTaskInteractions + INNER JOIN LabelingTasks + ON LabelingTasks.LabelingTaskId = LabelingTaskInteractions.LabelingTaskId + INNER JOIN Assignments + ON Assignments.AssignmentId = LabelingTasks.AssignmentId + WHERE Action = %s + AND Assignments.TaskDescription = %s + """ + return self.query(sql, ("OnboardingQuickCheck_submitClick", assignment_group,), + return_dataframe=return_dataframe) + + def fetch_quick_verification_results(self, task_description="VerificationImages_v2", assignment_group="TurkerTask", + return_dataframe=True): + """ + This method returns the results of QuickVerification tasks. + """ + sql = """SELECT Labels.LabelId, Labels.LabelGSVPanoramaId, +Assignments.AssignmentId, Assignments.AmazonTurkerId, +ValidationTaskResults.ValidationTaskResultId, ValidationTaskResults.LabelTypeId AS VerifiedLabelTypeId, +ValidationTaskResults.Timestamp +FROM TaskImages +INNER JOIN ValidationTaskResults +ON ValidationTaskResults.TaskImageId = TaskImages.TaskImageId +INNER JOIN ValidationTasks +ON ValidationTasks.ValidationTaskId = ValidationTaskResults.ValidationTaskId +INNER JOIN Assignments +ON Assignments.AssignmentId = ValidationTasks.AssignmentId +INNER JOIN Images +ON Images.ImageId = TaskImages.ImageId +INNER JOIN Labels +ON Labels.LabelId = Images.LabelId +WHERE Assignments.TaskDescription = %s AND +Assignments.Completed = 1 AND +TaskImages.TaskDescription = %s""" + return self.query(sql, (assignment_group, task_description), return_dataframe=return_dataframe) + + def fetch_skipped_tasks(self, task_description='UISTTurkTasks', assignment_group='TurkerTask', header=False, + return_dataframe=False, **kwargs): + """ + + """ + sql = """ + SELECT LabelingTasks.TaskGSVPanoramaId, Assignments.AmazonTurkerId, PanoYawDeg, Intersections.City, + LabelingTasks.LabelingTaskId, LabelingTasks.PreviousLabelingTaskId, LabelingTasks.NoLabel + FROM LabelingTasks + INNER JOIN Assignments + ON Assignments.AssignmentId = LabelingTasks.AssignmentId + INNER JOIN PanoramaProjectionProperties + ON PanoramaProjectionProperties.GSVPanoramaId = LabelingTasks.TaskGSVPanoramaId + INNER JOIN TaskPanoramas + ON LabelingTasks.TaskPanoramaId = TaskPanoramas.TaskPanoramaId + INNER JOIN Intersections + ON Intersections.NearestGSVPanoramaId = LabelingTasks.TaskGSVPanoramaId + WHERE LabelingTasks.NoLabel = 1 + AND Assignments.TaskDescription = %s + AND TaskPanoramas.TaskDescription = %s + """ + + if 'workers' in kwargs and (type(kwargs['workers']) == tuple or type(kwargs['workers']) == list) and \ + len(kwargs['workers']) > 0: + workers = ["'%s'" % worker for worker in kwargs['workers']] + sql += " AND (Assignments.AmazonTurkerId = " + \ + " OR Assignments.AmazonTurkerId = ".join(workers) + ")" + + if 'pano_ids' in kwargs and (type(kwargs['pano_ids']) == set or type(kwargs['pano_ids']) == list or + type(kwargs['pano_ids']) == tuple) and len(kwargs['pano_ids']) > 0: + pano_ids = ["'%s'" % pano_id for pano_id in kwargs['pano_ids']] + sql += " AND (LabelingTasks.TaskGSVPanoramaId = " + \ + " OR LabelingTasks.TaskGSVPanoramaId = ".join(pano_ids) + ")" + + if return_dataframe: + return self.query(sql, (assignment_group, task_description,), return_dataframe=return_dataframe) + else: + if header: + header_row = ("TaskGSVPanoramaId", "AmazonTurkerId", "PanoYawDeg", "City", "LabelingTaskId", + "PreviousLabelingTaskId", "NoLabel") + records = [header_row] + list(self.query(sql, (assignment_group, task_description,))) + return records + else: + return self.yield_query(sql, (assignment_group, task_description,)) + + def fetch_task_interaction(self, task_description='UISTTurkTasks', assignment_group='TurkerTask', header=False, + details=True, return_dataframe=True, **kwargs): + """ + This method returns all the task interactions + + :param kwargs: + Pass the types of interactions that you want to retrieve. If a key 'actions' is in kwargs + (i.e., kwargs['actions'] ) and the value is tuple. + """ + + # Cache + # http://stackoverflow.com/questions/146557/do-you-use-the-global-statement-in-python + + if details: + # provide details of the task + sql = """ + SELECT LabelingTaskInteractions.*, Assignments.*, LabelingTasks.PreviousLabelingTaskId, + LabelingTasks.NoLabel + FROM LabelingTaskInteractions + INNER JOIN LabelingTasks + ON LabelingTasks.LabelingTaskId = LabelingTaskInteractions.LabelingTaskId + INNER JOIN TaskPanoramas + ON TaskPanoramas.TaskPanoramaId = LabelingTasks.TaskPanoramaId + INNER JOIN Assignments + ON Assignments.AssignmentId = LabelingTasks.AssignmentId + WHERE Assignments.Completed = 1 + AND Assignments.TaskDescription = %s + AND TaskPanoramas.TaskDescription = %s + """ + else: + sql = """ + SELECT LabelingTaskInteractions.* FROM LabelingTaskInteractions + INNER JOIN LabelingTasks + ON LabelingTasks.LabelingTaskId = LabelingTaskInteractions.LabelingTaskId + INNER JOIN TaskPanoramas + ON TaskPanoramas.TaskPanoramaId = LabelingTasks.TaskPanoramaId + INNER JOIN Assignments + ON Assignments.AssignmentId = LabelingTasks.AssignmentId + WHERE Assignments.Completed = 1 + AND Assignments.TaskDescription = %s + AND TaskPanoramas.TaskDescription = %s + """ + + if 'actions' in kwargs and (type(kwargs['actions']) == tuple or type(kwargs['actions']) == list) and \ + len(kwargs['actions']) > 0: + actions = ["'%s'" % action for action in kwargs['actions']] + sql += " AND (LabelingTaskInteractions.Action = " + \ + " OR LabelingTaskInteractions.Action = ".join(actions) + ")" + if 'pano_id' in kwargs: + sql += " AND LabelingTaskInteractions.GSVPanoramaId = '%s'" % kwargs['pano_id'] + + if return_dataframe: + df = self.query(sql, (assignment_group, task_description), return_dataframe=return_dataframe) + + # Cache + # http://stackoverflow.com/questions/146557/do-you-use-the-global-statement-in-python + return df + else: + if header: + if details: + header_row = ("LabelingTaskInteractionId", "LabelingTaskId", "Action", "GSVPanoramaId", "Lat", "Lng", + "Heading", "Pitch", "Zoom", "Note", "Timestamp", "AssignmentId", "AmazonTurkerId", "AmazonHitId", + "AmazonAssignmentId", "InterfaceType", "InterfaceVersion", "Completed", "NeedQualification", + "TaskDescription", "DatetimeInserted") + else: + header_row = ("LabelingTaskInteractionId", "LabelingTaskId", "Action", "GSVPanoramaId", "Lat", "Lng", + "Heading", "Pitch", "Zoom", "Note", "Timestamp") + return [header_row] + list(self.query(sql, (assignment_group, task_description), + return_dataframe=return_dataframe)) + + else: + + df = self.query(sql, (assignment_group, task_description), return_dataframe=return_dataframe) + + return df + + def fetch_task_interactions_by_task_id(self, labeling_task_id, return_dataframe=False, **kwargs): + if 'task_description' in kwargs: + df = self.fetch_task_interaction(kwargs['task_description']) + return df[df.LabelingTaskId == labeling_task_id] + else: + sql = """SELECT * FROM LabelingTaskInteractions WHERE LabelingTaskInteractions.LabelingTaskId = %s""" + return self.query(sql, (str(labeling_task_id),), return_dataframe=return_dataframe) + + def fetch_labeling_task_attribute(self, labeling_task_id, attribute): + global _sidewalk_db_cache + key = 'SidewalkDB.fetch_labeling_task_attribute' + + if key in _sidewalk_db_cache: + df = _sidewalk_db_cache[key] + else: + sql = """SELECT LabelingTaskId, Attribute, Value FROM LabelingTaskAttributes""" + df = self.query(sql, return_dataframe=True) + _sidewalk_db_cache[key] = df + + if len(df[(df.LabelingTaskId == labeling_task_id) & (df.Attribute == attribute)]) > 0: + return df[(df.LabelingTaskId == 18565) & (df.Attribute == 'FalseNegative')].Value.iloc[0] + else: + return None + + + sql = """SELECT Attribute, Value FROM LabelingTaskAttributes + WHERE LabelingTaskId = %s AND Attribute = %s""" + records = self.query(sql, (str(labeling_task_id), attribute)) + if len(records) > 0: + return records[0][1] + else: + return None + + + + def fetch_task_counts(self, task_description='PilotTask_v2_MountPleasant', return_dataframe=False, **kwargs): + """ + + """ + sql = """SELECT * FROM LabelingTaskCounts + INNER JOIN TaskPanoramas + ON TaskPanoramas.TaskPanoramaId = LabelingTaskCounts.TaskPanoramaId + WHERE TaskDescription = %s""" + return self.query(sql, (task_description,), return_dataframe=return_dataframe) + + def fetch_task_panorama(self, panorama_id, task_description, return_dataframe=False, **kwargs): + """ + This method returns the specified task panorama + """ + sql = """SELECT * FROM TaskPanoramas + WHERE TaskPanoramas.TaskGSVPanoramaId = %s + AND TaskPanoramas.TaskDescription = %s""" + return self.query(sql, (panorama_id, task_description), return_dataframe=return_dataframe) + + def fetch_task_panoramas(self, task_description='PilotTask_v2_MountPleasant', header=False, **kwargs): + """ + This method returns the all the task panoramas + + :param task_description: + The task group name. I've been using PilotTask_v2_MountPleasant although the area is not limited to + Mt. Pleasant and is misleading. Because I am lazy to update all the code accordingly + :param header: + If header is true, it inserts the header too. + :return: + Task panoramas + """ + + if "with_id" in kwargs and kwargs["with_id"]: + if "VerificationImages" in task_description: + sql = """SELECT TaskPanoramas.TaskPanoramaId, TaskPanoramas.TaskGSVPanoramaId, Intersections.City, + TaskPanoramas.TaskDescription FROM TaskImages + INNER JOIN Images + ON Images.ImageId = TaskImages.ImageId + INNER JOIN Labels + ON Labels.LabelId = Images.LabelId + INNER JOIN LabelingTasks + ON LabelingTasks.LabelingTaskId = Labels.LabelingTaskId + INNER JOIN TaskPanoramas + ON TaskPanoramas.TaskPanoramaId = LabelingTasks.TaskPanoramaId + INNER JOIN Intersections + ON Intersections.NearestGSVPanoramaId = TaskPanoramas.TaskGSVPanoramaId + WHERE TaskImages.TaskDescription = %s + """ + else: + sql = """ + SELECT TaskPanoramas.TaskPanoramaId, TaskPanoramas.TaskGSVPanoramaId, Intersections.City, + TaskPanoramas.TaskDescription FROM TaskPanoramas + INNER JOIN Intersections + ON Intersections.NearestGSVPanoramaId = TaskPanoramas.TaskGSVPanoramaId + WHERE TaskDescription = %s + """ + if header: + header_row = ('TaskPanoramaId', 'GSVPanoramaId', 'City', 'TaskDescription') + records = [header_row] + list(self.query(sql, (task_description,))) + return records + else: + return self.query(sql, (task_description,)) + else: + if "VerificationImages" in task_description: + sql = """SELECT TaskPanoramas.TaskGSVPanoramaId, Intersections.City, TaskPanoramas.TaskDescription FROM TaskImages + INNER JOIN Images + ON Images.ImageId = TaskImages.ImageId + INNER JOIN Labels + ON Labels.LabelId = Images.LabelId + INNER JOIN LabelingTasks + ON LabelingTasks.LabelingTaskId = Labels.LabelingTaskId + INNER JOIN TaskPanoramas + ON TaskPanoramas.TaskPanoramaId = LabelingTasks.TaskPanoramaId + INNER JOIN Intersections + ON Intersections.NearestGSVPanoramaId = TaskPanoramas.TaskGSVPanoramaId + WHERE TaskImages.TaskDescription = %s + """ + else: + sql = """ + SELECT TaskPanoramas.TaskGSVPanoramaId, Intersections.City, TaskPanoramas.TaskDescription + FROM TaskPanoramas + INNER JOIN Intersections + ON Intersections.NearestGSVPanoramaId = TaskPanoramas.TaskGSVPanoramaId + WHERE TaskDescription = %s + """ + if header: + header_row = ('GSVPanoramaId', 'City', 'TaskDescription') + records = [header_row] + list(self.query(sql, (task_description,))) + return records + else: + return self.query(sql, (task_description,)) + + def fetch_tasks(self, task_description='PilotTask_v2_MountPleasant', return_dataframe=False, **kwargs): + """ + + """ + sql = """SELECT * FROM TaskPanoramas WHERE TaskDescription = %s""" + return self.query(sql, (task_description,), return_dataframe=return_dataframe) + + def get_last_row_id(self): + """ + This method returns the lastrowid + """ + return self._cur.lastrowid + + def insert_assignment(self, amazon_worker_id, amazon_hit_id='DefaultValue', amazon_assignment_id='DefaultValue', + interface_type='', interface_version='', completed='0', need_qualification='0', + assignment_description=''): + sql = "SELECT * FROM Assignments WHERE AmazonTurkerId = %s" + records = self.query(sql, (amazon_worker_id)) + if len(records) > 0: + assignment_id = records[0][0] + return assignment_id + else: + sql = """INSERT Assignments (AmazonTurkerId, AmazonHitId, AmazonAssignmentId, InterfaceType, + InterfaceVersion, Completed, NeedQualification, TaskDescription) + VALUES (%s, %s, %s, %s, %s, %s, %s, %s) + """ + self.query(sql, (amazon_worker_id, amazon_hit_id, amazon_assignment_id, interface_type, interface_version, + completed, need_qualification, assignment_description)) + return None + + def insert_detected_bounding_box(self, record, bin_description): + """ + This method takes a record + + :param record: + Record should have {'assignments': [...], 'labeling_tasks': [...], 'labels': [...], 'label_points': [...]} + :param bin_description: + Description of the detected bounding boxes. + """ + sql = """SELECT * FROM Assignments WHERE AmazonTurkerId=%s""" + for asg in record['assignments']: + amazon_turker_id = str(asg[0]) + cur_data = self.query(sql, (amazon_turker_id,)) + + if len(cur_data) == 0: + sql = """ + INSERT INTO Assignments + (AmazonTurkerId, AmazonHitId, AmazonAssignmentId, + InterfaceType, InterfaceVersion, Completed, TaskDescription) + VALUES (%s, %s, %s, %s, %s, %s, %s) + """ + self._cur.execute(sql, (asg)) + asg_id = self._cur.lastrowid + else: + asg_id = cur_data[0][0] + + sql = """SELECT * FROM LabelingTasks WHERE AssignmentId=%s AND TaskGSVPanoramaId=%s""" + for task in record['labeling_tasks']: + task_gsv_panorama_id = task[0] + records = self.query(sql, (asg_id, task_gsv_panorama_id)) + if len(records) == 0: + sql = """ + INSERT INTO LabelingTasks (AssignmentId, TaskGSVPanoramaId, NoLabel, Description) + VALUES (%s, %s, %s, %s) + """ + self._cur.execute(sql, (asg_id, task[0], task[1], task[2])) + task_id = self._cur.lastrowid + else: + task_id = records[0][0] + + # + # Go through all the labels and points and check if the exact same points exist or not. + # If so, do not enter results. + for label in record['labels']: + sql = """ + INSERT INTO Labels (LabelingTaskid, LabelGSVPanoramaId, LabelTypeId, Deleted) + VALUES (%s, %s, %s, %s) + """ + cur_data = self._cur.execute(sql, (task_id, label[0], label[1], label[2])) + label_id = self._cur.lastrowid + + sql = """ + INSERT INTO LabelPoints (LabelId, svImageX, svImageY, heading, pitch, zoom, labelLat, labelLng) + VALUES (%s, %s, %s, %s, %s, %s, %s, %s) + """ + for label_point in record['label_points']: + cur_data = self._cur.execute(sql, (label_id, label_point[0], label_point[1], label_point[2], label_point[3], + label_point[4], label_point[5], label_point[6])) + + sql = """INSERT INTO BinnedLabels (LabelBinId, LabelId) VALUES (%s, %s)""" + return + + def insert_intersections(self, data, force=False): + """ + This function takes data returned by get_nearest_pano_metadata() helper function + If force is True, insert the bus stop data even if the bus sotp already exists. + """ + # Check if the bus stop record already exist + if not force: + sql = """SELECT * FROM Intersections WHERE NearestGSVPanoramaId=%s""" + cur_data = self._cur.execute(sql, (data['data_properties']['pano_id'])) + if cur_data != 0: + raise RecordError('The record already exists.') + + num_links = str(len(data['links'])) + sql = """INSERT INTO Intersections (Lat, Lng, NearestGSVPanoramaId, NumberOfLinks) VALUES (%s, %s, %s, %s)""" + self._cur.execute(sql, (data['intersection']['lat'], data['intersection']['lng'], + data['data_properties']['pano_id'], num_links)) + return + + def insert_label(self, labeling_task_id, gsv_panorama_id, label_type_id, deleted): + """ + This method inserts a label into the Labels table + """ + sql = """INSERT Labels (LabelingTaskId, LabelGSVPanoramaId, LabelTypeId, Deleted) + VALUES (%s, %s, %s, %s)""" + self.query(sql, (labeling_task_id, gsv_panorama_id, label_type_id, deleted)) + return + + def insert_label_point(self, label_id, x, y, heading, pitch, zoom, lat, lng): + """ + This method inserts a label point into the LabelPoints table + """ + sql = """INSERT INTO LabelPoints (LabelId, svImageX, svImageY, heading, pitch, zoom, labelLat, labelLng) + VALUES (%s, %s, %s, %s, %s, %s, %s, %s)""" + self.query(sql, (label_id, x, y, heading, pitch, zoom, lat, lng)) + return + + def insert_label_confidence(self, label_id, confidence): + """ + This method inserts a label confidence record into the LabelConfidenceScores table + """ + sql = "INSERT INTO LabelConfidenceScores (LabelId, ConfidenceScore) VALUES (%s, %s)" + self.query(sql, (label_id, confidence)) + return + + def insert_labeling_task(self, assignment_id, task_panorama_id, gsv_panorama_id, no_label, description): + """ + This method inserts a LabelingTask record + """ + sql = "SELECT LabelingTaskId FROM LabelingTasks WHERE TaskPanoramaId = %s" + records = self.query(sql, (str(task_panorama_id))) + if len(records) > 0: + labeling_task_id = records[0][0] + return labeling_task_id + else: + sql = """INSERT LabelingTasks (AssignmentId, TaskPanoramaId, TaskGSVPanoramaId, NoLabel, Description) + VALUES (%s, %s, %s, %s, %s)""" + + self.query(sql, (str(assignment_id), str(task_panorama_id), gsv_panorama_id, str(no_label), description)) + labeling_task_id = self.get_last_row_id() + return None + + def insert_nearby_panorama(self, data, force=False): + """ + This function takes one item from the data returned by get_nearby_pano_ids() helper function + If force is True, insert the bus stop data even if the bus sotp already exists. + """ + # Check if the bus stop record already exist + if not force: + sql = """SELECT * FROM NearbyPanoramas WHERE TaskGSVPanoramaId=%s AND GSVPanoramaId=%s""" + cur_data = self._cur.execute(sql, (data['origin_pano_id'], data['pano_id'])) + if cur_data != 0: + raise RecordError('The record already exists.') + + sql = """INSERT INTO NearbyPanoramas (TaskGSVPanoramaId, GSVPanoramaId, StepSize) VALUES (%s, %s, %s)""" + self._cur.execute(sql, (data['origin_pano_id'], data['pano_id'], data['step_size'])) + return + + def insert_panoramas(self, data): + """ + This function takes data returned by get_nearest_pano_metadata() helper function + + :param data: An object that holds panoarama data + """ + # Check if the record already exists + sql = """SELECT * FROM Panoramas WHERE GSVPanoramaId=%s""" + cur_data = self._cur.execute(sql, (data['data_properties']['pano_id'])) + if cur_data != 0: + raise RecordError('The record already exists.') + + sql = """INSERT INTO Panoramas ( + GSVPanoramaId, + ImageWidth, + ImageHeight, + TileWidth, + TileHeight, + ImageDate, + NumZoomLevels, + Lat, + Lng, + OriginalLat, + OriginalLng, + ElevationWgs84M, + Copyright, + Text, + StreetRange, + Region, + Country) VALUES (%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s)""" + + # + if 'elevation_wgs84_m' not in data['data_properties']: + data['data_properties']['elevation_wgs84_m'] = 'undefined' + if 'image_date' not in data['data_properties']: + data['data_properties']['image_date'] = 'undefined' + print 'image_date undefined' + + self._cur.execute(sql, (str(data['data_properties']['pano_id']), + str(data['data_properties']['image_width']), + str(data['data_properties']['image_height']), + str(data['data_properties']['tile_width']), + str(data['data_properties']['tile_height']), + str(data['data_properties']['image_date']), + str(data['data_properties']['num_zoom_levels']), + str(data['data_properties']['lat']), + str(data['data_properties']['lng']), + str(data['data_properties']['original_lat']), + str(data['data_properties']['original_lng']), + str(data['data_properties']['elevation_wgs84_m']), + str(data['data_properties']['copyright']), + str(data['data_properties']['text']), + str(data['data_properties']['street_range']), + str(data['data_properties']['region']), + str(data['data_properties']['country']))) + return + + def insert_panorama_projection_properties(self, data): + """ + This function takes data returned by get_nearest_pano_metadata() helper function + """ + # Check if the record already exists + sql = """SELECT * FROM PanoramaProjectionProperties WHERE GSVPanoramaId=%s""" + cur_data = self._cur.execute(sql, (data['data_properties']['pano_id'])) + if cur_data != 0: + raise RecordError('The record already exists.') + + sql = """INSERT INTO PanoramaProjectionProperties ( + GSVPanoramaId, + ProjectionType, + PanoYawDeg, + TiltYawDeg, + TiltPitchDeg) VALUES (%s,%s,%s,%s,%s)""" + + self._cur.execute(sql, (str(data['data_properties']['pano_id']), + str(data['projection_properties']['projection_type']), + str(data['projection_properties']['pano_yaw_deg']), + str(data['projection_properties']['tilt_yaw_deg']), + str(data['projection_properties']['tilt_pitch_deg']))) + return + + def insert_panorama_links(self, data): + """ + This function takes data returned by get_nearest_pano_metadata() helper function + """ + sql = """SELECT * FROM PanoramaLinks WHERE SourceGSVPanoramaId=%s""" + cur_data = self._cur.execute(sql, (data['data_properties']['pano_id'])) + if cur_data != 0: + raise RecordError('The record already exists.') + + sql = """INSERT INTO PanoramaLinks ( + SourceGSVPanoramaId, + TargetGSVPanoramaId, + YawDeg, + RoadArgb, + Scene, + LinkText) VALUES (%s,%s,%s,%s,%s,%s)""" + + for link in data['links']: + self._cur.execute(sql, (str(data['data_properties']['pano_id']), + link['pano_id'], + link['yaw_deg'], + link['road_argb'], + link['scene'], + link['link_text'])) + return + + def insert_task_panorama(self, panorama_id, task_description): + """ + This function inserts each of data returned by get_task_panoramas into TaskPanorama + """ + sql = """SELECT * FROM TaskPanoramas WHERE TaskGSVPanoramaId=%s AND TaskDescription=%s""" + records = self.query(sql, (panorama_id, task_description)) + if len(records) > 0: + task_panorama_id = records[0][0] + return task_panorama_id + else: + sql = "INSERT TaskPanoramas (TaskGSVPanoramaId, TaskDescription) VALUES (%s, %s)" + self.query(sql, (panorama_id, task_description)) + return None + + def insert_labeling_task_counts(self, task_panorama_id): + sql = "INSERT LabelingTaskCounts (TaskPanoramaId) VALUES (%s)" + self.query(sql, (task_panorama_id,)) + return None + + def query(self, sql, data=None, return_dataframe=False): + """ + This function runs an arbitrary query. + + Returning dataframes + http://pandas.pydata.org/pandas-docs/stable/io.html#sql-queries + """ + try: + if data: + if return_dataframe: + data = tuple(["'%s'" % str(item) for item in data]) + sql = sql % data + return pdsql.frame_query(sql, self._conn) + else: + self._cur.execute(sql, data) + return self._cur.fetchall() + else: + if return_dataframe: + return pdsql.frame_query(sql, self._conn) + else: + self._cur.execute(sql) + return self._cur.fetchall() + except: + raise + + def yield_query(self, sql, data=None): + """ + This function runs an arbitrary query. + """ + if data: + self._cur.execute(sql, data) + else: + self._cur.execute(sql) + + # + # Lazy fetching + # https://github.com/PyMySQL/PyMySQL/blob/master/pymysql/cursors.py + # return self._cur.fetchall_unbuffered() + # + # Or actually, return a generator using fetchmany() + N = 3000 + while True: + records = self._cur.fetchmany(N) + yield records + if len(records) < N: + break + +if __name__ == '__main__': + print "SidewalkDB.py" + + from SidewalkUtilities import modify_city_name + with SidewalkDB() as db: + df = db.fetch_labeling_tasks('VerificationExperiment_VerifyHumanLabels_2', assignment_description='TurkerTask', + return_dataframe=True) + previous_task_ids = df.PreviousLabelingTaskId[~df.PreviousLabelingTaskId.isnull()].unique() + df = db.fetch_labeling_tasks_with_task_ids(map(int, list(previous_task_ids))) + print df diff --git a/crop_from_tohme.py b/crop_from_tohme.py new file mode 100644 index 0000000..fb99d49 --- /dev/null +++ b/crop_from_tohme.py @@ -0,0 +1,218 @@ + +import xml.etree.ElementTree +from utilities import * +import logging +from PIL import Image, ImageDraw +from random import randint + +logging.basicConfig(filename='tohme3.log', level=logging.DEBUG) +logging.info("Logging started!") +class BoundingBox: + def __init__(self, newId, newXmin, newXmax, newYmin, newYmax): + self.panoId = newId + self.xmin = newXmin + self.xmax = newXmax + self.ymin = newYmin + self.ymax = newYmax + + def __str__(self): + return "PanoID: " + self.panoId + " (" + str(self.xmin) + ", " + str(self.xmax) + ", " + str( + self.ymin) + ", " + str(self.ymax) + ")" + + def get_center(self): + return (self.xmin+float(self.xmax))/2.0, (self.ymin+float(self.ymax))/2.0 + + def get_width(self): + return self.xmax - self.xmin + + def get_height(self): + return self.ymax - self.ymin + + +def get_depth_at_location(path_to_depth_txt, xi, yi): + depth_location = path_to_depth_txt + + filename = depth_location + + # print(filename) + + with open(filename, 'rb') as f: + depth = loadtxt(f) + + depth_x = depth[:, 0::3] + depth_y = depth[:, 1::3] + depth_z = depth[:, 2::3] + + + val_x, val_y, val_z = interpolated_3d_point(xi, yi, depth_x, depth_y, depth_z) + # print 'depth_x, depth_y, depth_z', val_x, val_y, val_z + return val_x, val_y, val_z + + +def predict_crop_size_by_position(x, y, im_width, im_height): + print("Predicting crop size by panorama position") + dist_to_center = math.sqrt((x - im_width / 2) ** 2 + (y - im_height / 2) ** 2) + # Calculate distance from point to center of left edge + dist_to_left_edge = math.sqrt((x - 0) ** 2 + (y - im_height / 2) ** 2) + # Calculate distance from point to center of right edge + dist_to_right_edge = math.sqrt((x - im_width) ** 2 + (y - im_height / 2) ** 2) + + min_dist = min([dist_to_center, dist_to_left_edge, dist_to_right_edge]) + + crop_size = (4.0 / 15.0) * min_dist + 200 + + logging.info("Depth data unavailable; using crop size " + str(crop_size)) + + return crop_size + + +def predict_crop_size(x, y, im_width, im_height, path_to_depth_file): + """ + # Calculate distance from point to image center + dist_to_center = math.sqrt((x-im_width/2)**2 + (y-im_height/2)**2) + # Calculate distance from point to center of left edge + dist_to_left_edge = math.sqrt((x-0)**2 + (y-im_height/2)**2) + # Calculate distance from point to center of right edge + dist_to_right_edge = math.sqrt((x - im_width) ** 2 + (y - im_height/2) ** 2) + + min_dist = min([dist_to_center, dist_to_left_edge, dist_to_right_edge]) + + crop_size = (4.0/15.0)*min_dist + 200 + + print("Min dist was "+str(min_dist)) + """ + crop_size = 0 + try: + depth = get_depth_at_location(path_to_depth_file, x, y) + depth_x = depth[0] + depth_y = depth[1] + depth_z = depth[2] + + distance = math.sqrt(depth_x ** 2 + depth_y ** 2 + depth_z ** 2) + print("Distance is " + str(distance)) + if distance == "nan": + print("Distance is not a number.") + # If no depth data is available, use position in panorama as fallback + # Calculate distance from point to image center + dist_to_center = math.sqrt((x - im_width / 2) ** 2 + (y - im_height / 2) ** 2) + # Calculate distance from point to center of left edge + dist_to_left_edge = math.sqrt((x - 0) ** 2 + (y - im_height / 2) ** 2) + # Calculate distance from point to center of right edge + dist_to_right_edge = math.sqrt((x - im_width) ** 2 + (y - im_height / 2) ** 2) + + min_dist = min([dist_to_center, dist_to_left_edge, dist_to_right_edge]) + + crop_size = (4.0 / 15.0) * min_dist + 200 + + logging.info("Depth data unavailable; using crop size " + str(crop_size)) + else: + # crop_size = (30700.0/37.0)-(300.0/37.0)*distance + # crop_size = 2600 - 220*distance + # crop_size = (5875.0/3.0)-(275.0/3.0)*distance + crop_size = 2050 - 110 * distance + crop_size = 8725.6 * (distance ** -1.192) + if crop_size < 50: + crop_size = 50 + elif crop_size > 1500: + crop_size = 1500 + + # logging.info("Distance " + str(distance) + "Crop size " + str(crop_size)) + except IOError: + # If no depth data is available, use position in panorama as fallback + # Calculate distance from point to image center + dist_to_center = math.sqrt((x - im_width / 2) ** 2 + (y - im_height / 2) ** 2) + # Calculate distance from point to center of left edge + dist_to_left_edge = math.sqrt((x - 0) ** 2 + (y - im_height / 2) ** 2) + # Calculate distance from point to center of right edge + dist_to_right_edge = math.sqrt((x - im_width) ** 2 + (y - im_height / 2) ** 2) + + min_dist = min([dist_to_center, dist_to_left_edge, dist_to_right_edge]) + + crop_size = (4.0 / 15.0) * min_dist + 200 + + logging.info("Depth data unavailable; using crop size " + str(crop_size)) + + return crop_size + +def scale_coords(xi, yi, orig_width, orig_height, new_width, new_height): + scaled_x = (new_width / float(orig_width)) * xi + scaled_y = (new_height / float(orig_height)) * yi + return scaled_x, scaled_y + + +def parse_voc_xmls(path_to_voc_xmls): + """ + Takes a path to a folder of PASCAL VOC formatted xml annotation files and returns a list of BoundingBoxes found + in all of the xml files combined. + :param path_to_voc_xmls: Path to a folder containing PASCAL VOC xml annotations. + :return: List of BoundingBox objects + """ + bounding_boxes = [] + for filename in os.listdir(path_to_voc_xmls): + if filename.endswith(".xml"): + xml_path = os.path.join(path_to_voc_xmls, filename) + print(xml_path) + e = xml.etree.ElementTree.parse(xml_path).getroot() + + objs = e.findall('object') + + for ix, obj in enumerate(objs): + bbox = obj.find('bndbox') + # Make pixel indexes 0-based + x1 = float(bbox.find('xmin').text) + y1 = float(bbox.find('ymin').text) + x2 = float(bbox.find('xmax').text) + y2 = float(bbox.find('ymax').text) + + bounding_boxes.append(BoundingBox(filename[:-4], x1, x2, y1, y2)) + continue + else: + continue + return bounding_boxes + +nodepth = 0 +for bbox in parse_voc_xmls("/home/anthony/Downloads/Tohme_dataset/2048_VOCformat/Annotations"): + # Path to VOC data + voc_dir = "/home/anthony/Downloads/Tohme_dataset/2048_VOCformat" + + # Get depth at center of bounding box + bbox_center = bbox.get_center() + center_x = bbox_center[0] + center_y = bbox_center[1] + # Scale coordinates up to full-size panorama image + # Our images our 2048x1024; full size is 13312x6656 + center_scaled = scale_coords(center_x, center_y, 2048, 1024, 13312, 6656) + # Retrieve the depth at the scaled location + scrape_dir = "/mnt/umiacs/Panoramas/tohme" + pano_id = bbox.panoId + depth_file_path = os.path.join(scrape_dir, pano_id[:2], pano_id+".depth.txt") + + # Predict the crop size from the depth + depth = get_depth_at_location(depth_file_path, center_scaled[0], center_scaled[1]) + crop_size = predict_crop_size(center_scaled[0], center_scaled[1], 13312, 6656, depth_file_path) + # Scale back down to 2048x1024 + crop_size = (1.0/6.0)*crop_size + # Load the panorama + pano_path = os.path.join(voc_dir, "JPEGImages", pano_id+".jpg") + im = Image.open(pano_path) + # Make the crops + + try: + crop_width = crop_size + crop_height = crop_size + + top_left_x = center_x - crop_width / 2 + top_left_y = center_y - crop_height / 2 + cropped_square = im.crop((top_left_x, top_left_y, top_left_x + crop_width, top_left_y + crop_height)) + # Save the output + random_filename_padding = randint(1000,9999) + output_dir = "tohme_crops3" + output_filename = os.path.join(output_dir, pano_id+str(random_filename_padding)+".jpg") + cropped_square.save(output_filename) + print("Output file "+pano_id+str(random_filename_padding)+".jpg") + logging.info("Image: "+pano_id+str(random_filename_padding)+" Predicted crop "+str(crop_size)+"x"+str(crop_size)+"; actual "+str(bbox.get_width())+"x"+str(bbox.get_height())+"; depth "+str(depth)) + + print("Predicted crop "+str(crop_size)+"x"+str(crop_size)+"; actual "+str(bbox.get_width())+"x"+str(bbox.get_height())) + except (ValueError, IndexError) as e: + logging.info("Skipped crop for "+pano_id+"; no depth information at point.") + continue diff --git a/decode_depthmap b/decode_depthmap new file mode 100755 index 0000000..aea2dfd Binary files /dev/null and b/decode_depthmap differ diff --git a/reformat_tohme_panos.py b/reformat_tohme_panos.py new file mode 100644 index 0000000..beab489 --- /dev/null +++ b/reformat_tohme_panos.py @@ -0,0 +1,43 @@ +import os +import shutil +def get_immediate_subdirectories(a_dir): + return [name for name in os.listdir(a_dir) + if os.path.isdir(os.path.join(a_dir, name))] + +def reformat_panos(path_to_originals, output_dir): + print("Reformatting files in "+path_to_originals+" to "+output_dir) + subdirs = get_immediate_subdirectories(path_to_originals) + completed_count = 0 + for dirname in subdirs: + print("Moving items for panorama "+dirname) + # Get first two letters of panorama ID + first_two = dirname[:2] + # Check for folder named first_two in distination directory; create if needed + destination_folder = os.path.join(output_dir, first_two) + if not os.path.isdir(destination_folder): + os.makedirs(destination_folder) + + # Look for and copy depth.txt file + depth_file_orig = os.path.join(path_to_originals, dirname, "depth.txt") + if os.path.exists(depth_file_orig): + # Copy to output_dir/XX/.depth.txt + shutil.copy2(depth_file_orig, os.path.join(destination_folder, dirname+".depth.txt")) + + # Look for and copy depth.xml file + xml_file_orig = os.path.join(path_to_originals, dirname, "depth.xml") + if os.path.exists(xml_file_orig): + # Copy to output_dir/XX/.xml + shutil.copy2(xml_file_orig, os.path.join(destination_folder, dirname+".xml")) + + # Look for and copy panorama image + pano_jpeg_orig = os.path.join(path_to_originals, dirname, "images", "pano.jpg") + if os.path.exists(pano_jpeg_orig): + # Copy to output_dir/XX/.jpg + shutil.copy2(pano_jpeg_orig, os.path.join(destination_folder, dirname+".jpg")) + + completed_count += 1 + print("Completed "+str(completed_count)) + + + +reformat_panos("/mnt/umiacs/Dataset/2013_09_19 - Extended Google Street View panorama dataset/GSV", "/mnt/umiacs/Panoramas/tohme") \ No newline at end of file diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..846a430 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,32 @@ +backports.functools-lru-cache==1.2.1 +backports.ssl-match-hostname==3.5.0.1 +cached-property==1.3.0 +chardet==2.3.0 +cycler==0.10.0 +docker-compose==1.9.0 +docker-py==1.10.6 +docker-pycreds==0.2.1 +dockerpty==0.4.1 +docopt==0.6.2 +enum34==1.0.4 +ipaddress==1.0.16 +IPy==0.81 +jsonschema==2.5.1 +matplotlib==1.5.3 +numpy==1.11.3 +pandas==0.19.2 +Pillow==3.4.2 +PyMySQL==0.7.9 +pyparsing==2.1.10 +PySocks==1.5.6 +python-dateutil==2.6.0 +pytz==2016.10 +PyYAML==3.11 +redis==2.10.5 +repoze.lru==0.4 +requests==2.10.0 +scipy==0.18.1 +six==1.10.0 +texttable==0.8.1 +urllib3==1.15.1 +websocket-client==0.37.0 diff --git a/samples/getFullLabelList.sql b/samples/getFullLabelList.sql new file mode 100644 index 0000000..c2d42e9 --- /dev/null +++ b/samples/getFullLabelList.sql @@ -0,0 +1,6 @@ +select gsv_panorama_id, sv_image_x, sv_image_y, label_type_id, photographer_heading, heading, pitch from sidewalk.label_point +inner join sidewalk.label +on sidewalk.label.label_id = sidewalk.label_point.label_id +--where sidewalk.label.label_type_id=1 +--14159 +; \ No newline at end of file diff --git a/samples/sample_crop.jpg b/samples/sample_crop.jpg new file mode 100644 index 0000000..db5201d Binary files /dev/null and b/samples/sample_crop.jpg differ diff --git a/samples/sample_label.txt b/samples/sample_label.txt new file mode 100644 index 0000000..d0ac5e8 --- /dev/null +++ b/samples/sample_label.txt @@ -0,0 +1,5 @@ +sv_image_x 5110 +sv_image_y -688 + +7442 +-708 diff --git a/samples/sample_meta.xml b/samples/sample_meta.xml new file mode 100644 index 0000000..2b0dc8c --- /dev/null +++ b/samples/sample_meta.xml @@ -0,0 +1 @@ +© 2016 Google4500 Sheriff Rd NE4500Washington, District of ColumbiaUnited StatesSheriff Rd NESheriff Rd NESheriff Rd NEeJzt2glc1HX-x3HkUBQUVFAxzPIAxc37gGHm9_vNDzPFCCG8ylLzANTyKMgU-KFZ8s-Ov5rmYxWztlqtvMqT-Q4_j2VT10zXlU1xK5XVLG9pK9Jtf3MAM8Nvjt_Mb-Yz43yej0dhIjO_-b3evx84j0LnBQQEBjQJDUAIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCEfFR4OfQQIUFR4A-hjQR4XZdo_vJmec48UbEreg0TuEiXWX9oKgsW485iRbKJs9Le_AdHy2N-H2OtvdQM20mN_nxHlSH_LCdhLj_19RZSj_QU9dF_gUHrs7ysk9Rf-UWH-ewr292tDsL9fk9K_B_a_52B__4b9_Rv292_Y379hf__mY_379HlA0NodD-2fsL9_k9q_r8P53dm_9T2ygEg9yCNwY__WDWQ73HuyP-QIGve3mh_7yy7SXExMjINf2F6uI_BM_yflOtx7vr-DC2gvkOUIsD8kl_o3WoDCiSPA-z8kF_tbLEDhxAC8r3-vXr3GDBtm5ZPYP8LAfACGDwpFcrLUI_C2_r16ydi_tyA2INbA8WPwHPv9jbkjYmIizJj1z2zf_iGBD_fvqjdunBv610lKSopm2Q7moi1JOHUysNe_obff9U9Le1jQdYDu9ycOGzbT8OdM-_9Bb0DrAXrp6Y0f1Yn-nl3Ao57rP9Kg7_3mpPZ3_O1fp_o_YbX_xIkTA8T66x4e-zvWv6_Aq_s_Id4_rXud-v7Cr-v7d-8u3n-SU_1jY0eM6NkzIGD2swGJUntKJUP_QYMGMSb9GYYZNKjR88zxh_65lg_qSv-eiTo9dVJSYmNTUlzJbJXr_Z8R7_-ITv3D0EZW-zv09j_2l52b-j-S7J7-zdzWv56hf_1_AvQfcc_0z8nJmSLIyvKx_mPGmPzGeJ2hdWz1nzZt2miD1FTDv_23f5bB8HpP6Zj1n_u4jkf6S3v7v1c3E6mpkye3a1DfXxiFof9gnaGmBteZ5VT_aMv-Pb25_3Ni_bOyMjLG1nvcYILOY2amTp3qof4tW8bFmb6wsLAwpVL_q1ATJqG72ZRax2wcpgZbeFpHSn9hAAkCb--vv7gzdR8zMoT-GUb3OSJeILF_C-f6tzQICorTCQ0NM6PUCRXVOGwTS6IDaSI6CvH-_YwSRBm_1C35bfSP17PfXx8y0_jRPG1j06dPt9yA3f4tTPTv39-l_voJBNWvwHIHdrZgexcW-lkn3jmh0UM0PFtgYKDZQXYxCDW8pLqYbdq0CWxax3ShnfRCHOxPGYUIZtjtLx46PsTaNW_66EbZDf1H6bQQIZTv31xHQv9WRsLZC2xpKijIfAWiQ5AwBqcGYiW1SfJAy-wm5QWmL6iNQaB4f5MNhIQMHDjwwXpi_dX1OutR9cT6CxsRy9w5xIrs7M4WrCW3LK_TsaOE_k1amS5AfAPmK7C-BJf3YEegGbHnNilv3j6oob6N_mYb0K9AOJ2i_RsVUJuqq2bMbxmzTgtrn7Db20p5XXsdKf2bWC7A2gZEZmBvCuLjkBLZkrXHNQvfqL1ZfZMBiPW32ICwApH-utuCA3F08R2q6ASL8PXtpfdvPABbG7AyA6fWYLu1neRm4buYHIPIgZvVt3MDENtA_c3A6EEL1vq7KXzzRhra645VYn_xBdjeQN0MrOzAsU1IbC0W3uK5RA_Yor5j_UVHoJ_BQItzbzkGR-4McnU3S9-x4TCl5be6ALsbsHc_cHYZUq53m-0tb_3S-otvwMoObG5CbBvNXckuctmba6vjcH_rC3BsA5LuCNL7W69uK714fYn9bYzAgSHY5Fx0y_Qi7Y319YyvwG7_hgE0XoCUETh-V3Ctuf32Vuo709_2BpydgklvKV9mJ71pfUkv0MYtQGQEjmzA9iSkt5aQ3np9Z_s7NgKLPTh9b7CX3sbzOlnfoQW4PoIGEmtLSd-QX-wVuNBf-gqsk798iMiN3w0LcObbgXz9HXti6_UD7b8BIGEELqxA5vJy1JewANdX4KbyduvLcANoNAN39JcSXrb6EhfgygrcUt6B-vL2N3JiBna7S5uTTPH1JC7AuRlIC-_wrKz-2OfW_nUkfFcQry4teqP4ctTXcWYBEocg2_XeuL6dI3RjfxMu_4Agvb6ch-_CAsSnYJlStuYS63uqvz1yx5e1vo4cC7DB6cou1r-n-rstvp5bF-CO-I7Uv3f6u-_Sr2f7XWGv6S-lfqBMbwC4yvvjG7hrAWD1veQG4BPx9dyzAHnrS3pqH-_v0fh67vg2IGN8afV9ur9nr3wTsi8ArL7v9geLryfzTQAqvq_2B21vJOcC5Kjv5FP7XH9vaG8k2wRgLn3f6-9F7Q1k-j4AFd9bBuCT7Y3kmABYfN_o763p67k6AWfju9reW_pbH4DXp6_j0gTA4ntxf18Jb8LpDUC198b-vpfdlHM_DUhqL2d8b-l_b5G8AYjrHvu7laQbAVB6HIC7ObgCqPTY3yNa2ZsBTHjs73GtxKdgVt0j0bG_N_F0cuzvXbC_fwPt7x3_C6hf847-OAAo2N-_YX__hv39G_b3b7D98S8A0LykPw4ACPb3b9jfv2F__4b9_RtwfxwAMOzv37C_f8P-_g26P74DCAs6P94AYEHnx_6woPNjf1jQ-bE_LOj8OABY0PWxPyzo-tgfFnR97A8Lun4gDgAUdPxA7A8KOn4g9gcFHV8HBwAHur0O9ocD3V4H-8OBbq-HAwADnV4P-4OBTq-H_cFApzfAAUCBLm-A_aFAlzfA_lCgyxvhAIBAhzfC_kCgw9fBAcCA7l4H-8OA7l4PBwACOns97A8COnsDHAAE6OoNsD8E6OomcAAAoKObwgF4HnRzU9jf86Cbm8EBeBx0cnM4AE-DLm4O-3sadHELOAAPgw5uCQfgWdC9G8EBeBR07sZwAJ4EXVsEDsCDoGOLwQF4DnRrcbgAT4EubQUOwEOgQ1vTFBfgEdCdrcMBeAJ0ZVtwAe4H3dg2XIC7QRe2BxfgXtB9HYATcCPouA7Bvwy4DXRah-EE3AI6qxR4G5AfdFPJcAOygs7pHLwTyAW6pEtwBS6DTigPHIKzoMvJD783SAFdy7OamtO9fugAwKCLgIMOAAz69IODDgAM-vSDgw4ADPr0g4MOAAz69IODDgAM-vSDgw4ADPr0g4MOAAv67HsF6AhgoE-8FwiwMLL1t4rUeQHK6rsc37U4iflrn0WamSUVZTN_5_g10xOZoQOXkl8L_5TM_ofj39-tZPpFFPAtrvSnTrVbrrmzulj9SeZWzaox2cqlv3B8YdkA5od9NR-8UtJCu_gMx3d7TMEo9m0s37F2Lr8y-Lhy8czmLJf4X014-JHkQ3c4PvMHFdN2xhbyybj5ytLLHD9kyUDm3IoxJHgVrbz4E8efy0tkXrzVnyzbGq0adZzj912azow91lZboRWOdUmVMnb2DXVX9avkky7nVD_lFvHlB4uZbqNe0UQujFOlHeb47I0PM7FN_qg8f5bjy4urlM_-GM2-e4Oj_5u1XVmyeULZrJDD6tykPH7I9RXUzlu3NWka4fNHczSbr9RoVT2L-FaLWzPN4o4oDk0t1cZQHF8wahzzwV_z-Y_iS6iyaec0-xf2VJ9PKqVuvfY8_dB7J0j3RWo270oiaXcoXfmmcL60p9MY5l8cHTEiRJVe-nbZnC7R6hvao8rakxwdu7yKFMYksBlVBXTelc7lxw4UJ4-uuKhuExxAbtQElw8-WcjnDn-DaRvEkiObP03majm-9fgMZvChOFXHzzj644Iqcu3AErak-qJmWV-NKuxKEb-2tJbOUOaWP9X7JT7x0Urlv09uZyefiCYbEkcrd_3M8Rc7DGaOHkrh2XPd6Ccrl5EFXVLYkuNzyZHD87WnOY7fQr3FhDVrQkb-p5LkXBLOf_9OzM0ijl_4XLFq84CysidH7GOpp-aQFxOU1HfrFvI38waoh62fSy3dMp8OvnCKVP_4MrunNIhknn1Hmz-C4z-__TXd8a0bmpbR8eSUcPzvJC1ipnBz-ILsrVTawjZkaO-r6h7TYlQLyzg6cnEV-T5pJXt5-x_I5eVNqZqDhTw_7xYT8EYxdWzni3R19CkyQ9uGZaYt4F9oOppKzNmiWV32s_rZftHljz9UyP_6-mnlqp_L2XObblAnvp1M5504QOYcf4Y9dSmaXj83nR7LfE6qLu5Xt79cxO_fRlRTvo7U1J5fqZ4YcZ8ibtpg5dXfOP4Oncoc2PQeWRC5X5X1VRE_ZU45vfuB2VTgqfn054GVZM97W1hFxjqy__8CqHVXC_kFO1-go15rSq3bVsjvnHZaOfbp79ia6J1kV-ebyYOvcfzGUgXD_DSdr_jiLPXKyiRS26kvc-kuR2-PCCJbzgzd027HSnbNslh-SkEqXfPcp6TmWhI7fOcCOiwyvfzYI5uSp4aeUP-ZKaR_WBFdflMzO_nAOqIu--CEZtxHfbQTD3D82PhKusvocH5p6Bi634JdZGR8Syazskrz2NZibUghxyfP3cXM-nVS-Wd_fom_tLVSeWnkXXZkq63U8wVz-OulXyof1lxl992-n-xWfkP-Vs3xU6qXMJOucvQ3vQ8qT9-dWcbteFr9xAWOpuJvkVeHhSnCdm1jI3-7QCY_sE2luV3EL7wvWJ30wQJeNX4U9fUzH2uWTWyasmXv29qXkzm64uszRHvuJPsh9a6qexeO77PvjDJ3R7uUhNw3ybhbfZLPCtd375m5zPnvOHr4pDuEhEcpviyfxP77ixXkTt4_kyNrhL1SCmbctdfLPwzJ50dvO6ncm90y5fmFYXTo6jF0depeonm7L6MZPZs6eWI-HRVQSYYe3cuOuLqI_DF8nXZwknD9ZbdVb38hl58Q_HfqtUf7kieGZLHdXkkli66PUI38hOMnZF6hT312l6xfc0d748dCPmxIvvrv3y9Tfv8Nx896uUqZsXIW2yfiIrlz-xftT78X8pm1mer1P-ruh3_RfCkc__TYtUzI8mLSQbtKtWg4x1PVD6rfCJvHl2_-iAr6pQU50qKE3dPxFqntkJIcfYHjh8_IYV7dsoDftSKdit-zSROZ_gXbMWADmbd4s-pUM47Pn_ELM3Avx3dQUdrz-RWK9IoL7G_TnyonZ17iD5-uVD4ad5ZOeuNN8v87vko-KJwfrlUaEzPoIxK_cY522CKOH70kQb2p_XFCffqxcoNwzytoWsAE7C7k738plFqZOF6z9vJhVrM8T7V6HkfvHlVFOtYqUlb_LY-UvXWevPktx8d8-hZ9rc1iKuHlF_miUyeVvdp2SQkb1p58eJYifxHu75W_72YuXE8gnZLOK-cIr6dD8Xf0sn1rFbcqvlKd2FxEf788VF3bJpBczWym6iHc__vl3KCblISTqzHjVKPe5_gWxyvo3KjXyc7gDdrmvTk-vGCNesL1EurYqnx6X_w_CPdNNXugRxzJKVuhGpQlPN-mrVTXfxTwJe89QM18epGmacJ6drviK3L7tXDq_aJCPqNnhvrkQxHUyamF9KS9p8k7p6vZrqUrSOyGg9ojh4v4roOWqu_7E0drasZr84-uUUy8eZRNe_x5em3Ld8u57GDlF-Eb1RmdjiqZmfnayy9w9NqsUvp_qLg2qQeJzt2Q1wm_Vhx3EnJokTJ3XiQAiUkSUQICVq9f5CiMnz6P1dj97lA3TW22O9v7-22fXtrnTsOm69a7u2FEZfVlJeRkaPKy9JJNN1MFp6XV-gsOtuR6_cxih0vW3ddu0kK3YUS7IlS_Ij4d_nLja2rEd_-_uTcg4TzrGxrWNbJsYAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAARtSWLUyfABi0dcslTJ8FNtx4Y_-eZnBFo_4eEgZlvFX_LVu2VXV6iStaGeSZoW_GV-m_9ga2b2_ZHv1Hxlr9225gR9V29B9x4530XzmBHUvQf8SNd9p_YmJicQM7LoP-I66r_tU_Ozruz_R3Bp1A_01tJ_pvat3034b-7znov7mh_-aG_psb-m9u6L-5jVj_XbsmJ3fv3j2IS29O6L-5ddt_x1D0f68sYE_N-5g8wQD7N9y3b8d9T_ZfxNAJ0J9JDf3f19UMtvbrBBvTf7Jfx33P9-9wAVur-nIC9GdST_2bFjC1jhPg9Z9JPfZfsYCpdQxg-Prv3bt39759bW5E_62XGa9_cufi26mp6f3dnqB9_22M9L_yyj72v6qm_rb6fgit3f_ARePjtdu2rlTvv3Pn1VOLRrb_wUXXHBxA_yXVO03unpy8tr3az7mLH10frNX_wIFN2__91_3R9ddff-iPDx06tG_fvoMH61_X2L-e9fCRGxbdeGPzVdfRf2MXsHWN_kf71_-mmxdN31J37KIh6_-Bev-Dzf3rE2juf6QK_Tvrv3__9HD3v7V1_0Nj9X67J5f7Vy33v-qq1v2Pr6v_7qtYrNq9P_ihMXZ3NbvXh_7T09P7G_pzuFwOp-lxbu69_8TQ9-etvGhP_fnsKoFAIBSKRGKxSNRD5fZ67y9p3f-2E1XLl-Fya5_lcNr27-jXP_TvuwH1v-3EYPo3_-_fAf3-t_xh3_uvkv9if1ZD_9tHuv_Jkydnqu64Y8T6X_bb3akGhw-v1r-KqDt-vP528_YnG5yQ3lSz8u__Wzau_7Vd_VD2Xmx90TXXNHwoW-pfLV7PuvjL35FGcrlCoVBWqdbV_8DK_sKR60-S-_fXnu11Uxdvq7llWS2--phmg_prdVpt4zem1Wr37l36rwaTS3avajnq0gqWHGlHX9NN_zE221A17P2Npqpqf4qqhZ5e0XxqytzolhU27Pmv60itweRqVh9F-5Usad1_6csbH-moZZFhkbVqIPlX6X_1orX7L0Zefv4vu3qlVvlruuk_se7-E810l33U7im56hy6HUgXFzp6iaWBoZWlmFZry5uXNmTosH_9kcbHd1aNrdm_KXTdznbGm2w51rp_0z8ALeqif4voLe1qsNrL8_om0Y2mvwoOHDi6Qr1O7f3RFZto3Ill5d0urchqs11aRav-9a9v7lRbQHP_nTvtnWVuVg_YMv2KFTTModf-DVPatZoOhtDlQrq6xIEmTUUbyra9rf2d6ntp0__AgUuB21jj5tZfvWxlxbVnsGR9v_81k8lkq_bv0xTWqTl_C5cF7egey_dcet-if_2Wpid6K10EXz1-FyvoS3_Zsk43sEFjWH6YrnL2rmk-za_0_bFK-86W0Gt_WQtdj6Bfe2h_2T0NNngLi_kbbXj71XfQaf8OwvdpA_23p6X-Nl6-YPPjvK-dWsaO9tD67utaQLuXgs77r5F-CDfQun_LHbQex-oXWFPbAfRsfRNov4Ka1v0v3dJp_6EZQY_1erd9--AW0MMG1hpCm0l01X8YNsB0fodj-0VDuoH2I2j5UtD1AJgdAbPp67Zfrt_5L162xxVcmkK73wF6GQBTK2C8fcsF9D6E1hfsdQdrlO95ABu-AtkQtF91Ap2MYa179mMJnZXvzwA2ZAXLDzQM7TsawAZab_h-DmBAO2jxEEPQfugWULeu-H0eQF-GsMa1h6D98C1g3fEHNICuxtDtBZlPP2QL6Kn-4AcwAIynH6IF9Fp_JAfQ1wWsp_yQLKAP8TfzAnopPwQL6Ff90RxADYPdmV5AH-OP8AC6nEB_wzO4gD7XH-UByKoBGGjO4AL6H3-kB7Ahgdc20vFHeQBMh182yvFHeABMZ28wwvFHdwBMR7_cyMYf2QEwXXylUY0_qgtguncLo9l-RAfAdOzWRrH9aA6A6dJtjV77kRwA05lXM2LpR3IBTDdew0iVv4jppF1hOvDaRif8EqabdoPpup0ZiewNmK7aOabLdm67g-mqXWA6a8eYrtoFppt2hemuHbriCqardo7ppF1iOm0HFs_JdNZOMVtzHZiuu6algzJdthMMdlw_pgOvquGcTMddG1MFe8R041VcflCm-66OkXb9wXTmNprOyXTi1Wx8tX5iOnUrrc7JdOV2NjhX_zEdu0m7gzJdupUN7DQ4TAe_zCrnZDp2s40qNGhMR1-2xjmZ7n25jSizUZgOv6iDczLd_JJBF9loTMeXdXhOprvXDbIEU0ahfg3T7d-b9WtGoX4N6g_MCNSvQf2BGf74i1B_cIa__iLUH5yhj78I8Qdo2OPXIf4ADXv8OsQfoCFvfxHiD9Jwt78I8QdrmNsvQfsBG970y9B-8IYyfCOk3zjD0rwFlIcaZAcAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAABgjD-YMjm8s1KzPiEzzsqjMbNmPk1LHToZoSxyEjJabyPMBsLjE5foDCfM11kSXoNQp1KaAk5NzhnmER4HrfDYgnlpRpLIBiQFX9SUsRFsbcqoTZZknpzRYPSmwzy-lVTT8YxLW3C7pFmqFBFF-Fa3lHZQlDY5G4klND5CZi1kVCYuRyMIOChZcl7h4XKSSVOKk7aqFCpDwmgiBPMagyUU1GYTKik3Ezc4UlqPPy_JmL3pOYtfEjYS0ZIhQ6tVMUc-n60-0ryYr5XPFzwxhzUmJ7gpgqa5bNuctiDJl-xxnTed1kZmWVTQHcuGwkmvJeyIFgwOdT4iT0do1ryBiHIiAl7OIUyY5mXGQCYjomKBWGHOnmAXQ64Q2-7hiA0ykTjvypJKMismrWaR0lJMJkvsqIGbDwdNBTWdCcmLYr7G71SLCXshQOgsXkswJ87YQlKxzeyMJk1miTifC2l0bkOSLfOyfLNzWb5XFHCJKdWsf96nj9B8eTzNi3mKSqNAqaRDZo004JSWbFEXzyqT5SICrS6cL9rNvpiD8rHZhlK8KNOKnRYFh0jKROkSW-MWpOU6IhPmqwM2ik3Hubk5PUsv88s1Ok_MmrRlJF5-Xl6wOlSuTKpIGLgcdcwkzTjZ7ryYckgppS8VCus4ERVPqRPqo2ZqvpRVpwR2p0roMKlNNl2CG-K5qLgi4FAqYjFnJmHkOfUJ0k6TccKT9pujfKUvkAiZ3AWTsZAPeZ2E11aivXyWIJnheaMpNcebCuac-Xm5nFNyl7Jxm0TkigSLfI9HGKQlHhtdksdYNqM_7FTH3TKLJ6Smw_aQ1ln9hiOaQDSfUxvNyqhemBfPFbMJnjNWkuporSGcYLmEegkt87EzplS8YGS5kiZONu41RlJZelbryfDDeZKnZQmkipDFxSUtAiKQklACXVJG-_JZ0ZyapiOCpINOK6N2GztHOTKkn2fkaQX2MBkmpDxpUaHjSRXOfIBHp4M6gSdVXbiXT_ndCrFSmhVrtVGRJVAwztMcrTnqjil48ahWkwwJU4o5hTBBsTlJldkQKSiinoCYynHFMpfTRgWS5kLQXwxlsjlutCDUEUpbliB80oyFNz-rjkcEWSXBToadJZWRmFPNc9QKrkGhiMb9VrlIm5pLEIQtZJGzZHKpPTHrCKQsVCzBFhN5c7DotIcj-ozKK2bzVKZQwK6O01Z2Jmbk0hmyGPe7bI6SzhSO84k8K-NymH18vt4R55ssoTDbntYSxSDPblHpdQp9xs2VFOb4-XlRkaAzPEIqNtJR9mxa6ZXI0jKuwkmZiny12S5V5ozJYvWZ5AmwHAYLVTILldykoZSQByVSwhQ0irU5Vt4miNl9Qf6cMGEzp5J2Dy8jShpp_7zWXMpS4YQy5Z-neP60lC8TU_wQn2WxiaMRfkIj8MzGtBmPmkO40hlz3mXJFCOhEDGn15jNbjuHniuWYo6oPCUPhOTJojBjj6Z9ViJVfQGqflZABWiT0xK3GknKEKOshMgmd2j90ahcY-QWg7aIgpUiuAq1UJrxO0wyW1I76xGJoj4Pn3Ykqt9IPhMyJg3BiF_tm7PqzHO03MSTu10-nySRVBaD2disge_iEvko36VXFa3ReU1SHnDpFW7Srs4T8tnCXClgMZhUs5m8VWhMh0therYY8aeIgDSfUJh0yoCKN58l-YSbDJEFOqx38vNuWp_0ZxS03kmbXNKwxG40cqovXMJ4XDarVLs8FKvkNkiUtrzFLdIaONyE0M7XqClaaXSoVUQgI5RziIgkwGbRwTzlSlCRkNROyxTzah2P5WWzlaQ9r7HIqaywVJjzylI5pypctJAOE8XzJaIONy-ncoR1Gm1RapeHOQHam7JRKqtBFmJFncGQMC0vJg1ZIm00JMV0wCBw57R6XtzFIvJZM6WMZhz6iJwSEwJ_KKyI-CxJriovCLr1cbmAz5IpJSa3lC93JgRmSuwTJ50adiRJy4R8DcttSFh5AYLMSl1OboKSeVglq89G8bWkiE17CiGFRZjwx4yeYpwVUHsTDqWaUMq1Ng2tyMnNZEkc1BOkgcsT5K0padJOCWfnvO5CQSGSyezmnIcqSbmeQlwXcBmz1qw7Z5f7g4SQiGZUswTpcjpKWUVwXhVP8f1EggjnhHM8uYWSaDRFJZkUaAuzQSshMSQdpchsSG1MRf0xc8yWVIfEdNjvyfmjHJHYqggqCA6l8xs0Vvq6Hz1UfuGznyy_vv_Iuf8Wv1L-Vu6Bc4--eWVlbGzso9U_Y6_c_2b5PwV3ll_MvXn7Pzx3vHyLW1x-9Y0_K4u-_P7Kf8wIy85fTJx_8PMnKt994dsnLPGPEk_f-PDJl566kjT925YZwcmT5J3Wu4gnX9l7x_Ejd588duLj5WnLxIWb936pYnnnM5W_MXzywldP_EWl8ubz5N2-2y-UJY9XPnzw9AVasVCZij1YeewTWy78-KbTM6J7A-TbZ18iPv7arjv-7r758-onvBUr-ZHK2X89R55_21w5mfkO-Wqcc_Ifn_oJ8XP27yqP_Pk_kVLzG5WXp18lPyX9fuXlQ2-Sr_7uu5XPvfwHclz0jcpv_v1F8rmjT1eImR-S6Ud-X_acf5w8pPr1zIfd95Ayw0Pl6949Q_762X3lzxw5Q97-2_2k9K2_mnl82w_K1z2wbUH78FRl29e3Lnj-9JnysYlUxfHBXZXp7wUq255-tfKRN35K3vfGDyrPfvPH5B_2Kc-_9p090luf_yVxzPAEeXDHc6TtwhlSL_wR-dz0GZL77gly_-uPkzckv3GB-9LPKr-_8ycXJs-9VfnKC4-ee336dfLEV_zk2f97jHS-Nil94whBLqgmpX_ywH5yWn0PSf_9o2RoelL68Jf_h_jWF24jD_9s18LYl753SuTZuXA4ySqrv7Z34fZTu6VnPh0ln3z7gbLpQePCL4LfLN_1om_h8Pd_XlbemVi4hz5WufZvnySveHby_FdvfJe85qCqfM_91y489JvPliWxHQt_7Tc-85e8M6ee-VzswumrHiXdj_1q5gnzWfJr73xx5vFvP0t6P_FfM3dn7yfPPuKpqHT_W_npqQ9VDI_cVbmHOF3JH3_-wlu_1VfyL0yVr7z5fkL9-bPkU7E7yNOumZnA6T3krcldC0e_-BZxr2rXwj-L7iP2_MvOhWMv3Vr-wMc-VnnnV8fIw898euZpeZh0_dB57rYj95I3UIbz8a8rSY0hvnD9F-479cu9N174f9xC1qg \ No newline at end of file diff --git a/samples/sample_pano.jpg b/samples/sample_pano.jpg new file mode 100644 index 0000000..22a7d1e Binary files /dev/null and b/samples/sample_pano.jpg differ diff --git a/utilities.py b/utilities.py new file mode 100644 index 0000000..0e086d2 --- /dev/null +++ b/utilities.py @@ -0,0 +1,494 @@ +import numpy as np +import os +import scipy.stats as spstats + +from math import radians, cos, sin, asin, sqrt +from pylab import * + +try: + from xml.etree import cElementTree as ET +except ImportError, e: + from xml.etree import ElementTree as ET + +EARTH_RADIUS_M = 6371000 +EARTH_RADIUS_KM = 6371 + + +def basic_stats(value_array, verbose=False): + """ This method takes an array like object and returns median, mean, stdev, min and max values """ + if len(value_array) == 0: + return (0, 0, 0, 0, 0) + median = np.median(value_array) + mean = np.mean(value_array) + stdev = np.std(value_array) + min = np.min(value_array) + max = np.max(value_array) + total = len(value_array) + + if verbose: + print "Median:", "%.2f" % median + print "Mean:", "%.2f" % mean + print "Std:", "%.2f" % stdev + print "Min:", "%.2f" % min + print "Max:", "%.2f" % max + print "Total:", total + + return (median, mean, stdev, min, max, total) + + +def bilinear_interpolation(x, y, points): + '''Interpolate (x,y) from values associated with four points. + + The four points are a list of four triplets: (x, y, value). + The four points can be in any order. They should form a rectangle. + + >>> bilinear_interpolation(12, 5.5, + ... [(10, 4, 100), + ... (20, 4, 200), + ... (10, 6, 150), + ... (20, 6, 300)]) + 165.0 + + Code written by Raymond Hettinger. Check: + http://stackoverflow.com/questions/8661537/how-to-perform-bilinear-interpolation-in-python + + Modified by Kotaro. + In case four points have same x values or y values, perform linear interpolation + ''' + # See formula at: http://en.wikipedia.org/wiki/Bilinear_interpolation + + points = sorted(points) # order points by x, then by y + (x1, y1, q11), (_x1, y2, q12), (x2, _y1, q21), (_x2, _y2, q22) = points + + + if (x1 == _x1) and (x1 == x2) and (x1 == _x2): + if x != x1: + raise ValueError('(x, y) not on the x-axis') + if y == y1: + return q11 + return (q11 * (_y2 - y) + q22 * (y - y1)) / ((_y2 - y1) + 0.0) + if (y1 == _y1) and (y1 == y2) and (y1 == _y2): + if y != y1 : + raise ValueError('(x, y) not on the y-axis') + if x == x1: + return q11 + return (q11 * (_x2 - x) + q22 * (x - x1)) / ((_x2 - x1) + 0.0) + + + if x1 != _x1 or x2 != _x2 or y1 != _y1 or y2 != _y2: + raise ValueError('points do not form a rectangle') + if not x1 <= x <= x2 or not y1 <= y <= y2: + print "x, y, x1, x2, y1, y2", x, y, x1, x2, y1, y2 + raise ValueError('(x, y) not within the rectangle') + + return (q11 * (x2 - x) * (y2 - y) + + q21 * (x - x1) * (y2 - y) + + q12 * (x2 - x) * (y - y1) + + q22 * (x - x1) * (y - y1) + ) / ((x2 - x1) * (y2 - y1) + 0.0) + + +def chunks(l, n, lazy=True): + """ Yield successive n-sized chunks from l. + http://stackoverflow.com/questions/312443/how-do-you-split-a-list-into-evenly-sized-chunks-in-python + """ + for i in xrange(0, len(l), n): + yield l[i:i+n] + + +def distance_to_latlng(path, distance, heading): + """ + This function takes a path, GSV image point, and heading + http://www.movable-type.co.uk/scripts/latlong.html + """ + with open(path + 'meta.xml', 'rb') as xml: + tree = ET.parse(xml) + root = tree.getroot() + yaw_deg = float(root.find('projection_properties').get('pano_yaw_deg')) + yaw_deg = (yaw_deg + 180) % 360 + lat = float(root.find('data_properties').get('lat')) + lng = float(root.find('data_properties').get('lng')) + + + #R = 6371000 # Earch radius in meters + R = 6353000 #Wikipedia + #R = 6384000 + d = distance + # bearing = (heading + yaw_deg) % 360 + #bearing = math.radians(heading + yaw_deg) + bearing = math.radians(heading) + + lat_radian = math.radians(lat) + # lng_radian = math.radians(lng) + + plat = math.asin(math.sin(lat_radian) * math.cos(d/R) + math.cos(lat_radian) * math.sin(d/R) * math.cos(bearing)) + plng = math.atan2(math.sin(bearing) * math.sin(d / R) * math.cos(lat_radian), math.cos(d / R) - math.sin(lat_radian) * math.sin(plat)); + + plat = math.degrees(plat) + plng = lng + math.degrees(plng) + + return (plat, plng) + + +def ensure_dir(path, verbose=False): + """ + This function checkes if the given path exists. if not, it creates a new path + http://stackoverflow.com/questions/273192/python-best-way-to-create-directory-if-it-doesnt-exist-for-file-write + """ + if path[-1] != '/': + path = path + '/' + + d = os.path.dirname(path) + if not os.path.exists(d): + os.makedirs(d) + if verbose: + print 'The directory "' + path + '" does not exist. Created a new directory.' + else: + if verbose: + print 'Directory exists.' + return + + +def find_outliers(data, mode='quartile', m=2): + # Similar to reject_outliers function. Howerver, this function returns the indices (a boolean array) of the + # outliers instead of the actual values + if type(data) != np.array: + data = np.array(data) + + if mode == 'stdev': + return abs(data - np.mean(data)) > m * np.std(data) + else: + percentile_lower = spstats.scoreatpercentile(data, per=25) + percentile_upper = spstats.scoreatpercentile(data, per=75) + + upper = data > percentile_upper + lower = data < percentile_lower + return upper + lower + + +def frequency(arr): + """ + This function returns a frequency of items in an array + """ + frequency_holder = {} + for item in arr: + if not item in frequency_holder: + frequency_holder[item] = 0 + frequency_holder[item] += 1 + return frequency_holder + + +def haversine(lon1, lat1, lon2, lat2, unit="km"): + """ + Haversine formula + Calculate the great circle distance between two points + on the earth (specified in decimal degrees) + http://stackoverflow.com/questions/4913349/haversine-formula-in-python-bearing-and-distance-between-two-gps-points + + It returns the distance in km by default + """ + # convert decimal degrees to radians + lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2]) + # haversine formula + dlon = lon2 - lon1 + dlat = lat2 - lat1 + a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2 + c = 2 * asin(sqrt(a)) + #km = 6367 * c + km = EARTH_RADIUS_KM * c + + if unit == "m": + return km * 1000 + else: + return km + +def histogram(hist_data, header=None, x_axis_title=None, step_size=1): + """ + Takes a list of interval data and plot. + + !Deprecated! Use Chart.py. + """ + + med_val, mean_val, std_val, min_val, max_val, count_val = basic_stats(hist_data) + annotations = [ + 'Count: %d' % count_val, + 'Max: %.2f' % max_val, + 'Min: %.2f' % min_val, + 'Stdev: %.2f' % std_val, + 'Mean: %.2f' % mean_val, + 'Median: %.2f' % med_val + ] + + frequencies, bin_edges = np.histogram(hist_data, bins=arange(0, max(hist_data)+ step_size, step_size)) + + + fig = figure() + ax = fig.add_subplot(111, autoscale_on=False) + + # + # Define the domain and range + if max(hist_data) < 1: + ax.axis([0, max(hist_data), 0, max(frequencies) + 1]) + else: + ax.axis([1, max(hist_data), 0, max(frequencies) + 1]) + + # + # Put the stats about the data + for i, annotation in enumerate(annotations): + ax.annotate(annotation, xy=(-10, max(frequencies) - 15 * i - 10), + xycoords='axes points', + verticalalignment='top', + horizontalalignment='right', + fontsize=14) + + if header: title(header) + if x_axis_title: xlabel(x_axis_title) + + # + # Use rhist and rstyle to prettify the graph + rhist(ax, hist_data, label=None) + ax.legend() + rstyle(ax) + show() + return + + +def interpolated_3d_point(xi, yi, x_3d, y_3d, z_3d, scale=26): + """ + This function takes a GSV image point (xi, yi) and 3d point cloud data (x_3d, y_3d, z_3d) and + returns its estimated 3d point. + """ + xi = float(xi) / scale + yi = float(yi) / scale + xi1 = int(math.floor(xi)) + xi2 = int(math.ceil(xi)) + yi1 = int(math.floor(yi)) + yi2 = int(math.ceil(yi)) + + if xi1 == xi2 and yi1 == yi2: + val_x = x_3d[yi1, xi1] + val_y = y_3d[yi1, xi1] + val_z = z_3d[yi1, xi1] + else: + points_x = ((xi1, yi1, x_3d[yi1, xi1]), (xi1, yi2, x_3d[yi2, xi1]), (xi2, yi1, x_3d[yi1, xi2]), (xi2, yi2, x_3d[yi2, xi2])) + points_y = ((xi1, yi1, y_3d[yi1, xi1]), (xi1, yi2, y_3d[yi2, xi1]), (xi2, yi1, y_3d[yi1, xi2]), (xi2, yi2, y_3d[yi2, xi2])) + points_z = ((xi1, yi1, z_3d[yi1, xi1]), (xi1, yi2, z_3d[yi2, xi1]), (xi2, yi1, z_3d[yi1, xi2]), (xi2, yi2, z_3d[yi2, xi2])) + val_x = bilinear_interpolation(xi, yi, points_x) + val_y = bilinear_interpolation(xi, yi, points_y) + val_z = bilinear_interpolation(xi, yi, points_z) + + return (val_x, val_y, val_z) + + +def point_inside_polygon(x, y, poly): + """ + This function checks whether a given point (x, y) is in a polygon poly. + + From http://www.ariel.com.au/a/python-point-int-poly.html + """ + 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: + xinters = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x + if p1x == p2x or x <= xinters: + inside = not inside + p1x,p1y = p2x,p2y + + return inside + + +def points_to_latlng(path, points): + ''' + This function wraps point_to_latlng to get latlng coordinates of a list of points + ''' + latlngs = [point_to_latlng(path, point) for point in points] + return latlngs + + +def point_to_latlng(path, point): + """ + This function converts a 3D (x, y) coordinate on depth data provided on a GSV image into latlng coordinate. + :param path: + e.g., '../data/GSV/rP_WcfFFp3V23ESWa59p4Q/' + :param points: + e.g., (18.3720218935 -1.45833482249) + """ + # xml = open(path + 'meta.xml', 'rb') + with open(path, 'rb') as xml: + tree = ET.parse(xml) + root = tree.getroot() + yaw_deg = float(root.find('projection_properties').get('pano_yaw_deg')) + yaw_deg = (yaw_deg + 180) % 360 + lat = float(root.find('data_properties').get('lat')) + lng = float(root.find('data_properties').get('lng')) + yaw_radian = radians(yaw_deg) + rotation_matrix = array([[cos(yaw_radian), -sin(yaw_radian)], [sin(yaw_radian), cos(yaw_radian)]]) + + rotated_x, rotated_y = rotation_matrix.dot(array(point)) + + # + # http://www.movable-type.co.uk/scripts/latlong.html + #plat = lat + rotated_y * (0.00001 / 1.1132) + #plng = lng + rotated_x * (0.00001 / 1.1132) + R = EARTH_RADIUS_M # m + d = math.sqrt(rotated_x * rotated_x + rotated_y * rotated_y) + bearing = math.atan2(rotated_x, rotated_y) + # bearing = -bearing + + lat_radian = math.radians(lat) + lng_radian = math.radians(lng) + + plat = math.asin(math.sin(lat_radian) * math.cos(d/R) + math.cos(lat_radian) * math.sin(d/R) * math.cos(bearing)) + plng = math.atan2(math.sin(bearing) * math.sin(d / R) * math.cos(lat_radian), math.cos(d / R) - math.sin(lat_radian) * math.sin(plat)); + + plat = math.degrees(plat) + plng = lng + math.degrees(plng) + #plat = + + return (plat, plng) + + +def reject_outliers(data, mode='quartile', m=2): + # http://stackoverflow.com/questions/11686720/is-there-a-numpy-builtin-to-reject-outliers-from-a-list + # http://stackoverflow.com/questions/2374640/how-do-i-calculate-percentiles-with-python-numpy + # http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.scoreatpercentile.html#scipy.stats.scoreatpercentile + if type(data) != np.array: + data = np.array(data) + + if mode == 'stdev': + return data[abs(data - np.mean(data)) < m * np.std(data)] + else: + percentile_lower = spstats.scoreatpercentile(data, per=25) + percentile_upper = spstats.scoreatpercentile(data, per=75) + + data = data[data < percentile_upper] + data = data[data > percentile_lower] + return data + + +def split_list(alist, wanted_parts=1): + """ + http://stackoverflow.com/questions/752308/split-array-into-smaller-arrays + """ + length = len(alist) + return [ alist[i*length // wanted_parts: (i+1)*length // wanted_parts] + for i in range(wanted_parts) ] + + +""" + Styling matplotlib by Bicubic + http://messymind.net/2012/07/making-matplotlib-look-like-ggplot/ +""" +def rstyle(ax): + """Styles an axes to appear like ggplot2 + Must be called after all plot and axis manipulation operations have been carried out (needs to know final tick spacing) + """ + #set the style of the major and minor grid lines, filled blocks + ax.grid(True, 'major', color='w', linestyle='-', linewidth=1.4) + ax.grid(True, 'minor', color='0.92', linestyle='-', linewidth=0.7) + ax.patch.set_facecolor('0.85') + ax.set_axisbelow(True) + + #set minor tick spacing to 1/2 of the major ticks + ax.xaxis.set_minor_locator(MultipleLocator( (plt.xticks()[0][1]-plt.xticks()[0][0]) / 2.0 )) + ax.yaxis.set_minor_locator(MultipleLocator( (plt.yticks()[0][1]-plt.yticks()[0][0]) / 2.0 )) + + #remove axis border + for child in ax.get_children(): + if isinstance(child, matplotlib.spines.Spine): + child.set_alpha(0) + + #restyle the tick lines + for line in ax.get_xticklines() + ax.get_yticklines(): + line.set_markersize(5) + line.set_color("gray") + line.set_markeredgewidth(1.4) + + #remove the minor tick lines + for line in ax.xaxis.get_ticklines(minor=True) + ax.yaxis.get_ticklines(minor=True): + line.set_markersize(0) + + #only show bottom left ticks, pointing out of axis + rcParams['xtick.direction'] = 'out' + rcParams['ytick.direction'] = 'out' + ax.xaxis.set_ticks_position('bottom') + ax.yaxis.set_ticks_position('left') + + + if ax.legend_ <> None: + lg = ax.legend_ + lg.get_frame().set_linewidth(0) + lg.get_frame().set_alpha(0.5) + + +def rhist(ax, data, **keywords): + """Creates a histogram with default style parameters to look like ggplot2 + Is equivalent to calling ax.hist and accepts the same keyword parameters. + If style parameters are explicitly defined, they will not be overwritten + """ + + defaults = { + 'facecolor' : '0.3', + 'edgecolor' : '0.28', + 'linewidth' : '1', + 'bins' : 100 + } + + for k, v in defaults.items(): + if k not in keywords: keywords[k] = v + + return ax.hist(data, **keywords) + + +def rbox(ax, data, **keywords): + """Creates a ggplot2 style boxplot, is eqivalent to calling ax.boxplot with the following additions: + + Keyword arguments: + colors -- array-like collection of colours for box fills + names -- array-like collection of box names which are passed on as tick labels + + """ + + hasColors = 'colors' in keywords + if hasColors: + colors = keywords['colors'] + keywords.pop('colors') + + if 'names' in keywords: + ax.tickNames = plt.setp(ax, xticklabels=keywords['names'] ) + keywords.pop('names') + + bp = ax.boxplot(data, **keywords) + pylab.setp(bp['boxes'], color='black') + pylab.setp(bp['whiskers'], color='black', linestyle = 'solid') + pylab.setp(bp['fliers'], color='black', alpha = 0.9, marker= 'o', markersize = 3) + pylab.setp(bp['medians'], color='black') + + numBoxes = len(data) + for i in range(numBoxes): + box = bp['boxes'][i] + boxX = [] + boxY = [] + for j in range(5): + boxX.append(box.get_xdata()[j]) + boxY.append(box.get_ydata()[j]) + boxCoords = zip(boxX,boxY) + + if hasColors: + boxPolygon = Polygon(boxCoords, facecolor = colors[i % len(colors)]) + else: + boxPolygon = Polygon(boxCoords, facecolor = '0.95') + + ax.add_patch(boxPolygon) + return bp + +if __name__=='__main__': + pass + #ensure_dir('../data/DepthMap/test/')