1#!/usr/local/bin/python3.8 2# -*- coding: utf-8 -*- 3 4gpl = r""" 5 woa.py - warped overlap analysis 6 generate control points from warped overlap images 7 8 Copyright (C) 2011 Kay F. Jahnke 9 10 This program is free software: you can redistribute it and/or modify 11 it under the terms of the GNU General Public License as published by 12 the Free Software Foundation, either version 3 of the License, or 13 (at your option) any later version. 14 15 This program is distributed in the hope that it will be useful, 16 but WITHOUT ANY WARRANTY; without even the implied warranty of 17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18 GNU General Public License for more details. 19 20 You should have received a copy of the GNU General Public License 21 along with this program. If not, see <http://www.gnu.org/licenses/>. 22""" 23# @category Control Points 24# @name Warped Overlap Analysis 25# @api-min 2014.0 26# @api-max 2019.2 27 28# note that if you want to read the script, it's written bottom-up, so the 29# higher-level routines are towards the end. 30 31import hsi 32import os 33import sys 34import traceback 35import re 36import copy 37import math 38import itertools 39import subprocess 40 41# our very own exception: 42 43class WoaError ( Exception ) : 44 45 def __init__ ( self , text ) : 46 self.text = text 47 48 def __str__ ( self ) : 49 return self.text 50 51# verbose-print: print only if args.verbose is set 52 53def vpr ( *x ) : 54 if args.verbose : 55 print ( *x ) 56 57# we keep a global variable for arguments - if called from the command 58# line the argument parser will overwrite it, but since a plugin 59# can't receive arguments, we use a different approach then: if there 60# is a file 'woa.ini' in the user's home directory, the arguments are read 61# from it; otherwise sensible defaults are used. 62# Since the mechanism is only there to bridge the time until a GUI for 63# plugins becomes reality, it's sloppily programmed and does not check 64# the woa.ini file very much. 65# hugin's pwd is /tmp, which is no good location to look for woa.ini. 66# That's why it's expected in the user's home directory. 67 68class argobj : 69 70 def __init__ ( self , inifile = '~/woa.ini' ) : 71 72 self.focus = None 73 74 # if woa.ini can be found, read arguments from it 75 76 path = os.path.expanduser ( inifile ) 77 if os.path.exists ( path ) : 78 79 import ConfigParser 80 f = open ( path , 'r' ) 81 config = ConfigParser.ConfigParser() 82 config.readfp ( f ) 83 84 self.basename = config.get ( 'woa arguments' , 'basename' ) 85 self.ceiling = config.getfloat ( 'woa arguments' , 'ceiling' ) 86 self.cpg = config.get ( 'woa arguments' , 'cpg' ) 87 self.margin = config.getint ( 'woa arguments' , 'margin' ) 88 self.prolific = config.getboolean ( 'woa arguments' , 'prolific' ) 89 self.scale = config.getfloat ( 'woa arguments' , 'scale' ) 90 self.stop = config.get ( 'woa arguments' , 'stop' ) 91 self.threshold = config.getfloat ( 'woa arguments' , 'threshold' ) 92 self.verbose = config.getboolean ( 'woa arguments' , 'verbose' ) 93 94 else : 95 96 # no woa.ini - use sensible defaults 97 98 self.basename = 'woa' 99 self.ceiling = 1.0 100 self.cpg = 'cpfind' # changed from 'autopano-sift-c' KFJ 2011-06-15 101 self.margin = 0 102 self.prolific = False 103 self.scale = 0.25 104 self.stop = None 105 self.threshold = 0.01 106 self.verbose = False 107 108# we want to use this variable globally: 109 110args = None 111 112# require_program is used to check if a required helper program 113# can be used. The test is simple: try call it; if that fails, 114# raise a run time error which terminates woa. 115# command is the program's name. The output is sent to a pipe 116# but the pipe is never read, so it's effectively muted. 117 118def require_program ( command ) : 119 120 try : 121 122 vpr ( 'checking for availability of %s' % command ) 123 subject = subprocess.Popen ( args = command , 124 stdout = subprocess.PIPE , 125 stderr = subprocess.PIPE ) 126 ignore = subject.communicate() 127 128 except : 129 130 err = "can't find required program %s" % command 131 raise WoaError ( err ) 132 133# check_helpers makes sure all needed helper programs are available. 134# Note that this has to be called after the global args object is 135# initialized, because it is only known then which cpg is used. 136 137def check_helpers() : 138 139 require_program ( 'tiffdump' ) 140 require_program ( 'nona' ) 141 require_program ( args.cpg ) 142 143 144# when comparing positions, if they differ by at most POSITION_DELTA, 145# they will be considered equal 146 147POSITION_DELTA = 0.000001 148 149# I'm using a cartesian/polar coordinate transform modified so 150# it works with hugin's lon/lat system and degrees. The transform is 151# into a right handed system with (0,0) transforming to (1,0,0), z up 152 153def cartesian3_to_polar2 ( p3 ) : 154 155 azimuth = math.atan2 ( p3[1] , p3[0] ) 156 inclination = math.acos ( p3[2] ) 157 158 return ( math.degrees ( azimuth ) , 159 90.0 - math.degrees ( inclination ) ) 160 161# and back 162 163def polar2_to_cartesian3 ( p2 ) : 164 165 azimuth = math.radians ( p2.x ) 166 inclination = math.radians ( 90.0 - p2.y ) 167 x = math.sin ( inclination ) * math.cos ( azimuth ) 168 y = math.sin ( inclination ) * math.sin ( azimuth ) 169 z = math.cos ( inclination ) 170 171 return ( x , y , z ) 172 173def normalize ( cart ) : 174 175 sqmag = cart[0] * cart[0] + cart[1] * cart[1] + cart[2] * cart[2] 176 mag = math.sqrt ( sqmag ) 177 rv = ( cart[0] / mag , cart[1] / mag , cart[2] / mag ) 178 return rv 179 180# KFJ 2012-03-17 common subroutine for calling helper 181# programs. The helper's output is captured and returned 182# to the caller together with the return code. 183# per default, the command's output is not displayed, 184# which differs from the previous behaviour, where 185# the CPG's greeting lines etc. were littering the 186# output too much for my taste. 187 188def run_helper_program ( args , mute = True ) : 189 190 # vpr ( 'calling %s %s' % ( args[0] , args[1:] ) ) 191 192 try : 193 p = subprocess.Popen ( args , 194 stdout = subprocess.PIPE , 195 stderr = subprocess.STDOUT) 196 except OSError as problem : 197 message = ( 'OS error: "%s"\nwhen trying to execute:\n %s\n' % 198 ( problem.strerror , args ) ) 199 raise WoaError ( message ) 200 201 p.wait() 202 output=p.stdout.readlines() 203 204 if not mute: 205 for line in output: 206 print ( line.strip() ) 207 208 # vpr ( '%s returned %d' % ( args[0] , p.returncode ) ) 209 210 return p.returncode, output 211 212# get_tiff_offset is needed to find the offsets in cropped TIFFS. 213# This is done with a call to tiffdump, a utility program that comes 214# with libtiff. I'm not entirely happy with this solution, but I found 215# no more straightforward way to extract the data. 216# Basially, tiffdump is called, it's output is scanned into a dictionary 217# and the needed values to calculate the offsets are taken from the 218# dictionary. A tuple with xoffset and yoffset in pixels is returned. 219# If it's not a 'cropped TIFF' the offsets are returned as 0,0 220 221def get_tiff_offset ( image_file ) : 222 223 224 # KFJ 2012-03-17 now using run_helper_program() to execute tiffdump 225 226 returncode, output = run_helper_program ( [ 'tiffdump' , image_file ] ) 227 228 # to find the offset values, we have to pick apart tiffdump's 229 # console output: 230 231 gleaned = dict() 232 233 for tag in output : 234 t = tag.strip() 235 encoding = 'utf-8' 236 t = t.decode(encoding) 237 fieldname = re.match ( r'[^\(]+' , t ) . group(0) . strip() 238 fieldcontent = re.match ( r'([^<]+<)([^>]+)' , t ) 239 if fieldcontent: 240 gleaned [ fieldname ] = fieldcontent.group(2) 241 242 # nothing found at all? 243 244 if not gleaned : 245 err = 'could not find any TIFF information in %s' % image_file 246 print ( err ) 247 raise WoaError ( err ) 248 249 # we try to read out the XPosition and YPosition fields in 250 # the TIFF file, If they aren't set, it's simply not a 'cropped 251 # TIFF' and we take them to be zero. 252 253 xpos = float ( gleaned.get ( 'XPosition' , 0 ) ) 254 ypos = float ( gleaned.get ( 'YPosition' , 0 ) ) 255 xres = float ( gleaned.get ( 'XResolution' ) ) 256 yres = float ( gleaned.get ( 'YResolution' ) ) 257 258 xoff = xpos * xres 259 yoff = ypos * yres 260 261 return xoff , yoff 262 263# some global variables to reduce overhead - we just reuse the same 264# FDiff2D objects every time we call the transforms 265 266transform_in = hsi.FDiff2D(0,0) 267transform_inter = hsi.FDiff2D(0,0) 268transform_out = hsi.FDiff2D(0,0) 269 270# transform point coordinates using panotools transforms 271# using the global variables above for source and target 272 273def transform_point ( p , tf ) : 274 275 transform_in.x = p.x 276 transform_in.y = p.y 277 278 if tf.transform ( transform_out , transform_in ) : 279 return point ( transform_out.x , transform_out.y ) 280 281 vpr ( 'transform failed for %s' % str ( p ) ) 282 return None 283 284# when transforming from one image to the other, two transforms 285# are needed: img0 to pano and pano to img1 286 287def double_transform_point ( p , tf1 , tf2 ) : 288 289 transform_in.x = p.x 290 transform_in.y = p.y 291 292 if tf1.transform ( transform_inter , transform_in ) : 293 if tf2.transform ( transform_out , transform_inter ) : 294 return point ( transform_out.x , transform_out.y ) 295 296 vpr ( 'transform failed for %s' % str ( p ) ) 297 return None 298 299# next is a set of stripped-down geometrical primitives that are just 300# sufficient to do the maths we need for masking the overlapping regions 301# of the images 302 303class point : 304 305 def __init__ ( self , ix , iy ) : 306 307 self.x = ix 308 self.y = iy 309 310 def distance ( self , other ) : 311 312 dx = other.x - self.x 313 dy = other.y - self.y 314 return math.sqrt ( dx * dx + dy * dy ) 315 316 def delta ( self , other ) : 317 318 dx = other.x - self.x 319 dy = other.y - self.y 320 return ( dx , dy ) 321 322 def __str__ ( self ) : 323 324 return '%s (%.5f, %.5f)' % ( id(self) , self.x , self.y ) 325 326# a segment contains it's two end points and some auxilliary data 327# to facilitate mathematics on the segment. Note that our segment 328# notion is planar and uses image coordinates. 329 330class segment : 331 332 # we keep some more state in the segment object which is needed 333 # repeatedly in various places 334 335 def __init__ ( self , ia , ib ) : 336 337 self.a = ia 338 self.b = ib 339 340 self.length = ia.distance ( ib ) 341 342 if ia.x == ib.x : 343 self.minx = self.maxx = ia.x 344 self.dx = 0.0 345 else : 346 self.minx = min ( ia.x , ib.x ) 347 self.maxx = max ( ia.x , ib.x ) 348 self.dx = ib.x - ia.x 349 if ia.y == ib.y : 350 self.miny = self.maxy = ia.y 351 self.dy = 0.0 352 else : 353 self.miny = min ( ia.y , ib.y ) 354 self.maxy = max ( ia.y , ib.y ) 355 self.dy = ib.y - ia.y 356 357 def __str__ ( self ) : 358 359 return '[ %s %s ]' % ( self.a , self.b ) 360 361 # we detect only 'true' intersections, so if the segment 362 # lies fully on the horizontal/vertical, this is not considered 363 # an intersection. 364 365 def intersect_with_vertical ( self , x ) : 366 367 if abs ( self.dx ) > POSITION_DELTA : 368 if self.minx <= x <= self.maxx : 369 return ( self.a.y 370 + ( x - self.a.x ) 371 * ( self.b.y - self.a.y ) 372 / ( self.b.x - self.a.x ) ) 373 374 return None 375 376 def intersect_with_horizontal ( self , y ) : 377 378 if abs ( self.dy ) > POSITION_DELTA : 379 if self.miny <= y <= self.maxy : 380 return ( self.a.x 381 + ( y - self.a.y ) 382 * ( self.b.x - self.a.x ) 383 / ( self.b.y - self.a.y ) ) 384 385 return None 386 387 # segments are used with z-values that range from 0 (at point a) 388 # to 1 (at point b). Locate will yield the point corresponding to z. 389 390 def locate ( self , z ) : 391 392 assert 0.0 <= z <= 1.0 393 394 if z == 0.0 : 395 return self.a 396 elif z == 1.0 : 397 return self.b 398 else : 399 return point ( self.a.x + z * self.dx , self.a.y + z * self.dy ) 400 401 # contains tests if point p is on this segment. 402 # a certain delta is deemed acceptable 403 404 def contains ( self , p , delta = POSITION_DELTA ) : 405 406 # our test works like this: 407 # all points on the line segment satisfy L = A + z ( dx , dy ) 408 # with 0 <= z <= 1 409 # we assume P is on L. So we calculate a z for it and test if 410 # 1. 0 <= z <= 1 411 # 2. P = A + z ( dx , dy ) 412 # if dx has greater magnitude than dy, we calculate z by way of 413 # the horizontal distance, otherwise by way of the vertical distance. 414 # We assume that the case where dx == dy == 0 never occurs. 415 416 if abs ( self.dx ) >= abs ( self.dy ) : 417 # horizontal distance between p and self.a 418 dpx = p.x - self.a.x 419 # in ratio to self.dx yields z 420 z = dpx / self.dx 421 # if z is roughly between 0 and 1, the match is possible 422 if - delta <= z <= 1.0 + delta : 423 # going z * dy in vertical direction from self.a 424 # yields a point on the line segment 425 yy = self.a.y + z * self.dy 426 if abs ( yy - p.y ) < delta : 427 # if this value coincides with p.y, we've come 428 # full circle and the point must be on the segment. 429 # print ( self , 'contains' , p , z ) 430 return True 431 else : 432 # the symmetric case, where we start out with the vertical 433 # distance. 434 dpy = p.y - self.a.y 435 z = dpy / self.dy 436 if - delta <= z <= 1.0 + delta : 437 xx = self.a.x + z * self.dx 438 if abs ( xx - p.x ) < delta : 439 # print ( self , 'contains' , p , z ) 440 return True 441 442 # print ( self , 'does not contain' , p ) 443 return False 444 445# routine to make a list of corner points from a SrcPanoImage object. 446# The point list can be used to create a polygon object. 447# We also pass back left, top, right and bottom for good measure. 448 449def img_to_point_list ( img ) : 450 451 w = img.getWidth() 452 h = img.getHeight() 453 454 # note that we're in center-relative coordinates, so this symmetric 455 # setup is subpixel-correct. 456 457 left = - w / 2.0 - args.margin 458 bottom = - h / 2.0 - args.margin 459 right = w / 2.0 + args.margin 460 top = h / 2.0 + args.margin 461 462 return ( [ point ( left , top ) , 463 point ( right , top ) , 464 point ( right , bottom ) , 465 point ( left , bottom ) ] , 466 left , top , right , bottom ) 467 468# class polygon contains a set of segments. These segments are meant 469# to form a closed polygon, but this is not checked or enforced. 470# The implementation is minimal and geared towards the intended use, 471# this is not library-grade code. 472 473class polygon : 474 475 def __init__ ( self , 476 ipoints = [] , 477 reverse = False , 478 close = True ) : 479 480 # the points may be used in reverse order 481 482 if reverse : 483 ipoints = [ p for p in reversed ( ipoints ) ] 484 485 # keep number of segments as state, this has to be updated 486 # if new segments are taken in. A sacrifice to efficiency 487 # instead of calling len() always. 488 489 self.count = len ( ipoints ) 490 491 # the polygon is closed; the last segment 492 # joins last and first point 493 494 self.segments = [] 495 if self.count >= 2 : 496 a = ipoints[0] 497 for p in ipoints[1:] : 498 self.segments.append ( segment ( a , p ) ) 499 a = p 500 if close : 501 self.segments.append ( segment ( a , ipoints[0] ) ) 502 503 def __str__ ( self ) : 504 505 result = '' 506 nr = 0 507 for s in self.segments : 508 result += 'SEG %02d %s\n' % ( nr , s ) 509 nr += 1 510 511 return result 512 513 # take in a point. This point has to be on one of the segments 514 # and the effect is to split the segment in two where the point is. 515 # if this fails, False is retuned, otherwise True. 516 517 def take_in ( self , p ) : 518 519 new_seglist = [] 520 inserted = False 521 for s in self.segments : 522 if ( not inserted ) and s.contains ( p ) : 523 inserted = True 524 new_seglist.append ( segment ( s.a , p ) ) 525 new_seglist.append ( segment ( p , s.b ) ) 526 else : 527 new_seglist.append ( s ) 528 if inserted : 529 self.segments = new_seglist 530 self.count = len ( new_seglist ) 531 532 return inserted 533 534# class dot acts as a cursor (or iterator) on a polygon. It will 535# sample the polygon circumference, producing sample points every 536# so often, and will also stop at all corner points. 537 538class dot : 539 540 def __init__ ( self , ipolygon , 541 istage = 0 , iwind = 0 , 542 istride = 0.1 , iz = 0.0 , 543 idelta = POSITION_DELTA ) : 544 545 self.track = ipolygon 546 self.stage = istage 547 self.wind = iwind 548 self.stride = istride 549 self.z = iz 550 self.delta = idelta 551 self.current_segment = self.track.segments [ self.stage ] 552 553 # z is the normalized distance from the beginning of the current 554 # segment to the current sample location: 555 556 def position ( self ) : 557 558 return self.current_segment.locate ( self.z ) 559 560 def __str__ ( self ) : 561 562 return 'dot track %d stage %d wind %d z %f %s' % ( 563 id ( self.track ) , 564 self.stage , 565 self.wind , 566 self.z , 567 self.position() ) 568 569 # advance to next sample location, or to next corner point 570 # if that is nearer 571 572 def advance ( self ) : 573 574 self.z += self.stride / self.current_segment.length 575 576 # if z is at (within delta) or past the end of the current 577 # segment, advance to the beginning of the next segment. 578 579 if self.z + self.delta >= 1.0 : 580 self.stage += 1 581 if self.stage >= self.track.count : 582 self.stage = 0 583 self.wind += 1 584 self.current_segment = self.track.segments [ self.stage ] 585 self.z = 0.0 586 587 # termination criterion. If the dot is copied at the beginning 588 # of a circumambulation to 'start', passed ( start ) will yield 589 # True if the start point has been passed again 590 591 def passed ( self , start ) : 592 593 if self.wind > start.wind : 594 if self.stage < start.stage : 595 return False 596 elif self.stage == start.stage : 597 if self.z <= start.z : 598 return False 599 return True 600 return False 601 602# shift an initial part of the outline to it's end 603# so that the outline starts with an 'in' type crossover point 604 605def in_point_to_tip ( outline ) : 606 607 head = [] 608 tail = iter ( outline ) 609 for p in tail : 610 if p.direction == 'in' : 611 outline = [ p ] 612 break 613 head.append ( p ) 614 615 outline.extend ( tail ) 616 outline.extend ( head ) 617 618 return outline 619 620# produce a set of include and exclude masks depending on image overlap. 621# This routine is central to the script, and most code up to here serves 622# it. The masking of the overlapping regions is non-trivial: 623# A proper mathematical description like y=f(x) is only available for 624# either of the image margins in coordinates of the image itself, where 625# these are trivially horizontals and verticals. When transformed to pano 626# space or even to the other image's coordinates, these straight lines become 627# curves, and their shape depends on many factors: image size and orientation, 628# lens correction coefficients, projection types. The result of all these 629# factors can be obtained for a single point at a time with some precision 630# by using the panotools transforms, but this all we have. To match the 631# curved lines, iterations have to be used, and for some points in the 'other' 632# image there may not even exist a valid transform to coordinates of 'this' 633# image. I have taken the approach to move only 'inside' 'this' image and 634# avoid looking at points 'outside' apart from their status as being 'outside'. 635# The resulting code does the trick, but I haven't tested all corner cases - 636# I'm slightly concerned with images that have touching opposite margins 637# (e.g. 360X180 equirects) - the code may break here, I'll have to test. 638 639def mask_nonoverlaps ( pano , img0 , img1 ) : 640 641 # for consistency throughout the script, 642 # variables ending with '0' pertain to image 0, 643 # which is also refered to as 'this' image 644 # variables ending in '1' to image 1 645 # which is also refered to as the 'other' image. 646 # first we have to get the image metrics and the 647 # transforms 648 649 pano_options = pano.getOptions() 650 651 # we create transforms from pano coordinates to image coordinates 652 653 tf0 = hsi.Transform() 654 tf0.createTransform ( img0 , pano_options ) 655 656 #tf1 = hsi.Transform() 657 #tf1.createTransform ( img1 , pano_options ) 658 659 # and the reverse transform, image to pano 660 661 #rtf0 = hsi.Transform() 662 #rtf0.createInvTransform ( img0 , pano_options ) 663 664 rtf1 = hsi.Transform() 665 rtf1.createInvTransform ( img1 , pano_options ) 666 667 # we also generate a polygon from 'this' image and the 'other' 668 # image as well. 669 670 points , left , top , right , bottom = img_to_point_list ( img0 ) 671 672 # later on we'll insert crosover points into this polygon with 673 # direction 'in' or 'out', but the points we put in now are 674 # not crosover points and have no 'direction' 675 676 for p in points : 677 p.direction = None 678 679 p0 = polygon ( points ) 680 681 points , left1 , top1 , right1 , bottom1 = img_to_point_list ( img1 ) 682 p1 = polygon ( points ) 683 684 # we will need 'this' image's true extent and the x and y offset 685 # from the center 686 687 w0 = img0.getWidth() 688 h0 = img0.getHeight() 689 690 # note, again, we're using center-realtive symmetric coordinates 691 692 true_left = - w0 / 2.0 693 true_bottom = - h0 / 2.0 694 695 # the next section defines a bunch of local functions that 696 # use the current metrics. This reduces parameter passing 697 # and makes the code easier to read. 698 699 # min_distance_from_margin will calculate the distance from the 700 # nearest image boundary. 701 702 def min_distance_from_margin ( p ) : 703 704 return min ( abs ( p.x - left ) , 705 abs ( right - p.x ) , 706 abs ( p.y - bottom ) , 707 abs ( top - p.y ) ) 708 709 # 'inside' criterion: True if point p is within this images's 710 # boundaries. Note that here the routine works on points in 711 # 'this' image's coordinates. 712 # We return the minimal distance from the margin to allow 713 # for smaller strides near the edge. Negative proximity 714 # signals a failed transform. Points very near the edge 715 # are taken to be inside. 716 717 def inside ( p ) : 718 719 if not p : 720 return False , -1.0 721 722 if ( ( left <= p.x <= right ) 723 and ( bottom <= p.y <= top ) ) : 724 return True , min_distance_from_margin ( p ) 725 726 md = min_distance_from_margin ( p ) 727 if md <= POSITION_DELTA : 728 return True , 0.0 729 else : 730 return False , md 731 732 # the double projection from one image to another is quite imprecise. 733 # We need the crossover point to be precisely on 'this' image's boundary 734 # for the inclusion in take_in to succeed. 735 736 # initially I was using a threshold of .001 in the comparisons 737 # below, but this didn't work sometimes. I've now increased the 738 # threshold to .01, which seems to work, but TODO: I'd like to 739 # figure out a better way KFJ 2011-06-13 740 # lowered further from .01 to .1 KFJ 2011-09-15 741 742 def put_on_margin ( p ) : 743 744 x = None 745 y = None 746 747 if abs ( p.x - left ) <= 0.1 : 748 x = left 749 y = p.y 750 elif abs ( p.x - right ) <= 0.1 : 751 x = right 752 y = p.y 753 754 if abs ( p.y - top ) <= 0.1 : 755 y = top 756 x = p.x 757 elif abs ( p.y - bottom ) <= 0.1 : 758 y = bottom 759 x = p.x 760 761 if x is None or y is None : 762 raise WoaError ( 'cannot put %s on margin' % p ) 763 764 return point ( x , y ) 765 766 # project_in yields a point from the other image 767 # projected into this one 768 769 def project_in ( p ) : 770 771 return double_transform_point ( p , rtf1 , tf0 ) 772 773 # localize_crossing uses an iteration to approach the point where the 774 # contours cross. This is superior to intersecting a line segment between 775 # two margin points from the 'other' image. it approaches the true 776 # intersection point to arbitrary precision - within the precision 777 # constraints of the double panotools transform. 778 779 def localize_crossing ( cseg ) : 780 781 # first, find the point that is inside 782 if cseg.a.inside : 783 pin = cseg.a 784 pout = cseg.b 785 direction = 'out' 786 else : 787 pin = cseg.b 788 pout = cseg.a 789 direction = 'in' 790 791 stride = 0.5 792 delta = pin.delta ( pout ) 793 794 # we go on from the current position, but if we've 795 # gone too far, we omit the step. Then we repeat with 796 # half the step, and so on. 797 798 while stride >= POSITION_DELTA : 799 candidate = point ( pin.x + stride * delta[0] , 800 pin.y + stride * delta[1] ) 801 twin = project_in ( candidate ) 802 ins , md = inside ( twin ) 803 if ins : 804 pin = candidate 805 stride /= 2.0 806 807 # pin is now very near the actual crossing, and in coordinates 808 # of the other image. We return this point as the crossing point. 809 810 pin.direction = direction 811 pin.twin = project_in ( pin ) 812 813 # we need to locate pin.twin precisely on the image boundary, 814 # so that the insertion routine succeeds. 815 816 pin.twin = put_on_margin ( pin.twin ) 817 818 # vpr ( 'crossover i0:' , pin , 'i1:' , pin.twin ) 819 pin.inside = True 820 pin.proximity = 0.0 821 return pin 822 823 # now we're set up and start the routine proper: 824 # we already have p1, a polygon from the 'other' image's margin, 825 # and we sample it with standard_stride pixel steps. 826 827 standard_stride = 100.0 828 d = dot ( p1 , istride = standard_stride ) 829 outline = [] 830 total_coincidence = True 831 832 # we walk the whole contour of the other image, saving each 833 # point along the way, amended with additional information 834 835 start = copy.copy ( d ) 836 while not d.passed ( start ) : 837 # get position of point (in other image's coordinates) 838 p = d.position() 839 # try getting the corresponding coordinates in this image 840 p.twin = project_in ( p ) 841 # note that p.twin may be None if the transform failed, 842 # in which case inside() returns False, 0.0 843 p.inside , p.proximity = inside ( p.twin ) 844 # if any of the points are further than POSITION_DELTA 845 # away from this image's boundary, it's not a total coincidence. 846 if p.proximity >= POSITION_DELTA : 847 total_coincidence = False 848 p.direction = None 849 # here we might modify the stride depending on proximity 850 # for now we just carry on with the same stride 851 outline.append ( p ) 852 d.advance() 853 854 if total_coincidence : 855 vpr ( 'images boundaries coincide totally' ) 856 points , left , top , right , bottom = img_to_point_list ( img0 ) 857 include_mask = [ ( p.x , p.y ) for p in points ] 858 img0.exclude_masks = [] 859 img0.include_masks = [ include_mask ] 860 861 # TODO: using .000001 instead of plain 0.0 is a workaround 862 # as the current reverse transform fails for 0.0 863 864 img0.overlap_center = point ( 0.000001 , 0.000001 ) 865 866 # TODO: maybe check here if rotation wouldn't be better 867 868 img0.rotate_overlap_pano = False 869 return True 870 871 # now we inspect the outline for 'crossings', where one point 872 # is 'inside' and the next isn't - or the other way round. 873 874 crossings = [] 875 amended_outline = [] 876 previous = outline[-1] 877 hop_threshold = 0.7 * min ( w0 , h0 ) # estimate only 878 879 for current in outline : 880 if current.inside != previous.inside : 881 # we have found a crossing, localize it precisely 882 xp = localize_crossing ( segment ( previous , current ) ) 883 amended_outline.append ( xp ) 884 crossings.append ( xp ) 885 elif current.inside : # both points are inside 886 # check for a crossing to the other side of a very-wide-angle image 887 if current.distance ( previous ) >= hop_threshold : 888 raise WoaError ( 'hop threshold exceeded, jump to other image margin' ) 889 # TODO: this isn't dealt with yet, we'd have to create two 890 # crossover points, one for going out and one for coming 891 # back in again. 892 amended_outline.append ( current ) 893 previous = current 894 895 if not crossings : 896 if not outline[0].inside : 897 vpr ( 'other image totally outside this one' ) 898 # ... but it might be totally enclosing it, in which case 899 # the symmetric call will deal with it. 900 return False 901 vpr ( 'other image totally within this one' ) 902 # note that the special case of total coincidence is dealt with 903 # separately above. 904 # we need to construct a mask with a hole. 905 # this is topologically simple, but to construct a mask that 906 # will work for hugin takes some trickery. I would have moved 907 # this code into a separate routine, but I need the infrastructure 908 # here, like the transforms and extents. 909 910 # first, find the leftmost point and extents of 911 # the 'other' image's outline. 912 913 minx = 1000000000.0 914 miny = 1000000000.0 915 maxx = -1000000000.0 916 maxy = -1000000000.0 917 918 for p in outline : 919 if p.twin.x < minx : 920 leftmost = p 921 minx = p.twin.x 922 elif p.twin.x > maxx : 923 maxx = p.twin.x 924 if p.twin.y < miny : 925 miny = p.twin.y 926 elif p.twin.y > maxy : 927 maxy = p.twin.y 928 929 # we make a list of outline points with 'leftmost' at it's tip 930 931 source = iter ( outline ) 932 head = [] 933 for p in source : 934 if p is leftmost : 935 inner_circuit = [ p ] 936 break 937 head.append ( p ) 938 inner_circuit.extend ( source ) 939 inner_circuit.extend ( head ) 940 941 # this happens to be the include mask, let's store it 942 include_mask = [ ( p.twin.x , p.twin.y ) 943 for p in inner_circuit ] 944 945 # we add start point again at the end to use this bit 946 # for the exclude mask 947 948 inner_circuit.append ( leftmost ) 949 950 # for the mask, we need these points, in reverse order, 951 # transformed to img0 coordinates, as tuples: 952 953 inner_circuit_transformed = [ ( p.twin.x , p.twin.y ) 954 for p in reversed ( inner_circuit ) ] 955 956 points , left , top , right , bottom = img_to_point_list ( img0 ) 957 958 # the 'mask with hole' starts with img0's corner points in 959 # clockwise order, plus 'leftmost' projected to the left 960 # margin 961 962 exclude_mask = [ ( left , top ) , 963 ( right , top ) , 964 ( right , bottom ) , 965 ( left , bottom ) , 966 ( left , leftmost.twin.y ) ] 967 968 969 # now the 'inner circuit' is appended (counterclockwise) 970 # and finally the inner circuit is linked again with 971 # the outline of img0 972 973 exclude_mask.extend ( inner_circuit_transformed ) 974 exclude_mask.append ( ( left , leftmost.twin.y ) ) 975 976 img0.exclude_masks = [ exclude_mask ] 977 img0.include_masks = [ include_mask ] 978 979 # TODO: using .000001 instead of plain 0.0 is a workaround 980 # as the current reverse transform fails for 0.0 981 982 img1.overlap_center = point ( 0.000001 , 0.000001 ) 983 984 img0.overlap_center = project_in ( img1.overlap_center ) 985 986 if ( maxx - minx ) < ( maxy - miny ) : 987 img0.rotate_overlap_pano = True 988 else : 989 img0.rotate_overlap_pano = False 990 991 # finally, we have to add img1's outline to img1 as it's include 992 # mask - it doesn't have an exclude mask. Note it is only in this 993 # special case where the symmetric case won't produce any masks, 994 # since it does rely on it being done here. 995 996 img1.include_masks = [ [ ( p.x , p.y ) 997 for p in outline ] ] 998 img1.exclude_masks = [] 999 1000 # now we're done with the 'mask with hole' and return True 1001 1002 return True 1003 1004 # the special case is dealt with, we continue with the ordinary: 1005 1006 # section removed, 'center' is calculated further down 1007 ## later on we need the 'center' of the overlap. We make an 1008 ## informed guess here by taking the mean of all crossover 1009 ## points. Notice that the mean is calculated in real 3D 1010 ## coordinates and reprojected onto the sphere. 1011 1012 #xm = 0.0 1013 #ym = 0.0 1014 #zm = 0.0 1015 #for xp in crossings : 1016 #xxp = transform_point ( xp , rtf1 ) 1017 #cart = polar2_to_cartesian3 ( xxp ) 1018 #xm += cart[0] 1019 #ym += cart[1] 1020 #zm += cart[2] 1021 #xm /= len ( crossings ) 1022 #ym /= len ( crossings ) 1023 #zm /= len ( crossings ) 1024 #cart = normalize ( ( xm, ym, zm ) ) 1025 #pol = cartesian3_to_polar2 ( cart ) 1026 #center = point ( pol[0] , pol[1] ) 1027 #center = transform_point ( center , tf0 ) 1028 1029 #print ( '******* old center calculation: %s' % center ) 1030 1031 outline = in_point_to_tip ( amended_outline ) 1032 1033 # now we generate separate lists for the parts of the 'other' outline 1034 # that lie inside this image. This is where we change to points 1035 # in 'this' image's coordinates. 1036 1037 chains = [] 1038 record = False 1039 for p in outline : 1040 if p.direction == 'in' : 1041 previous = None 1042 record = True 1043 if record : 1044 pin = p.twin 1045 pin.direction = p.direction 1046 if previous : 1047 previous.next = pin 1048 else : 1049 chains.append ( pin ) 1050 previous = pin 1051 if p.direction == 'out' : 1052 record = False 1053 1054 for xp in crossings : 1055 success = p0.take_in ( xp.twin ) 1056 if not success : 1057 raise WoaError ( 'failed to insert crossover point %s' % 1058 xp.twin ) 1059 1060 outline = [] 1061 for s in p0.segments : 1062 outline.append ( s.a ) 1063 1064 # of this contour, we use both parts - those that are 1065 # within the other image and those that aren't 1066 1067 inchains = [] 1068 outchains = [] 1069 previous = outline[-1] 1070 for p in outline : 1071 # print ( 'outline:' , p ) 1072 # note that when the OTHER image's contour comes IN 1073 # THIS image's contour goes OUT since we've traced 1074 # the contours both clockwise. Here we can append the 1075 # contour points without transformation. We make a doubly 1076 # linked list, since we need to go both ways 1077 if p.direction == 'in' : 1078 inchains.append ( p ) 1079 elif p.direction == 'out' : 1080 outchains.append ( p ) 1081 previous.right = p 1082 p.left = previous 1083 previous = p 1084 1085 # finally, we have to connect the out nodes 1086 # we walk the paths in different ways, depending on purpose. 1087 # first we generate the exclude mask. This is done by turning 1088 # 'left' at out nodes. 1089 1090 for p in inchains : 1091 # print ( 'inchain:' , p ) 1092 p.used = False 1093 1094 exclude_masks = [] 1095 1096 for p in inchains : 1097 if p.used : 1098 continue 1099 start = p 1100 mask = [] 1101 exclude_masks.append ( mask ) 1102 while True : 1103 # print ( p ) 1104 assert ( p.x , p.y ) not in mask 1105 mask.append ( ( p.x , p.y ) ) 1106 if p.direction == 'out' : 1107 keep = 'left' 1108 elif p.direction == 'in' : 1109 p.used = True 1110 keep = 'next' 1111 if keep == 'left' : 1112 p = p.left 1113 else : 1114 p = p.next 1115 if p == start : 1116 break 1117 1118 # next the include masks 1119 1120 include_masks = [] 1121 1122 for p in inchains : 1123 p.used = False 1124 1125 for p in inchains : 1126 if p.used : 1127 continue 1128 start = p 1129 mask = [] 1130 include_masks.append ( mask ) 1131 1132 while True : 1133 mask.append ( ( p.x , p.y ) ) 1134 if p.direction == 'out' : 1135 keep = 'right' 1136 elif p.direction == 'in' : 1137 p.used = True 1138 keep = 'next' 1139 if keep == 'right' : 1140 p = p.right 1141 else : 1142 p = p.next 1143 if p == start : 1144 break 1145 1146 # we use the include masks to find the maximal extents 1147 1148 minx = 1000000000.0 1149 maxx = -1000000000.0 1150 miny = 1000000000.0 1151 maxy = -1000000000.0 1152 1153 for m in include_masks : 1154 for p in m : 1155 minx = min ( p[0] , minx ) 1156 miny = min ( p[1] , miny ) 1157 maxx = max ( p[0] , maxx ) # KFJ fixed typo min() -> max() 1158 maxy = max ( p[1] , maxy ) # KFJ fixed typo min() -> max() 1159 1160 dx = maxx - minx 1161 dy = maxy - miny 1162 1163 if dx < dy : 1164 rotate_pano = True 1165 else : 1166 rotate_pano = False 1167 1168 # finally we attach all the mask data to 'this' image as additional 1169 # attributes to conveniently pass them around instead of having to 1170 # pass around numerous return values. 1171 1172 img0.exclude_masks = exclude_masks 1173 img0.include_masks = include_masks 1174 1175 # KFJ 2012-03-18 center is now calculated here: 1176 1177 center = calculate_overlap_center ( img0 ) 1178 1179 img0.overlap_center = center 1180 img0.rotate_overlap_pano = rotate_pano 1181 1182 # we return True to signal that masks have been created 1183 return True 1184 1185# attach_exclude_masks will aply the previously calculated masks to exclude 1186# the non-overlapping parts from the image in question. 1187 1188def attach_exclude_masks ( img ) : 1189 1190 # the mask data are calculated with the origin in the center of the 1191 # image, we have to change that to corner-relative pixel coordinates. 1192 1193 wh = img.getWidth() / 2.0 1194 hh = img.getHeight() / 2.0 1195 1196 for m in img.exclude_masks : 1197 Mask = hsi.MaskPolygon() 1198 for p in m : 1199 # KFJ 2012-03-17 added '- 0.5', because the coordinate system 1200 # used by PT for pixel coordinates, where the coordinates run 1201 # from -0.5 to w-0.5 and -0.5 to h-0.5 1202 cpx = p[0] + wh - 0.5 1203 cpy = p[1] + hh - 0.5 1204 # this fix I had to use earlier (to work around the assumed 1205 # mask bug, which isn't there after all) can now be taken out: 1206 #if abs ( cpx ) < POSITION_DELTA : 1207 #cpx = mask_offset 1208 Mask.addPoint ( hsi.FDiff2D ( cpx , cpy ) ) 1209 img.addMask ( Mask ) 1210 1211# to determine if it's worth our while to consider the overlap between 1212# two images for CP detection, we need the size of the overlapping area. 1213# This can be easily computed from the 'include' masks: 1214 1215def calculate_overlap_ratio ( img ) : 1216 1217 area2 = 0.0 1218 for m in img.include_masks : 1219 previous = m[-1] 1220 for current in m : 1221 dy = previous[1] - current[1] 1222 xsum = current[0] + previous[0] 1223 area2 += dy * xsum 1224 previous = current 1225 1226 overlap_area = area2 / 2.0 1227 vpr ( 'overlap area:' , overlap_area , 'pixels' ) 1228 w = img.getWidth() 1229 h = img.getHeight() 1230 img.overlap_ratio = overlap_area / ( w * h ) 1231 vpr ( 'overlap ratio:' , img.overlap_ratio ) 1232 1233# a routine to calculate the overlap's center from the include masks. 1234# the 'center' is found by calcuating the center of gravity off all 1235# mask's outlines, for simplicity's sake. 1236 1237def calculate_overlap_center ( img ) : 1238 1239 xsum = 0.0 1240 ysum = 0.0 1241 lsum = 0.0 1242 for m in img.include_masks : 1243 previous = m[-1] 1244 for current in m : 1245 # length of the edge 1246 dx = previous[0] - current[0] 1247 dy = previous[1] - current[1] 1248 l = math.sqrt ( dx * dx + dy * dy ) 1249 # center of the edge 1250 mx = ( previous[0] + current[0] ) / 2.0 1251 my = ( previous[1] + current[1] ) / 2.0 1252 # weighted sum of edge centers, and also total of 1253 # lengths, to calcualte the final outcome 1254 lsum += l 1255 xsum += l * mx 1256 ysum += l * my 1257 previous = current 1258 1259 return point ( xsum / lsum, ysum / lsum ) 1260 1261# the next routine will calculate the overlap of the images a and b 1262# in the panorama pano and generate control points for the overlapping 1263# area using the warped overlap method. 1264# overlap_threshold defines from what size overlap the images will be 1265# processed at all. The value refers to the normalized overlap area, 1266# that's the ration of pixels in overlapping areas to total pixels. 1267# scaling_factor defines by how much the images are scaled compared to 1268# the 'optimal size' that hugin would calculate for the panorama that 1269# generates the warped images. 1270 1271def img_from_file ( filename ) : 1272 1273 i = hsi.SrcPanoImage() 1274 i.setFilename ( filename ) 1275 if i.readEXIF(): 1276 return i 1277 return None 1278 1279def cps_from_overlap ( pano , a , b ) : 1280 1281 # hugin doesn't like wrongly-ordered numbers... 1282 1283 if a > b : 1284 help = a 1285 a = b 1286 b = help 1287 1288 # we make a subset panorama with only these two images and 1289 # set it to 360 degrees X 180 degrees equirect with width = 360 1290 # and height = 180, so that we can use panorama coordinates 1291 # directly as lat/lon values 1292 1293 subset = hsi.UIntSet ( [ a , b ] ) 1294 subpano = pano.getSubset ( subset ) 1295 1296 # if the images aren't active, nona will not process them. 1297 1298 subpano.activateImage ( 0 ) 1299 subpano.activateImage ( 1 ) 1300 1301 # what were image a and b in the original pano are now 0 and 1 1302 1303 img0 = subpano.getImage ( 0 ) 1304 img0.exclude_masks = [] 1305 img0.include_masks = [] 1306 img0.overlap_center = None 1307 img0.rotate_overlap_pano = False 1308 1309 img1 = subpano.getImage ( 1 ) 1310 img1.exclude_masks = [] 1311 img1.include_masks = [] 1312 img1.overlap_center = None 1313 img1.rotate_overlap_pano = False 1314 1315 pano_options = subpano.getOptions() 1316 pano_options.setProjection ( hsi.PanoramaOptions.EQUIRECTANGULAR ) 1317 pano_options.setHFOV ( 360.0 ) 1318 pano_options.setWidth ( 360 ) 1319 pano_options.setHeight ( 180 ) 1320 1321 subpano.setOptions ( pano_options ) 1322 1323 # now we call the routine to set exclude masks for the parts 1324 # of the images which don't overlap. Notice that this process 1325 # works only on the y, p, and r values found in the pto 1326 # we're working on, it may be totally wrong if these values 1327 # aren't roughly correct. With x degrees hov and 1/3 overlap, 1328 # an error of, say, x/20 degrees won't make to much difference. 1329 1330 mask_nonoverlaps ( subpano , img0 , img1 ) 1331 calculate_overlap_ratio ( img0 ) 1332 mask_nonoverlaps ( subpano , img1 , img0 ) 1333 calculate_overlap_ratio ( img1 ) 1334 1335 if ( img0.overlap_ratio <= args.threshold 1336 and img1.overlap_ratio <= args.threshold ) : 1337 vpr ( 'overlap is below or equal specified threshold' ) 1338 return 0 1339 1340 if ( img0.overlap_ratio > args.ceiling 1341 and img1.overlap_ratio > args.ceiling ) : 1342 vpr ( 'overlap is above specified ceiling' ) 1343 return 0 1344 1345 attach_exclude_masks ( img0 ) 1346 attach_exclude_masks ( img1 ) 1347 1348 # if we're here, the images have qualified for being submitted 1349 # to the warped overlap method. 1350 # uncomment the next three lines to write this intermediate pto 1351 # to disk for debugging purposes: 1352 1353 # ofs = hsi.ofstream ( 'subset0.pto' ) 1354 # subpano.writeData ( ofs ) 1355 # del ofs 1356 1357 # we transform the panorama so that: 1358 # - the assumed center of the overlap is in the center of the pano 1359 # - the larger extent of the overlap is parallel to the pano's horizon 1360 # these transformations assure that the distortions of the warped 1361 # images we create from the overlap are minimal. 1362 1363 # we need the reverse transform, from image to pano, for img0, 1364 # so we can shift the center: 1365 1366 rtf0 = hsi.Transform() 1367 rtf0.createInvTransform ( img0 , pano_options ) 1368 1369 # we transform our center point to pano coordinates, which amount 1370 # to lon/lat since we adjusted the panorama dimensions to 360X180 1371 1372 center = transform_point ( img0.overlap_center , rtf0 ) 1373 1374 # now we can use center's coordinates directly to specify the 1375 # needed rotations to put center at the pano's origin: 1376 1377 rp = hsi.RotatePanorama ( subpano , -center.x , 0.0 , 0.0 ) 1378 rp.runAlgorithm() 1379 rp = hsi.RotatePanorama ( subpano , 0.0 , center.y , 0.0 ) 1380 rp.runAlgorithm() 1381 1382 # if the masking routine has determined that we'd be better off 1383 # with the pano rotated 90 degrees, we do that as well 1384 1385 if img0.rotate_overlap_pano : 1386 rp = hsi.RotatePanorama ( subpano , 0.0 , 0.0 , 90.0 ) 1387 rp.runAlgorithm() 1388 1389 # we get the transforms from panospace to image coordinates 1390 # for this panorama, since later on our control points from 1391 # the warped imaged will have to be mapped back to the original 1392 # images with this very transform. We also keep a record of 1393 # the images' dimensions, expressed as an offset, since later 1394 # on we have to convert from image coordinates, which are 1395 # with origin in the image center, to pixel coordinates used 1396 # by the control points which put the origin in a corner. 1397 1398 pano_options = subpano.getOptions() 1399 img0 = subpano.getImage ( 0 ) 1400 img1 = subpano.getImage ( 1 ) 1401 1402 w0 = img0.getWidth() 1403 h0 = img0.getHeight() 1404 xcenter0 = w0 / 2.0 1405 ycenter0 = h0 / 2.0 1406 1407 w1 = img1.getWidth() 1408 h1 = img1.getHeight() 1409 xcenter1 = w1 / 2.0 1410 ycenter1 = h1 / 2.0 1411 1412 tf0 = hsi.Transform() 1413 tf0.createTransform ( img0 , pano_options ) 1414 1415 tf1 = hsi.Transform() 1416 tf1.createTransform ( img1 , pano_options ) 1417 1418 # subpano is the panorama we pass to nona for generation of 1419 # the warped overlap images. We need to increase the size to 1420 # near-original (like the 'optimal size' button in hugin) 1421 # here we may choose a different scale to further influence 1422 # the outcome of the CPG, depending on purpose: 1423 # to get as many CPs as possible, like for calibration, using 1424 # the original size or even a larger size should be optimal, 1425 # whereas if it's only a matter of gaining well-distributed CPs, 1426 # the size can be lower. 1427 1428 scale_calc = hsi.CalculateOptimalScale ( subpano ) 1429 scale_calc.run() 1430 optimal_width = scale_calc.getResultOptimalWidth() 1431 1432 # we make sure we use an even number for the width - and it has to 1433 # be integral as well: 1434 1435 subpano_width = 2 * int ( optimal_width * ( args.scale / 2.0 ) ) 1436 1437 vpr ( 'setting subpano width to' , subpano_width ) 1438 1439 # we set the width of the panorama 1440 1441 pano_options.setWidth ( subpano_width ) 1442 subpano.setOptions ( pano_options ) 1443 1444 # and write it to disk so we can pass it to nona 1445 1446 if args.prolific : 1447 warped_image_basename = '%s_%03d_%03d' % ( args.basename , a , b ) 1448 nona_file_name = warped_image_basename + '.pto' 1449 warped_pto = warped_image_basename + '_warped.pto' 1450 else : 1451 warped_image_basename = args.basename 1452 nona_file_name = '%s_base.pto' % args.basename 1453 warped_pto = '%s_warped.pto' % args.basename 1454 1455 ofs = hsi.ofstream ( nona_file_name ) 1456 subpano.writeData ( ofs ) 1457 del ofs 1458 1459 # we let nona execute this panorama to yield the warped overlaps. 1460 # the warped images are saved under the name specified 1461 # by warped_image_basename. 1462 # Notice that we assume nona is in the system path. 1463 1464 vpr ( 'creating warped overlap images with' , nona_file_name ) 1465 1466 command = [ 'nona' , 1467 '-o' , 1468 warped_image_basename , 1469 nona_file_name ] 1470 1471 # KFJ 2012-03-17 using run_helper_program() instead of subprocess.call() 1472 1473 # retval = subprocess.call ( command ) 1474 1475 retval , output = run_helper_program ( command ) 1476 1477 if retval or output: 1478 print ( 'nona returned %d, output:' % retval ) 1479 for line in output : 1480 print ( line , end = "" ) 1481 1482 # we check if nona really made the two images to process 1483 # and only proceed then - there may be corner cases when an 1484 # overlap was detected for (a,b) but not (b,a) - in this case 1485 # only one image might be made by nona and assuming there are two 1486 # (with the intended names) will trip the CPG which won't find 1487 # the second image 1488 1489 nona_output = [ '%s0000.tif' % warped_image_basename , 1490 '%s0001.tif' % warped_image_basename ] 1491 1492 # KFJ 2012-03-14 now cleaning up all transient files per default 1493 1494 cleanup_list = [] 1495 1496 for f in nona_output : 1497 if not os.path.exists ( f ) : 1498 # AFAICT this happens if the image boundaries do overlap 1499 # but there is no output because the overlap is masked 1500 vpr ( 'no intermediate file %s, no CPs added' % f ) 1501 return 0 1502 cleanup_list.append ( f ) 1503 1504 # next generate CPs from the warped images. What I do currently is 1505 # use autopano-sift-c with hardcoded parameters. apsc isn't easy 1506 # to handle for the purpose, so it has to be tricked into doing 1507 # the job. First, even though the images are really 360X180 equirect 1508 # in cropped TIFF format, as made by nona, apsc doesn't expect 1509 # such input and fails spectacularly if we pass 360 hfov and 1510 # RECTILINEAR as input. I think what happens is that it makes funny 1511 # assumptions and transforms the images somehow, resulting in 1512 # zero CPs. 1513 # so we trick autopano-sift-c to take the images, which are really 1514 # 360x180 equirect, in cropped tiff format, as 30 degree hfov 1515 # rectilinear so that it doesn't warp them itself. 1516 1517 # KFJ 2011-06-02 Finally I'm happy with a set of settings 1518 # for cpfind. The sieve2size of 10 is quite generous and I 1519 # suppose this way up to 250 CPs per pair can be found; if 1520 # even more are needed this is the parameter to tweak. 1521 # --fullscale now works fine and is always used, since 1522 # the scaling is done as desired when the warped images 1523 # are created. 1524 1525 # KFJ 2012-01-12 I found sieve2size of 10 is too generous, 1526 # so I set it down to 5, which is now my favourite setting. 1527 1528 # KFJ 2012-03-21 added --ransacmode hom. If warped images 1529 # are of different size due to masking in the original 1530 # pano, default ransac mode rpy fails, as the arbitrary 1531 # fov of 50 is applied to both warped images. rpy mode 1532 # needs precise fov data. 1533 1534 mute_command_output = True 1535 1536 if args.cpg == 'cpfind' : 1537 1538 cpfind_input = '_%s' % warped_pto 1539 cleanup_list.append ( cpfind_input ) 1540 1541 ofs = hsi.ofstream ( cpfind_input ) 1542 warp = hsi.Panorama() 1543 wi0 = img_from_file ( nona_output[0] ) 1544 wi1 = img_from_file ( nona_output[1] ) 1545 warp.addImage ( wi0 ) 1546 warp.addImage ( wi1 ) 1547 warp.writeData ( ofs ) 1548 del ofs 1549 1550 command = [ 'cpfind' , 1551 '--fullscale' , 1552 '--sieve2size' , '5' , 1553 '--ransacmode' , 'hom' , 1554 '-o' , warped_pto , 1555 cpfind_input ] 1556 1557 elif args.cpg == 'autopano-sift-c' : 1558 1559 # apsc works just fine like this: 1560 1561 vpr ( 'making pto for CP detection: %s' % warped_pto ) 1562 1563 command = [ 'autopano-sift-c' , 1564 '--projection' , '%d,30' % hsi.PanoramaOptions.RECTILINEAR , 1565 '--maxdim' , '%d' % subpano_width , 1566 '--maxmatches' , '0' , 1567 warped_pto ] + nona_output 1568 1569 # when using apsc, we let it display it's license warning etc: 1570 1571 mute_command_output = False 1572 1573 else : 1574 1575 raise WoaError ( 'only cpfind and autopano-sift-c are allowed as CPG' ) 1576 1577 vpr ( 'cp detection: %s %s' % ( command[0] , command[1:] ) ) 1578 1579 # KFJ 2012-03-17 using run_helper_program() instead of subprocess.call() 1580 1581 # retval = subprocess.call ( command ) 1582 1583 retval , output = run_helper_program ( command , 1584 mute_command_output ) 1585 1586 ## now the CPG has made a panorama for us, and the only thing we 1587 # care for from that is the CPs it's generated. Let's get them: 1588 1589 # KFJ 2012-03-18 check if the CPG actually produced output; 1590 # if not, assume there are no CPs. 1591 1592 if not os.path.exists ( warped_pto ) : 1593 cpv = hsi.CPVector() 1594 else : 1595 ifs = hsi.ifstream ( warped_pto ) 1596 warp = hsi.Panorama() 1597 warp.readData ( ifs ) 1598 del ifs 1599 1600 cleanup_list.append ( warped_pto ) 1601 1602 # we take the lot as a CPVector: 1603 1604 cpv = warp.getCtrlPoints() 1605 1606 vpr ( 'CPG found %d CPs' % len ( cpv ) ) 1607 1608 # We want to put the remapped CPs into this CPVector: 1609 1610 tcpv = hsi.CPVector() 1611 1612 # if we're prolific, we delete the CPs present in subpano 1613 # so only newly generated points will be seen 1614 1615 if args.prolific: 1616 subpano.setCtrlPoints ( tcpv ) 1617 1618 # the images have been saved by nona as cropped tiff. We need 1619 # to know the absolute positions of the CPs, so we have to 1620 # get the offsets. I do this using tiffdump, maybe there's 1621 # a better way? 1622 # Anyway, cropped TIFF is a fine thing, it saves processing 1623 # time and disk space, so we accept the little extra work: 1624 1625 xoff1 , yoff1 = get_tiff_offset ( nona_output[0] ) 1626 xoff2 , yoff2 = get_tiff_offset ( nona_output[1] ) 1627 1628 # we need to get back from absolute panorama coordinates to 1629 # 360X180 degrees, so we can map the CPs back to original image 1630 # coordinates. 1631 1632 scale = 360.0 / float ( subpano_width ) 1633 cps_added = 0 1634 1635 for cp in cpv : 1636 # the coordinates in cpv refer to the warped images. Since 1637 # these are, really, 360x180 equirect in cropped TIFF, all we 1638 # need to do with the coordinates the CPG has generated is 1639 # - add the TIFF crop offsets 1640 # KFJ 2012-03-17: 1641 # - add 0.5 te get from pixel coordinate to image-center-relative 1642 # - scale to 360X180 1643 # - shift the origin to the center by subtracting 180, or 90 1644 # vpr ( 'cp: ' , cp.x1 , cp.y1 , cp.x2 , cp.y2 ) 1645 x1 = scale * ( cp.x1 + xoff1 + 0.5 ) - 180.0 1646 x2 = scale * ( cp.x2 + xoff2 + 0.5 ) - 180.0 1647 y1 = scale * ( cp.y1 + yoff1 + 0.5 ) - 90.0 1648 y2 = scale * ( cp.y2 + yoff2 + 0.5 ) - 90.0 1649 # vpr ( 'angles' , x1 , y1 , x2 , y2 ) 1650 # - transform these angles into original image coordinates 1651 cp0 = transform_point ( point ( x1 , y1 ) , tf0 ) 1652 cp1 = transform_point ( point ( x2 , y2 ) , tf1 ) 1653 # vpr ( 'in img0:' , cp0 , 'in img1' , cp1 ) 1654 # - convert from image-center-relative coordinates 1655 # to corner origin pixel coordinates 1656 # KFJ 2012-03-17, added '- 0.5': 1657 x1 = cp0.x + xcenter0 - 0.5 1658 y1 = cp0.y + ycenter0 - 0.5 1659 x2 = cp1.x + xcenter1 - 0.5 1660 y2 = cp1.y + ycenter1 - 0.5 1661 1662 # KFJ 2011-06-13 the inverse transformation sometimes failed 1663 # on me producing nan (for example when reverse-transforming 1664 # (0,0)) - so I test here to avoid propagating the nan values 1665 # into the panorama. 1666 # I've changed the places where (0,0) would be fed into the 1667 # transform, so I avoid the problem, but maybe it can still 1668 # crop up somehow. 1669 # TODO: remove test if certain that outcome is never nan. 1670 1671 if ( math.isnan(x1) 1672 or math.isnan(x2) 1673 or math.isnan(y1) 1674 or math.isnan(y2) ) : 1675 1676 # Instead of raising an exception I just ignore 1677 # the problem for now, skip the point and emit a message. 1678 1679 # raise WoaError ( 'coordinate turned out nan' ) 1680 1681 print ( 'failed to convert pano coordinates %s\n' % cp ) 1682 print ( 'back to image coordinates. point is ignored' ) 1683 continue 1684 1685 # - generate a ControlPoint object if the point is inside 1686 # the image - sometimes the CPGs produce CPs outside the 1687 # overlap area, I haven't found out why 1688 # KFJ 2012-03-18 I think I know now, but I'm not sure 1689 # TODO later remove the print statement 1690 if ( x1 < -0.5 or x1 > w0 - 0.5 1691 or y1 < -0.5 or y1 > h0 - 0.5 1692 or x2 < -0.5 or x2 > w1 - 0.5 1693 or y2 < -0.5 or y2 > h1 - 0.5 ) : 1694 vpr ( '****** CP outside image: (%f,%f) (%f,%f)' 1695 % ( x1, y1 , x2, y2 ) ) 1696 continue 1697 1698 tcp = hsi.ControlPoint ( a , x1 , y1 , b , x2 , y2 , 1699 hsi.ControlPoint.X_Y ) 1700 cps_added += 1 1701 # vpr ( 'as CP:' , tcp.x1 , tcp.y1 , tcp.x2 , tcp.y2 ) 1702 # - and put that into the panorama 1703 # tcpv.append ( tcp ) 1704 pano.addCtrlPoint ( tcp ) 1705 if args.prolific: 1706 subpano.addCtrlPoint ( tcp ) 1707 1708 # if we're prolific, we write the nona file with added CPs 1709 # so the yield of the CPG can be inspected 1710 1711 if args.prolific : 1712 ofs = hsi.ofstream ( nona_file_name ) 1713 subpano.writeData ( ofs ) 1714 del ofs 1715 else : 1716 for f in cleanup_list : 1717 # vpr ( 'removing temporary file %s' % f ) 1718 os.remove ( f ) 1719 1720 # finally we return the number of control points we've added 1721 # to the panorama 1722 1723 return cps_added 1724 1725# from now on it's basically administration: determine which image pairs 1726# are looked at and what the parameters are. 1727 1728def process_image_set ( pano , image_set ) : 1729 1730 vpr ( 'processing image set' , image_set ) 1731 1732 # the user may have specified a 'focus' image, in which case 1733 # only combinations of the focus image with the other images 1734 # are looked at. Otherwise, all combinations 1735 1736 if args.focus : 1737 pairs = [ ( args.focus , image ) 1738 for image in image_set 1739 if args.focus != image ] 1740 else : 1741 pairs = list ( itertools.combinations ( image_set , 2 ) ) 1742 1743 # finally we use the workhorse routine on all pairs: 1744 1745 total_cps_added = 0 1746 1747 for pair in pairs : 1748 1749 # first check if a stop file is set and present 1750 1751 if args.stop and os.path.exists ( args.stop ) : 1752 1753 print ( 'found stop file %s' % args.stop ) 1754 print ( 'ending woa, %d new CPs' % total_cps_added ) 1755 break 1756 1757 vpr ( 'examining image pair' , pair ) 1758 try: 1759 cps_added = cps_from_overlap ( pano , *pair ) 1760 total_cps_added += cps_added 1761 if cps_added: 1762 print ( 'image pair %s: %3d control points added' 1763 % ( '(%3d, %3d)' % pair , cps_added ) ) 1764 except Exception: 1765 1766 # KFJ 2012-03-18 took over catch-all exception handling 1767 # from hpi.py so that now there's a stack trace printout 1768 1769 # first we retrieve the details of the exception: 1770 1771 exc_type, exc_value, exc_traceback = sys.exc_info() 1772 1773 # emit a message with the type and text of the exception: 1774 1775 print('matching pair %s failed with %s: "%s"' 1776 % ( pair, 1777 exc_type.__name__, 1778 exc_value) ) 1779 1780 # and, for good measure, a stack traceback to point 1781 # to the precise location and call history of the error: 1782 1783 print ( 'stack traceback:' ) 1784 traceback.print_tb ( exc_traceback ) 1785 1786 return total_cps_added 1787 1788# entry() is what is called directly when the hugin plugin interface calls 1789# this script. Since there is no parameter passing yet, the routine figures 1790# out which image pairs to process by looking at which images are 'active'. 1791 1792def entry ( pano ) : 1793 1794 # we emulate a parameter set as it would have been produced 1795 # by the argument parser in the CLI version: 1796 1797 global args 1798 args = argobj() 1799 1800 # see if all required helper programs are accessible 1801 1802 check_helpers() 1803 1804 vpr ( 'entry: parameters used for this run:' ) 1805 vpr ( 'ceiling:' , args.ceiling ) 1806 vpr ( 'cpg:' , args.cpg ) 1807 vpr ( 'focus: ' , args.focus ) 1808 vpr ( 'margin: ' , args.margin ) 1809 vpr ( 'scale: ' , args.scale ) 1810 vpr ( 'stop: ' , args.stop ) 1811 vpr ( 'thresh: ' , args.threshold ) 1812 vpr ( 'prolific' , args.prolific ) 1813 vpr ( 'verbose:' , args.verbose ) 1814 1815 # we restrict the search to the active images 1816 # in the panorama. This is working around the fact that 1817 # as of this writing plugins can't receive parameters, 1818 # so we use the 'active' state of the images for selection. 1819 1820 image_set = [ image for image in pano.getActiveImages() ] 1821 1822 # we doublecheck: if only one image is active, we change 1823 # our logic and take the only active image as our focus 1824 # and all images as image set: 1825 1826 if len ( image_set ) == 1 : 1827 args.focus = image_set[0] 1828 image_set = [ int ( image ) 1829 for image 1830 in range ( pano.getNrOfImages() ) ] 1831 1832 total_cps_added = process_image_set ( pano , image_set ) 1833 1834 # KFJ 2012-03-18 I used to return total_cps_added. 1835 # This resulted in the 'script returned XXX' popup 1836 # turining up, usually, hidden behind other windows 1837 # so now I return 0 and there is only a popup if an 1838 # error occurred. 1839 1840 # return total_cps_added 1841 return 0 1842 1843# main() is called if the program is called from the command line. 1844# In this case we have a proper set of arguments to process: 1845 1846def main() : 1847 1848 # when called from the command line, we import the argparse module 1849 1850 import argparse 1851 1852 # and create an argument parser 1853 1854 parser = argparse.ArgumentParser ( 1855 formatter_class=argparse.RawDescriptionHelpFormatter , 1856 description = gpl + ''' 1857 woa will look at every possible pair of images that can be made 1858 from the image numbers passed as -i parameters, or from all images 1859 in the pto if no -i parameters are passed. If the images in a pair 1860 overlap significantly, the overlapping parts will be warped and 1861 a CPG run on them. The resulting CPs will be transferred back into 1862 the original pto file. 1863 ''' ) 1864 1865 parser.add_argument('-c', '--ceiling', 1866 metavar='<overlap ceiling>', 1867 default = 1.0 , 1868 type=float, 1869 help='ignore overlaps above this value (0.0-1.0)') 1870 1871 # this next paramater does nothing per se - if image sets are 1872 # specified with -i or -x and no other parameter with a leading 1873 # '-' follows naturally, put in -e to finish the sequence 1874 # TODO: see if this can't be done more elegantly 1875 1876 parser.add_argument('-e' , '--end', 1877 help='dummy: end a group of image numbers') 1878 1879 parser.add_argument('-g', '--cpg', 1880 metavar='<cpfind or autopano-sift-c>', 1881 default = 'cpfind', # KFJ 2011-06-13, was apsc 1882 type=str, 1883 help='choose which CP generator to use') 1884 1885 parser.add_argument('-f' , '--focus' , 1886 metavar='<image number>', 1887 default=None, 1888 type=int, 1889 help='only look for overlaps with this image') 1890 1891 parser.add_argument('-i', '--images', 1892 metavar='<image number>', 1893 default = [] , 1894 nargs = '+', 1895 type=int, 1896 help='image numbers to process') 1897 1898 parser.add_argument('-m', '--margin', 1899 metavar='<margin>', 1900 default = 0 , 1901 type=int, 1902 help='widen examined area [in pixels]') 1903 1904 parser.add_argument('-o', '--output', 1905 metavar='<output file>', 1906 default = None, 1907 type=str, 1908 help='write output to a different file') 1909 1910 parser.add_argument('-b', '--basename', 1911 metavar='<base name>', 1912 default = 'woa', 1913 type=str, 1914 help='common prefix for intermediate files') 1915 1916 parser.add_argument('-p', '--prolific', 1917 action='store_true', 1918 help='keep all intermediate files') 1919 1920 parser.add_argument('-t', '--threshold', 1921 metavar='<overlap threshold>', 1922 default = 0.0 , 1923 type=float, 1924 help='ignore overlaps below this threshold (0.0-1.0)') 1925 1926 parser.add_argument('-s', '--scale', 1927 metavar='<scaling factor>', 1928 default = 0.25 , 1929 type=float, 1930 help='scaling for warped overlap images') 1931 1932 parser.add_argument('-v' , '--verbose', 1933 action='store_true', 1934 help='produce verbose output') 1935 1936 parser.add_argument('-x', '--exclude', # KFJ 2011-06-02 1937 metavar='<image number>', 1938 default = [] , 1939 nargs = '+', 1940 type=int, 1941 help='image numbers to exclude') 1942 1943 parser.add_argument('-z' , '--stop', 1944 metavar = '<stop file>' , 1945 type = str , 1946 default = None , 1947 help = 'stop woa if this file is present' ) 1948 1949 parser.add_argument('input' , 1950 metavar = '<pto file>' , 1951 type = str , 1952 help = 'pto file to be processed' ) 1953 1954# we may add logging at some point... 1955 1956## parser.add_argument('-l', '--log', 1957## metavar='<log file>', 1958## type=argparse.FileType('w'), 1959## help='write log file') 1960 1961 # if the argument count is less than two, there aren't any arguments. 1962 # Currently this is considered an error and triggers help output: 1963 1964 if len ( sys.argv ) < 2 : 1965 parser.print_help() 1966 return 1967 1968 # we parse the arguments into a global variable so we can access 1969 # them from everywhere without having to pass them around 1970 1971 global args 1972 args = parser.parse_args() 1973 1974 # see if all required helper programs are accessible 1975 1976 check_helpers() 1977 1978 vpr ( 'main: parameters used for this run:' ) 1979 vpr ( 'ceiling:' , args.ceiling ) 1980 vpr ( 'exclude:' , args.exclude ) 1981 vpr ( 'focus: ' , args.focus ) 1982 vpr ( 'images: ' , args.images ) 1983 vpr ( 'margin: ' , args.margin ) 1984 vpr ( 'output: ' , args.output ) 1985 vpr ( 'scale: ' , args.scale ) 1986 vpr ( 'stop: ' , args.stop ) 1987 vpr ( 'thresh: ' , args.threshold ) 1988 vpr ( 'prolific' , args.prolific ) 1989 vpr ( 'ptofile:' , args.input ) 1990 vpr ( 'verbose:' , args.verbose ) 1991 1992 # see if we can open the input file 1993 1994 ifs = hsi.ifstream ( args.input ) 1995 if not ifs.good() : 1996 raise WoaError ( 'cannot open input file %s' % args.input ) 1997 1998 pano = hsi.Panorama() 1999 success = pano.readData ( ifs ) 2000 del ifs 2001 if success != hsi.DocumentData.SUCCESSFUL : 2002 raise WoaError ( 'input file %s contains invalid data' % args.input ) 2003 2004 # let's see if the image numbers passed are correct if there are any 2005 2006 ni = pano.getNrOfImages() 2007 vpr ( 'found %d images in panorama' % ni ) 2008 if ni < 2 : 2009 raise WoaError ( 'input file %s contains less than two images' 2010 % args.input ) 2011 2012 for nr in range ( ni ) : 2013 img = pano.getImage ( nr ) 2014 name = img.getFilename() 2015 # check if we can access this image 2016 # TODO: doublecheck in folder of pto file 2017 if not os.path.exists ( name ) : 2018 raise WoaError ( "cannot access image %s (try cd to pto's directory)" % name ) 2019 vpr ( 'image %d: %s' % ( nr , name ) ) 2020 2021 for nr in args.images : 2022 if nr < 0 or nr >= ni : 2023 raise WoaError ( 'no image number %d in input file' % nr ) 2024 2025 if args.focus : 2026 if not ( 0 <= args.focus < ni ) : 2027 raise WoaError ( 'invalid focus image %d' % args.focus ) 2028 2029 if args.threshold > args.ceiling : 2030 raise WoaError ( 'threshold must be below ceiling' ) 2031 2032 # if a separate output file was chosen, we open it now to avoid 2033 # later disappointments 2034 2035 if args.output: 2036 ofs = hsi.ofstream ( args.output ) 2037 if not ofs.good() : 2038 raise WoaError ( 'cannot open output file %s' % args.output ) 2039 2040 if len ( args.images ) : 2041 2042 # if specific image numbers have been passed as arguments, 2043 # we look for CPs no matter if they're active or not 2044 2045 image_set = [ int ( image ) for image in args.images ] 2046 2047 else : 2048 2049 # otherwise, we look at all images. Note that this behaviour 2050 # is different to the plugin version which looks at the 'active' 2051 # state of the images. 2052 2053 image_set = [ int ( image ) 2054 for image 2055 in range ( pano.getNrOfImages() ) ] 2056 2057 # if exclusions are given, they are taken away from the set 2058 # of images to be processed now (KFJ 2011-06-02) 2059 2060 if args.exclude : 2061 2062 orig = set ( image_set ) 2063 sub = set ( args.exclude ) 2064 image_set = orig - sub 2065 2066 # all set up. We process the image set 2067 2068 total_cps_added = process_image_set ( pano , image_set ) 2069 2070 # if no different output file was specified, overwrite the input 2071 2072 if not args.output : 2073 ofs = hsi.ofstream ( args.input ) 2074 if not ofs.good() : 2075 raise WoaError ( 'cannot open file %s for output' % args.input ) 2076 2077 success = pano.writeData ( ofs ) 2078 del ofs 2079 if success != hsi.DocumentData.SUCCESSFUL : 2080 raise WoaError ( 'error writing pano to %s' % args.input ) 2081 2082 # done. 2083 2084 return total_cps_added 2085 2086# finally the test for the invocation mode: 2087# are we main? This means we've been called from the command line 2088# as a standalone program, so we have to go through the main() function. 2089 2090if __name__ == "__main__": 2091 2092 try: 2093 2094 total_cps_added = main() 2095 2096 except WoaError as e : 2097 2098 print ( 'Run Time Error: %s' % e ) 2099 sys.exit ( -1 ) 2100 2101 sys.exit ( total_cps_added ) 2102