Merge branch 'auto_ground2' into 'main' (79c42a94) · Commits · aflab / astrogeology / Autocnet

autocnet/control/control.py

+30 −13

Original line number	Diff line number	Diff line
		@@ -6,11 +6,12 @@ import geopandas as gpd
		from shapely.geometry import Point

		from autocnet.matcher import subpixel as sp
		from autocnet.transformation.spatial import reproject

		from plio.io.io_controlnetwork import to_isis, write_filelist
		from plio.utils import covariance

		def compute_covariance(df, latsigma, lonsigma, radsigma, radius):
		def compute_covariance(df, dem, latsigma, lonsigma, radsigma):
		"""
		Compute the covariance matrices for constrained or fixed points.

		@@ -19,6 +20,9 @@ def compute_covariance(df, latsigma, lonsigma, radsigma, radius):
		df : pd.DataFrame
		with columns pointtype, adjustedY, and adjustedX

		dem : ~autocnet.spatial.surface.EllipsoidDem or ~autocnet.spatial.surface.GdalDem
		Digital Elevation Model (DEM) object described the target body

		latsigma : int/float
		The estimated sigma (error) in the latitude direction

		@@ -27,19 +31,31 @@ def compute_covariance(df, latsigma, lonsigma, radsigma, radius):

		radsigma : int/float
		The estimated sigma (error) in the radius direction

		radius : int/float
		The body semimajor radius
		"""
		def compute_covar(row, latsigma, lonsigma, radsigma, radius):

		semi_major = dem.a
		semi_minor = dem.c

		def compute_covar(row, latsigma, lonsigma, radsigma, semi_major, semi_minor):
		if row['pointtype'] == 3 or row['pointtype'] == 4:
		return covariance.compute_covariance(row['adjustedY'],
		row['adjustedX'],

		if semi_minor is None:
		semi_minor = semi_major

		x,y,z = row[['aprioriX', 'aprioriY', 'aprioriZ']]
		lon, lat, _ = reproject([x,y,z],
		semi_major, semi_minor,
		'geocent', 'latlon')
		# compute_covariance requires the radius, autocnet operates in height
		radius = dem.get_radius(lat, lon)

		return covariance.compute_covariance(lat,
		lon,
		radius,
		latsigma=latsigma,
		lonsigma=lonsigma,
		radsigma=radsigma,
		semimajor_axis=radius)
		semimajor_axis=semi_major)
		return []

		df['aprioriCovar'] = df.apply(compute_covar,
		@@ -47,7 +63,8 @@ def compute_covariance(df, latsigma, lonsigma, radsigma, radius):
		args=(latsigma,
		lonsigma,
		radsigma,
		radius))
		semi_major,
		semi_minor))
		return df

		def identify_potential_overlaps(cg, cn, overlap=True):

autocnet/graph/network.py

+58 −19

Original line number	Diff line number	Diff line
		@@ -7,6 +7,7 @@ import os
		from shutil import copyfile
		import threading
		from time import gmtime, strftime, time, sleep
		import time
		import logging
		from itertools import combinations

		@@ -1476,8 +1477,8 @@ class CandidateGraph(nx.Graph):
		df["serialnumber"] = serials

		# only populate the new columns for ground points. Otherwise, isis will
		#recalculate the control point lat/lon from control measures which where
		#"massaged" by the phase and template matcher.
		# recalculate the control point lat/lon from the average location of the
		# control measures projection to ground after autocnet matching.
		for i, group in df.groupby('point_id'):
		zero_group = group.iloc[0]
		apriori_geom = np.array(point_info(self.nodes[zero_group.image_index]['data'].geodata.file_name, zero_group.x, zero_group.y, 'image')['BodyFixedCoordinate'].value) * 1000
		@@ -2105,10 +2106,12 @@ class NetworkCandidateGraph(CandidateGraph):

		def to_isis(self,
		path,
		dem = None,
		flistpath=None,
		latsigma=10,
		lonsigma=10,
		radsigma=15,
		ground_xyz=None,
		**db_kwargs):
		"""
		Write a NetworkCandidateGraph to an ISIS control network
		@@ -2118,6 +2121,9 @@ class NetworkCandidateGraph(CandidateGraph):
		path : str
		Outpath to write the control network

		dem : ~autocnet.spatial.surface.EllipsoidDem or ~autocnet.spatial.surface.GdalDem
		Digital Elevation Model (DEM) object described the target body

		flishpath : str
		Outpath to write the associated file list. If None (default),
		the file list is written alongside the control network
		@@ -2131,8 +2137,10 @@ class NetworkCandidateGraph(CandidateGraph):
		radsigma : int/float
		The estimated sigma (error) in the radius direction

		radius : int/float
		The body semimajor radius
		ground_xyz: str
		Path to the file that determined image coordinates of ground points,
		if different than dem argument. This is the file typically used in
		the image registration step of ground points creation.

		db_kwargs : dict
		Kwargs that are passed to the io.db.controlnetwork.db_to_df function
		@@ -2146,15 +2154,26 @@ class NetworkCandidateGraph(CandidateGraph):
		of paths to the images being included in the control network

		"""

		if dem is None:
		dem_file = None
		log.warning(f'No dem argument passed; covariance matrices will be computed for the points.')
		else:
		if isinstance(dem, EllipsoidDem): # not sure about this
		dem_file = f'EllipsoidDem a={dem.a} b {dem.b} c={dem.c}'
		elif isinstance(dem, GdalDem):
		dem_file = dem.dem.file_name

		# Read the cnet from the db
		df = io_controlnetwork.db_to_df(self.engine, **db_kwargs)
		df = io_controlnetwork.db_to_df(self, ground_radius=dem_file, ground_xyz=ground_xyz, **db_kwargs)

		# Add the covariance matrices to ground measures
		if dem is not None:
		df = control.compute_covariance(df,
		dem,
		latsigma,
		lonsigma,
		radsigma,
		self.config['spatial']['semimajor_rad'])
		radsigma)

		if flistpath is None:
		flistpath = os.path.splitext(path)[0] + '.lis'
		@@ -2241,7 +2260,7 @@ class NetworkCandidateGraph(CandidateGraph):

		return obj

		def add_from_filelist(self, filelist, clear_db=False):
		def add_from_filelist(self, filelist, clear_db=False, exist_check=False):
		"""
		Parse a filelist to add nodes to the database.

		@@ -2265,17 +2284,37 @@ class NetworkCandidateGraph(CandidateGraph):
		self.clear_db()

		total=len(filelist)
		db_images = []
		for cnt, f in enumerate(filelist):
		# Create the nodes in the graph. Really, this is creating the
		# images in the DB
		log.info('loading {} of {}'.format(cnt+1, total))
		self.add_image(f)

		image_name = os.path.basename(f)
		node = NetworkNode(image_path=f, image_name=image_name)
		node.parent = self

		i = node.create_db_element(exist_check=exist_check)
		db_images.append(i)
		break

		if cnt%1000 == 0:
		log.info('uploading 1000 images to database...')
		with self.session_scope() as session:
		session.add_all(db_images)
		db_images = []

		with self.session_scope() as session:
		log.info(f'uploading final {len(db_images)} images to database.')
		session.add_all(db_images)



		self.from_database()
		# Execute the computation to compute overlapping geometries
		self._execute_sql(sql.compute_overlaps_sql)

		def add_image(self, img_path):
		def add_image(self, img_path, exist_check=True, add_keypoints=False):
		"""
		Upload a single image to NetworkCandidateGraph associated DB.

		@@ -2292,7 +2331,7 @@ class NetworkCandidateGraph(CandidateGraph):
		image_name = os.path.basename(img_path)
		node = NetworkNode(image_path=img_path, image_name=image_name)
		node.parent = self
		node.populate_db()
		node.populate_db(exist_check=exist_check, add_keypoints=add_keypoints)
		return node['node_id']

		def copy_images(self, newdir):
		@@ -2771,7 +2810,7 @@ class NetworkCandidateGraph(CandidateGraph):
		# return missing image id pairs
		return [e for e in all_edges if e not in graph.edges]

		def cluster_propagate_control_network(self,
		def cluster_propagate_ground_points(self,
		base_cnet,
		walltime='00:20:00',
		chunksize=1000,
		@@ -2910,8 +2949,8 @@ class NetworkCandidateGraph(CandidateGraph):
		def subpixel_regiter_mearure(self, measureid, **kwargs):
		subpixel.subpixel_register_measure(self.Session, measureid, **kwargs)

		def propagate_control_network(self, control_net, **kwargs):
		cim.propagate_control_network(self.Session,
		def propagate_ground_points(self, control_net, **kwargs):
		cim.propagate_ground_points(self.Session,
		self.config,
		self.dem,
		control_net)

autocnet/graph/node.py

+73 −34

Original line number	Diff line number	Diff line
		@@ -3,6 +3,8 @@ import itertools
		import os
		import logging

		import time

		from csmapi import csmapi
		import numpy as np
		import pandas as pd
		@@ -507,7 +509,9 @@ class NetworkNode(Node):
		# If this is the first time that the image is seen, add it to the DB
		self.job_status = defaultdict(dict)

		def populate_db(self):

		def create_db_element(self, exist_check=False, add_keypoints=False):
		if exist_check:
		with self.parent.session_scope() as session:
		res = session.query(Images).filter(Images.path == self['image_path']).first()
		if res:
		@@ -516,18 +520,19 @@ class NetworkNode(Node):

		# If the geodata is not valid, do no create an assocaited keypoints file
		# One instance when invalid is during testing.
		kps = None
		if add_keypoints:
		if hasattr(self.geodata, 'file_name'):
		kpspath = io_keypoints.create_output_path(self.geodata.file_name)
		# Create the keypoints entry
		kps = Keypoints(path=kpspath, nkeypoints=0)
		else:
		kps = None

		try:
		fp, cam_type = self.footprint
		fp, cam_type = self.generate_footprint(exist_check=exist_check)
		except Exception as e:
		log.warning('Unable to generate image footprint.\n{}'.format(e))
		fp = cam_type = None

		# Create the image
		i = Images(name=self['image_name'],
		path=self['image_path'],
		@@ -537,6 +542,11 @@ class NetworkNode(Node):
		serial=self.isis_serial,
		cam_type=cam_type)

		return i

		def populate_db(self, exist_check=False, add_keypoints=False):
		i = self.create_db_element(exist_check=exist_check, add_keypoints=add_keypoints)

		with self.parent.session_scope() as session:
		session.add(i)

		@@ -663,22 +673,25 @@ class NetworkNode(Node):
		self._camera = plugin.constructModelFromState(res.camera)
		return self._camera

		@property
		def footprint(self):

		def footprint_from_database(self):
		with self.parent.session_scope() as session:
		res = session.query(Images).filter(Images.id == self['node_id']).first()

		# not in database, create footprint
		if res is None:
		# get ISIS footprint if possible
		if utils.find_in_dict(self.geodata.metadata, "Polygon"):
		return None, None
		else:
		footprint_latlon = res.geom
		cam_type = res.cam_type
		return footprint_latlon, cam_type

		def footprint_from_isis(self):
		footprint_latlon = shapely.wkt.loads(self.geodata.footprint.ExportToWkt())
		if isinstance(footprint_latlon, shapely.geometry.Polygon):
		footprint_latlon = shapely.geometry.MultiPolygon(list(footprint_latlon))
		cam_type = 'isis'
		return footprint_latlon, cam_type
		# Get CSM footprint
		else:

		def footprint_from_csm(self):
		boundary = generate_boundary(self.geodata.raster_size[::-1]) # yx to xy
		footprint_latlon = generate_latlon_footprint(self.camera,
		boundary,
		@@ -686,10 +699,36 @@ class NetworkNode(Node):
		footprint_latlon.FlattenTo2D()
		cam_type = 'csm'
		return footprint_latlon, cam_type

		def generate_footprint(self, exist_check=True):
		footprint_latlon = cam_type = None

		if exist_check:
		print('ERRORRRORORORORRR!!!')
		footprint_latlon, cam_type = self.footprint_from_database()

		# not in database, create footprint
		if footprint_latlon is None and cam_type is None:
		if utils.find_in_dict(self.geodata.metadata, "Polygon"):
		# get ISIS footprint if possible
		t1 = time.time()
		footprint_latlon, cam_type = self.footprint_from_isis()
		t2 = time.time()
		print(f'time to get ISIS footprint: {t2 - t1}')
		else:
		# in database, return footprint
		footprint_latlon = res.footprint_latlon
		return footprint_latlon
		# Get CSM footprint
		t1 = time.time()
		footprint_latlon, cam_type = self.footprint_from_csm()
		t2 = time.time()
		print(f'time to get CSM footprint: {t2 - t1}')

		return footprint_latlon, cam_type

		@property
		def footprint(self):
		if not hasattr(self, '_footprint'):
		self._footprint = self.generate_footprint()
		return self._footprint

		@property
		def points(self):

autocnet/io/db/controlnetwork.py

+35 −21

Original line number	Diff line number	Diff line
		@@ -9,20 +9,32 @@ from autocnet.io.db.model import Measures
		from autocnet.spatial.isis import isis2np_types
		from ... import sql

		def db_to_df(engine, sql = sql.db_to_df_sql_string):
		def db_to_df(ncg, ground_radius=None, ground_xyz=None, sql=sql.db_to_df_sql_string):
		"""
		Given a set of points/measures in an autocnet database, generate an ISIS
		compliant control network.

		Parameters
		----------
		engine : Object
		sqlalchemy engine object to read from
		ncg : Object
		A NetworkCandidateGraph object connected to the target database.

		ground_radius : str
		Description of file used to generate radius values. This value
		can take the form of a file path to a DEM or a general 'EllipsoidDem'.

		ground_xyz: str
		Path to the file that determined image coordinates of ground points,
		if different than dem argument. This is the file typically used in
		the image registration step of ground points creation.

		sql : str
		The sql query to execute in the database.
		The sql query to execute in the database. Default grabs information
		for all non-ignored (including jigsaw ignored) points and measures
		that exist on images with 3 or more measures.
		"""
		df = pd.read_sql(sql, engine)

		df = pd.read_sql(sql, ncg.engine)

		# measures.id DB column was read in to ensure the proper ordering of DF
		# so the correct measure is written as reference
		@@ -37,26 +49,28 @@ def db_to_df(engine, sql = sql.db_to_df_sql_string):
		df['aprioriX'] = 0
		df['aprioriY'] = 0
		df['aprioriZ'] = 0
		df['adjustedX'] = 0
		df['adjustedY'] = 0
		df['adjustedZ'] = 0
		df['aprioriCovar'] = [[] for _ in range(len(df))]
		df['aprioriSurfPointSource'] = 0
		df['aprioriSurfPointSourceFile'] = None
		df['aprioriRadiusSource'] = 0
		df['aprioriRadiusSourceFile'] = None

		# only populate the new columns for ground points. Otherwise, isis will
		#recalculate the control point lat/lon from control measures which where
		#"massaged" by the phase and template matcher.
		# recalculate the control point lat/lon from the average location of the
		# control measures projection to ground after autocnet matching.
		for i, row in df.iterrows():
		if row['pointtype'] == 3 or row['pointtype'] == 4:
		if row['apriori']:
		apriori_geom = swkb.loads(row['apriori'], hex=True)
		row['aprioriX'] = apriori_geom.x
		row['aprioriY'] = apriori_geom.y
		row['aprioriZ'] = apriori_geom.z
		if row['adjusted']:
		adjusted_geom = swkb.loads(row['adjusted'], hex=True)
		row['adjustedX'] = adjusted_geom.x
		row['adjustedY'] = adjusted_geom.y
		row['adjustedZ'] = adjusted_geom.z
		row['aprioriZ'] = apriori_geom.z # this is a height measurement
		if ground_radius is not None:
		row['aprioriRadiusSource'] = 5 # corresponds to DEM in plio AprioriSource protobuf Enum
		row['aprioriRadiusSourceFile'] = ground_radius
		if ground_xyz is not None:
		row['aprioriSurfPointSource'] = 6 # corresponds to Basemap in plio AprioriSource protobuf Enum
		row['aprioriSurfPointSourceFile'] = ground_xyz
		df.iloc[i] = row

		return df

autocnet/matcher/cpu_extractor.py

+1 −0

Original line number	Diff line number	Diff line
		@@ -111,6 +111,7 @@ def extract_most_interesting(image, extractor_method='orb', extractor_parameters
		The keypoints row with the higest variance. The row has 'x' and 'y' columns to
		get the location.
		"""

		kps, desc = extract_features(image,
		extractor_method=extractor_method,
		extractor_parameters=extractor_parameters)