Take me to pcigeomatics.com

PCI Geomatics Help Center

How can we help you today?

Iterative GCP Collection (Until Success)

PCI Geomatics -

Overview


 

Creating highly-automated orthorectification workflows is one of the many strengths of developing with Geomatica and python. In this example, we will demonstrate a concept to further improve the success of automatic GCP collection through iterative and intelligent parameter setting.

Data Package


 

  • Click Here to download the data required for this tutorial

Watch this Tutorial


 

<Embed Video Here>

Benefits


 

Geomatica’s python library includes very powerful orthorectification functions for automation, such as, automatic GCP collection, automatic bundle adjustment, adjust ortho and more. While these functions are fast and generally succeed with the default settings, some parameters can affect the success of the process with certain images. As such, intelligently iterating through key settings can further improve the success of orthorectification when working on a batch.

Prerequisites


 

  • Geomatica 2015 or later installed
  • Python 101 (basic understanding of python coding concepts)
    • try, except, finally statements
    • How to access modules and classes
    • How to call a function
    • If-conditions
    • For-loops and while-loops
  • Some experience programming with Geomatica (recommended)

Try It First!


 

Try running the script provided with the data package and see what happens. You only need to update the working_dir parameter

By monitoring the terminal, you should notice that the first image runs successfully in the first iteration, but the second image does not succeed until the 3rd iteration.

Tutorial


 

To demonstrate how to create an intelligent workflow to iteratively change automatic GCP collection settings based on the type of failure, we are going to perform GCP collection and orthorectify two Rapideye scenes with different positional errors. In the first iteration, the image with a smaller positional error will succeed, but the image with the larger positional error will fail. We will create a set of rules that will evaluate why the second image failed and then determine which parameters in the automatic GCP collection algorithm should be changed and how, to improve the chances of success in the next run.

1. Import necessary modules & setup inputs/outputs


 

The first step is to import the Geomatica modules we will be calling in this tutorial and setting up the variables that will point to our input file(s), working file(s) (temporary) and output file(s).

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
import os
import fnmatch
from pci.fimport import fimport
from pci.autogcp import autogcp
from pci.gcprefn import gcprefn
from pci.nspio import Report, enableDefaultReport
from pci.exceptions import PCIException

# set input and output files for workflow
working_dir = r'C:\downloads\iter_gcp_coll'  # set this path to the dir where you extracted the tutorial data
input_dir = os.path.join(working_dir, 'Raw')
file_filter = 'K3*AUX.xml'
output_dir = os.path.join(working_dir, 'output')
ref_image = os.path.join(working_dir, r'Reference\ref_LC80160282015300LGN00_ORTHO_PAN.pix')
dem = os.path.join(working_dir, r'Reference\dem_Ottawa_10m.pix')
rep_file = os.path.join(output_dir, 'gcp_report.txt')

# check that the output folder exists. If not, create it!
if not os.path.exists(output_dir):
    os.makedirs(output_dir)

 

In the above code block, lines 1 and 2 import standard python modules os and fnmatch. Lines 3 to 7 import the different Geomatica modules that we will be working with:

  • fimport – Geomatica’s algorithm for importing raw satellite into a .pix This is required for photogrammetric operations.

 

Note: fimport can be replaced with link algorithm, which is faster and does not copy the data.

 

  • autogcp – Geomatica’s algorithm for automatically finding ground control points (GCPs) based on image matching techniques
  • gcprefn – Geomatica’s algorithm for automatically detecting bad GCP matches and removing them from the model
  • Report, enableDefaultReport – are two functions/classes that allow dumping output notifications to a report file
  • PCIException – is a class for handling errors specific to Geomatica’s algorithms

Line 10 points to a variable called working_dir that contains the path to the working directory. You will need to change this to the root folder of the tutorial data you downloaded i.e. r’C:\...\iter_gcp_coll’. Lines 11 to 16 are variable to hold sub directories or output filenames. These variables should not be changed.

Lines 19 and 20 do a quick check to see if the output directory, held in the output_dir exists. If not, it will be created.

 

2. Search input directory for all valid images


 

In this step, we will create a simple generator that will search the input directory for all valid input files and prepare them for processing

23
24
25
26
27
28
29
30
# Generator to get all images in a folder and prepare them for batch processing
def get_batch(in_dir, img_filter):
    for r, d, f in os.walk(in_dir):
        for in_file in fnmatch.filter(f, img_filter):
            yield os.path.join(r, in_file)


raw_images = get_batch(input_dir, file_filter)

 

On line 24 we define a function called get_batch, which takes in the in_dir variable and img_filter variable. The in_dir variable points to the input directory that it should recursively check for input files and the img_filter variable is a string that is used to determine which files are valid. Lines 25 to 27 loop through the directory and sub directories (if exist) and search for all files that match the img_filter string, using fnmatch.

On line 27, the yield operator is used to return every item (valid input file) on-the-fly.

On line 30, we output each iteration from the get_batch generator and temporarily store it in the raw_images variable.

 

3. Loop through raw images and import them to .pix


 

In this step, we will create a for-loop that will loop through each input image one at a time and perform all processing (import àiterative GCP collection). However, at this point in the tutorial, we will only focus on the import part of the loop.

32
33
34
35
for raw_img in raw_images:
    base_name = os.path.splitext(os.path.basename(raw_img))[0]
    pix_img = os.path.join(output_dir, base_name + '_ingested.pix')
    fimport(fili=raw_img, filo=pix_img)

 

On line 32 we create the for-loop that all remaining code will reside in. In the for-loop we will get each raw_img item that is yielded from raw_images.

On line 33, we create a variable called base_name and set it equal to the filename of the input file. The os.path.splitext() function is used to retrieve the filename without the path or file extension. Line 34 defines a variable called pix_image, and uses the function os.path.join() to add the output path to the base_name and appends a suffix with the .pix file extension.

On line 35 we run Geomatica algorithm, fimport() and pass the following input keyword arguments:

  • fili=raw_imgraw_img (input image) variable is assigned to the fili keyword
  • filo=pix_imgpix_img (output image) variable is assigned to the filo keyword

 

4. Setup key GCP Collection parameters for first iteration


 

Here we are simply creating variables for key parameters used by the AUTOGCP and GCPREFN Geomatica algorithms. We will set them equal to the ideal settings we would like them to use when searching for GCPs. However, if the first iteration fails to find or accept sufficient number of GCPs, these settings will be changed by code we will create later in this tutorial. What parameters are altered and how they are altered will depend on what type of failure occurs; we will discuss this in more details later.

37
38
39
40
i = 0
    radius = [50]  # 50 pixel search radius to begin with
    gcp_refine = [5, 1, 1]
    corr_score = [0.8]

 

On line 37 of the above code block, we create a variable i and set it equal to 0, this variable will keep track of the number of iterations we try. By not managing this properly, we run the chance of creating an infinite loop.

Lines 38 to 40 are the initial settings for the first GCP collection iteration:

  • radius = [50] – This will define the search radius to 50 pixels
  • gcp_refine = [5,1,1] – This defines the refinement method and extra parameters. The first value (5), tells the GCPREFN algorithm to remove GCPs until a certain RMSE is achieved. The second and third values define the RMSE in the x and y direction, respectively
  • corr_score = [0.8] – this defines the minimum correlation score that must be achieved for a feature to be accepted as a candidate GCP

5. Create while-loop and setup output parameters


 

In this step we will create a while-loop which will ensure we do not accidently create an infinite loop. The while loop will break after the variable i has increased to a certain value, 3 in this case. We will also create variables that will be used to collect output information about the success of AUTOGCP and GCPREFN.

41
42
43
44
45
while i < 3:
        print '\nStarting iteration {0} on image: {1}'.format(i + 1, base_name)
        gcps_found = []
        gcps_accepted = [0]
        refined_gcp_seg = []

 

In the above code block, line 41 defines the while-loop, While i < 3: which simply says that once i = 3, break the loop. In other words, if the processing is not successful after 3 iterations, it will fail the image and move on to the next image.

Lines 43 to 45 are variables that will be used to hold the output parameters from AUTOGCP and GCPREFN:

  • gcps_found = [] – empty list that will be populated by the AUTOGCP algorithm with the number of GCPs it successfully found
  • gcps_accepted = [0] – list with a value of 0 in it. This list will be modified by the GCPREFN algorithm with the number of GCPs accepted. It is initially set with a value of 0 in case no GCPs are found and therefore GCPREFN is not run.
  • refined_gcp_seg = [] – empty list that will be populated by the GCPREFN algorithm with the segment number that contains the refined GCPs

 

6. Run AUTOGCP and GCPREFN algorithms


 

In this step we will setup the report file so that the output notifications from AUTOGCP and GCPREFN are dumped to a report file. We will then run the AUTOGCP and GCPREFN Geomatica algorithms.

 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
try:
            Report.clear()
            enableDefaultReport(rep_file)
            gcp_seg = autogcp(mfile=pix_img,
                              dbic=[1],
                              filref=ref_image,
                              filedem=dem,
                              searchr=radius,
                              numgcps=gcps_found  # output parameter that lets us know if sufficient GCPs found
                              )

            if gcps_found[0] != 0:
                gcprefn(file=pix_img,
                        dbgc=gcp_seg,
                        modeltyp='RF',
                        reject=gcp_refine,
                        imstat=gcps_accepted,  # output parameter that lets us know if sufficient GCPs were accepted
                        lasc=refined_gcp_seg
                        )
        except PCIException, e:
            print e
        finally:
            enableDefaultReport('term')  # this will close the report file

 

We start on line 47 with a try: statement, which tells python to first try running the code inside and if an error occurs, handle it afterwards. This is known as EAFP (Easier to ask forgiveness than permission).

On lines 48 and 49, we initialize the report by running Report.clear() and enableDefaultReport(rep_file). This means that all output notifications will be written to the file held in the rep_file variable until the report is set back to term.

On line 50 to 56 we call the AUTOGCP algorithm. We set it equal to gcp_seg, as the algorithm returns the gcp segment that holds the GCPs that were collected. Check out the AUTOGCP documentation to understand what each parameter does.

On line 58 creates an if-condition that checks to see if the number of GCPs found by AUTOGCP is not equal to 0. If AUTOGCP failed to collect GCPs than there is no point to running GCPREFN. If the condition on line 58 is met, than the GCPREFN algorithm is run (lines 59 to 65).

Lines 66 and 67 simply handle any PCI exception that may have occurred by printing the error message and moving on (program will not crash). Lines 68 and 69 close the report and set the reporting to ‘term’.

 

7. Evaluate GCP results of current iteration and determine next steps


 

In this last block of code, we will evaluate the results of the GCP collection from the current iteration. If it was successful, we will break the loop and move on to the next image. If it was not successful, we will investigate why by mainly considering whether not enough GCPs were found or if not enough GCPs were accepted. We will then change the key parameters in a manner that will improve our chances of finding more GCPs or finding better GCPs.

 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
if gcps_accepted[0] >= 10:
            print 'Sufficient GCPs accepted! Found {0} GCPs and accepted {1} GCPs'.format(gcps_found[0],
                                                                                          gcps_accepted[0])
            break
        else:
            if gcps_found[0] <= 10:
                print 'Insufficient GCPs found! Found {0} GCPs'.format(gcps_found[0])
                radius[0] += 10
                print 'Increasing search radius to {0}'.format(radius[0])

                if corr_score >= 0.75:
                    corr_score[0] -= 0.03
                    print 'Decreasing Min Correlation Threshold to {0}'.format(corr_score[0])
            else:
                print 'Insufficient GCPs accepted! Found {0} GCPs and accepted only {1} GCPs'.format(gcps_found[0],
                                                                                                     gcps_accepted[0])
                gcp_refine[1] += 0.25
                gcp_refine[2] += 0.25

                if corr_score <= 0.85:
                    corr_score[0] += 0.03

                print 'Increasing Min Correlation Threshold to {0}' \
                      ' Increasing Max RMSE threshold to {1}'.format(corr_score[0], gcp_refine[1])

        i += 1

 

Note: None of the print statements will be explained. They are just included so that you can clearly see what is happening when you run the script.

On line 71 we perform our first condition to check and see if 10 or more GCPs were accepted. If so, we consider this a success and break the loop (line 74) and move on to the next image.

On line 75 we use the else keyword to handle the condition when less than 10 GCPs are accepted, which is considered a failure.

On line 76 we create another if-statement nested inside the original one to check if the reason we failed was because not enough GCPs were initially found. If this condition is true, we will increase the search radius by 10 pixels on line 78 (because perhaps the search radius was too small). Line 81 will check to see if the minimum correlation score is greater than or equal to 0.75, if so, it will lower it by 0.03 (line 82).

On line 84 we use the else condition in the case that enough GCPs were found, but not enough GCPS were accepted. This means we need to either find better GCPs are lower our expectations of the accuracy of the GCPs. In this case we do a little of both. First, on lines 87 and 88 we slightly increase the minimum acceptable RMSE by 0.25 pixels and finally on line 90 and 91 we increase the minimum correlation score, only if it is not already too high (above 0.85). We increase the correlation score so that only GCPs with better correlations are retained in the hopes that they will have higher accuracy.

Lastly, on line 96 we increase the value of i by 1 i += 1, which is required so our while-loop does not have a chance to fail for all of eternity!

Have more questions? Submit a request

Comments

Powered by Zendesk