ChartGeoreferencer
Written by Rob Truxler Friday, 09 April 2010 22:17
For centuries cartographers have been charting coast lines to ensure the expeditious and safe passage of ships at sea. Today, military, merchant marines, and recreational boaters use nautical charts for the same reasons. Additionally, these charts are used today by civil engineers, by scientists, and by emergency responders.
Nautical charts change rapidly. As the periodic push and pull of the ocean's tides erode our coastlines, depth changes and new hazards arise that require representation on these charts. As the magnetic north pole's deviation from Earth's true north pole changes annually, the mariner must know the local magnetic declination in order to use his compass. As bridges and other structures change, so might the commercial shipping lanes that pass underneath. For these reasons, much effort is invested in keeping nautical charts accurate and up to date.
Analyzing the changes between charts may be valuable to historians, climatologists, and ecologists. There is no guarantee, however, that any two nautical charts that cover the same area share the same boundaries. In fact, the National Oceanic and Atmospheric Association (NOAA) maintains and distributes the 2,229 chart panels that vary in scale, projection, and purpose. One chart might be a high-level informational chart that covers the entire Great Lakes region, another chart may cover all of Lake Michigan, one chart for all of Chicago, and maybe five charts for detailed harbor charts within the city. To analyze the spatial and temporal changes between a pair of charts, they first have to be projected into the same viewing space. That is, both charts have to be in the same projection, the same scale, and overlapping so that similar features would line up.
In order to reproject nautical charts from their rasterized primitive form to a common projection in which all charts can be viewed, there first must be knowledge about the latitude and longitude of a number of points in the image. For modern charts this data is available if the charts are distributed as GeoTIFFs,or some other format that supports geographic meta-data. For historical charts, some dating back to 1749, geo-reference data is simply unavailable.
In order to generate the pixel-to-latitude/longitude mappings for a given image, one would have to find the lines of latitude and longitude in the nautical chart, find their points of intersection, and label their coordinates. Finding the "lines" of latitude and longitude is no trivial task. A scanned nautical chart can be of such high resolution that even opening the it in a traditional image editing program would require more memory than most consumer computers have (>6GB).
The Solution?
Proposed is a user interface which simplifies the process of generating geographic meta-data for a nautical chart. The system would extract information from the image to make a guess at which locations in the image represent lines of latitude or longitude. Valid lines of latitude or longitude could then be selected by the user from the hints provided by the system. Once the user selects a single line, the system could then search through its set of possible lines, shrinking the search space smaller to only include lines that are parallel and on the same grid as the first line. Once a sufficient number of parallel lines is selected, the user can use the image data at the end points of each line segment to identify the actual latitude or longitude. Since the line segments in a nautical chart tend to terminate at the chart's margins and since each line usually has a label, a sub-image can be extracted at the line segments ends and the user can use that image data to properly label the line.
Once the above process is completed for more than two lines of latitude, the same process can be repeated for lines of longitude. Once two or more lines of longitude have been labeled by the user, geo-references can be computed at the intersection of each line of latitude and longitude. This data can then be used to project the original image onto a map.
A system was designed to support the user interface described above.
In order to compute the initial guesses at which lines represent lines of latitude and longitude, the system used the Hough transform (1)(2). This particular implementation of the Hough transform quantized values for R (a line's distance to the image origin) and theta (a line's angle relative to the x-axis). R values were quantized by using a bin size of 2 and theta values were quantized with a bin size of 5 degrees. These values were measured at each "edge" location in the image.
Pre-processing
In order to find edges, a number of pre-processing steps were applied to the source image.
Sub-sampling
Since the scanned images of nautical charts tend to be very large (some up to 18,000 pixels wide), a sub-sampled version was required in order to efficiently compute the Hough transform. Each chart was subsampled by a factor of 1000/(the chart's maximum dimension). If the original chart was 15000x20000 pixels, it was shrunk by a factor of (1/20) so that its sub-sampled version had dimensions 750x1000.
Pixel Intensity Thresholding
The sub-sampled image was then thresholded so that the result was a binary image that would hopefully include all of the graticule from latitude and longitude lines. The grayscale value of each RGB pixel was computed and then compared to a threshold value; if that grayscale value was greater than the threshold, the result pixel would be black; if the grayscale value was less than the threshold value then the resulting pixel would be white.
Edge detection
The binary image was then convolved with a Laplacian of a Gaussian operator so that zero-crossings in the 2nd derivative could be found. The Laplacian of a Gaussian used was a 9x9 pixel operator with sigma=1.125. Pixels were flagged as white if their pixel neighborhood in the 2nd derivative image contained both positive and negative values. The rest of the edge image was left black. The Hough transform was then only computed at pixel locations that were flagged as edges.
Figure 1. 9x9 Laplacian of a Gaussian with sigma=1.125.
Theta estimation
The Hough transform was modified to only consider theta values at a point if they were close to the angle represented by the gradient of the image at that location. In order to determine the angle of the gradient at a location, the sub-sampled image was blurred with a 5x5 Gaussian kernel (with sigma=1.125), and then the X and Y components of the gradient were found for each edge location by multiplying the 3x3 Sobel operators by the edge-pixel's neighborhood from the Gaussian-blurred image. The estimated angle of the gradient was then computed as the arctan(G_y/G_x) where G_y was the Y-component of the gradient, and G_x was the X-component of the gradient. The Hough transform was then optimized to consider only angles that were within some similarity threshold of the gradient angle.
Figure 2. 9x9 Gaussian with sigma=1.125.
Line Thresholding
Instead of using the traditional vote-threshold for determining which lines most likely exist in an image, the system used other knowledge of lines of latitude and longitude. If an image is landscape-aligned (so that lines of latitude extend longer than vertical lines of longitude), then lines of latitude would get more votes than lines of longitude simply because they are longer. Since both lines of latitude and longitude were to be captured, vote-thresholding was not sufficient on its own. In total, the following threshold techniques were employed.
Vote thresholding
The Hough transform first computes values R and Theta, which together define a line, at each point in the image. After these values were found for every line at every location in the image, they were put in a two-dimensional accumulator (often referred to as Hough-space). An image of Hough-space is an image whose x-axis represents different values of R and whose y-axis represents different values of Theta (from 0-Pi/2). The Hough-space image is brighter where R-Theta pairs have high counts in the source image. Vote thresholding is a way of thresholding the values of the Hough-space image to find lines that are strongly represented in the source image. A very low value for vote thresholding was used in order to filter out noise and avoid filtering short lines.
Line segment density
The traditional Hough transform detects lines, not line segments. The Hough transform was modified so that instead of just counting votes of R-Theta pairs, it kept track of all point locations on a line that contribute to the count. To store the location of each point, a value, b, is found that represents the positive and negative distance from the point on the line that is closest to the image origin. A b value of 10 means that the point is 10 pixels along the line away from the point on the line that is closest to the image origin. A b value of -10 means that the point is 10 pixels away from that closest point but in the opposite direction. The minimum and maximum b values define the start and end points of the detected line.
With this knowledge, a line segment's "density" was measured. The "density" is defined as the number of points along the line over the length of the line. Therefore ifa line contains three points, two of which are right next to each other, and the third is on the other side of the image, then the density is very low. Alternatively, if the points that contribute to the detection of a line are evenly distributed one pixel apart across the entire line, then the density is close to 1.
Line length threshold
Since bounded line segments were found, a segment length threshold was applied to reduce the amount of lines detected from noise in the source image. A line length threshold of 1/4 the minimum of the width and height of the source image was used. If the source image was 100x300 pixels then, a line had to be at least 25 pixels long to be considered.
Parallel Line Sets
In the proposed user interface, the user will select a set of parallel lines that represent the grid of latitude lines or longitude lines. The Hough transform at this stage has only detected independent lines. Since latitude and longitude lines tend to be periodic in a nautical chart, this knowledge can be used to construct a supplementary adaptation of the Hough transform. That is, a new sort of the Hough transform can be implemented: one that finds sets of parallel, evenly spaced, line sets.
In this new transform, the parameters of a set of parallel, evenly spaced, line set included the traditional R and theta parameters, as well as a third parameter, I. Here, I, represented the pixel spacing (or the difference in R values) between any two adjacent parallel lines in the line set. Unlike the traditional Hough transform, the input to the parallel line set transform was the output from the first Hough transform. So the parallel-line-Hough transform used lines as its input.
The output of the parallel-line-Hough transform is a count. For the three input values, R, theta, and I, the count represents the number of parallel lines in the set.
In implementation, the parallel-line-Hough transform looped through all of the lines that survived the above described thresholding metrics. For each line, it recorded the R value of that line and its theta value. An array of R values was constructed for each parallel line whose R value was greater than the first line. The values in the array were the differences between R values of each line with the first line. So if the first line's R value was 150 and the first parallel line has an R value of 200, then the first value in the array was 200-150=50. For each value in the array, I was set to that value. The array was looped over again, keeping a count of how many lines have an R value that is a distance from the first line's R value that is a multiple of I. This count was then kept in an accumulator (the keys for the accumulator were R, Theta, and I).
The maximum of this second Hough transform represents a group of parallel lines that has the maximum number of lines in it.
The performance of the second Hough transform is much quicker than the traditional Hough transform since values R and theta are gathered from the definition of each line, and I is computed from difference between pairs of lines. Unlike with the traditional Hough transform, this new Hough transform does not exhaustively score all parameter values.
Datasets
Test-Hough
Figure 3. A image used to test different threshold values used to filter the output of the Hough transform.
The remaining datasets were obtained through the National Oceanic and Atmospheric Association's Historical Map and Chart collection (3).
Hawaii 1794
Figure 4. Hawaii, or the "Sandwich Islands" as charted in 1794.
Boston 1857
Figure 5. Boston Harbor, 1857, Longitude measured relative to state house.
Boston 1878
Figure 6. Boston Harbor, 1878.
Pacific NW Coast 1791
Figure 7. Pacific NW Coast from 1791.
New York 1845
Figure 8. Approaches to NY Harbor, 1845.
New Jersey 1778
Figure 9. The Province of New Jersey at the time of the ratification of the United States Constitution. Longitude relative to the city of Philadelphia.
Experiments and Results
Finding a good threshold value
At first an arbitrary threshold value of 225 was chosen to create a binarized image. This value worked well for the Hawaii-1794 nautical chart, but it was too high or other images since it washed out all of the latitude and longitude lines. Closer investigation of eachimage's histogram showed that if a threshold value was chosen where the derivative of the histogram first reached above 0 before the histogram's peak, then just the white would be washed out and thelight gray lines of latitude and longitude would remain intact. Figures 10, 11, and 12 show the histograms and results of thresholding for the Hawaii andBoston datasets.
Figure 10. top: the histogram of the Hawaii 1794 image; bottom: the Hawaii 1794 image thresholded where the first derivative is above 18%.
Figure 11. top: the histogram of the Boston 1857 image; bottom: the Boston 1857 image thresholded where the first derivative is above 18%.
Figure 12. top: the histogram of the Boston 1878 image; bottom: the Boston 1878 image thresholded where the first derivative is above 18%.
Initial Hough transform results
Figure 13 (right) shows all detected lines as the result of the standard Hough transform with the Test-Hough image. Here a voting threshold of 20was used. Even though we were just computing the Hough transform at edge locations (see the approximation of the 2nd-derivative image in Figure 13 (left)), it was clear that there is nospecific vote threshold value that would return just the lines on the edges of the crosses in the image.
Figure 13. left: approximation of the 2nd derivative of Test-Hough; right: detected lines from Hough transform with vote threshold of 20.
By simply increasing the voting threshold, Figure 14 shows that the horizontal lines get filtered out of the result, but not all of the noisy horizontal lines are filtered out.This is the result of the horizontal lines being longer and therefore having more votes than the vertical lines.
Figure 14. Test-Hough with vote threshold of 60.
Introduce line-length thresholding
By computing the line-segment length of each resulting line, and thresholding based on that property instead of just the vote threshold, we can see from Figure 15 that most of the noisy lines are removed from the output. Some lines erroneously connect to the circle since at this point the Hough transform is consideringall values of Theta at each edge location.
Figure 15. Test-Hough with a line-segment length threshold of 50, vote threshold of 40.
Gradient Angle Hinting
By limiting the number of angles checked by the Hough transform, both performance and accuracy can be increased. Figure 16 shows the difference betweenthe regularHough-space image and the Hough-space image that results from only considering angles close to the gradient direction. Since the brightest spots in the Hough-space images correspond to most likely lines in the source image, Figure 16 right shows that by considering the gradient direction we have vastlyreduced the solution space without excluding any actual solutions. Furthermore, Figure 17 shows the resulting set of lines in the Test-Hough image.
Figure 16. left: The Hough-space image when considering all angles; right: the Hough-space image when only considering angles within Pi/16 of the measured gradient direction.
Figure 17. Test-Hough with a gradient-angle-deviation threshold of Pi/16. All angles considered by the Hough-transform must be within Pi/16 of the measured gradient direction.
Since the criteria formulated above works well for the simple Test-Hough case, it can now be tested on actual nautical charts. Figure 18 shows Hough-space results and lines drawn on the Hawaii-1794 image. Although a number of noisy diagonals still appear in the line image, the important pointto notice is that no lines of longitude or latitude were filtered out by the thresholds. Since the Hough-transform is going to be used as an aid to the user's selectionof latitude and longitude line sets, these results should be sufficient. In the proposed interface, the mouse can be used to cue the system to which set of lines to use.
Figure 18. top: the Hough-space image for lines in the Hawaii-1794 image; bottom: lines drawn on the Hawaii-1794 image's edges
Constructing a User Interface
When the user first starts the system, the source image is displayed and the user is cued to highlight sets of latitude lines (Figure 19). As the user mouses over different lines in the image,sets of periodic (evenly spaced and contiguous) parallel lines are highlighted (Figure 20). When the user is satisfied that the displayed lines highlight at least two lines of latitude in thenautical chart, the user clicks to begin labeling. Labeling can be completed with help from the system. Since the end points of each line segment are known, images of the neighborhoods ofeach end point can be constructed since they should overlap with the margins (and thus labels) of the lines of latitude. Labeling the latitude should then just bean issue of transcribing what appears on the margins (Figure 21). If the end-point images do not contain the latitude labels, the user can mouse over the original image to find the appropriate place in the margin. By hitting the "NEXT" button the user saves his results and moves on to the next line. By hitting "SKIP" the user ignores a given lineand that line does not get saved to the list of geo-references. Once this process is completed for lines of latitude, it is repeated for lines of longitude (Figure 22).
Figure 19. The user is cued to highlight sets of latitude lines
Figure 20. The user mouses over a set of latitude lines and they are highlighted by the system.
Figure 21. The user can label each line of latitude with help from the sub-images around the end point of each line of latitude. If the sub-images to not containlabels for the lines of latitude or longitude, the user can mouse over any area of the image to see a magnified view of that area.
Figure 22. The labeling process is repeated for lines of longitude.
Figure 23. Once lines of latitude and longitude are labeled, the geo-references can be saved.
The saved geo-references define a mapping from the source image coordinate space to world latitude-longitude space. Using this information the original image can be reprojectedso that it can be displayed using mapping toolkits like Google Maps. Using this tool to create geo-references, the Boston-1878 image was reprojected so that it could be displayedin Google maps (Figure 24). The value in doing this is illustrated by Figure 25, as it is now possible to compare Boston's coastline in 1878 with today's coastline as represented byGoogle's satellite layer (Figure 25).
Figure 24. The Boston-1878 image reprojected and displayed in Google maps.
Figure 25. Decreasing the opacity of the Boston-1878 image from top (100%) to bottom (0%) as it is displayed on top of Google's satellite layer. The development of South Boston, Winthrop, and the airport are made evident by this comparison.
Results with each dataset
For each dataset, success was measured by how easy it is to label the lines of latitude and longitude. The ease at which a chart can be labeled was decomposed into two operations:selecting parallel line sets, and labeling them. The ease of selecting parallel line sets was defined as the number of visible lines that when moused over, subsequentlyhighlighted a correct set of two or more lines of latitude or longitude. So if a set of four lines of latitude that were visible on a chart could only be highlighted whenone line was moused over, then its measure of ease would be 1/4=25%. A "correct" line set is one that contains at least two lines of latitude or longitude thatoverlapped with real lines of latitude or longitude in the image. Theease of labeling was defined as the percent of lines in a line set whose end point images contain latitude or longitude labels. Finally, the overall ease of labelinga nautical chart was defined as the multiplication of all measures.
| Dataset | Ease of selecting latitude | Ease of labeling latitude | Ease of selecting longitude | Ease of labeling longitude | Ease |
| Hawaii 1794 | 25% | 75% | 80% | 80% | 12% |
| *Boston 1857 | 87.5% | 87.5% | 75% | 100% | 57.4% |
| Boston 1878 | 100% | 37.5% | 81.8% | 17% | 5.2% |
| Pacific NW Coast 1791 | 100% | 80% | 100% | 66.7% | 53.4% |
| New York 1845 | 0% | 0% | 0% | 0% | 0% |
| *New Jersey 1778 | 33% | 0% | 0% | 50% | 0% |
Discussion
The traditional Hough transform alone, with vote thresholding, was not sufficient to detect lines of latitude and longitude in nautical charts. It was clear from early experimentsthat the nature of high-resolution nautical charts with very thin and often faint lines of latitude and longitude would be an obstacle to automatic detection. Since lines occur anywherein an image and since these lines are oriented with virtually any angle (some nautical charts are oriented diagonally so that their lines of latitude and longitude are diagonal),it was clear that a semi-automatic approach to detecting and labeling lines must be taken. For this reason the output from the Hough transform was used in conjunction with user input(via mouse gestures) to select lines.
The interface proposed and designed here worked well on average. Table 1 shows how easy, in some cases, it was to highlightlatitude and longitude. The Hawaii, Boston, and Pacific Northwest charts performed particularly well when it came selecting lines of latitude and longitude. Of these,the 1857 Boston chart fared the best since even labeling was relatively easy. On the other hand, the line segments detected for latitude and longitude in theHawaii and 1878 Boston datasets tended to not extend to the chart's margins, or extended past them, causing the end point images to not be helpful. The user interface could be adapted to generate the end point images from the intersection of each line with the chart's margins instead of just the line segment's end points. This mightmake it easier to label latitude and longitude, but it would introduce another problem: how to find the chart's margin.
The New York and New Jersey datasets fared particularly badly in their easiness metrics. The New York image had a very low contrast ratio which made it hard for the user to evensee where lines of latitude or longitude were. The New Jersey image, though not actually a nautical chart, had very few lines of latitude and longitude to serve as candidates tothe Hough transform. For this reason, it scored a low 0% on the selection of lines of longitude. In fact, only one line of longitude could be seen by the Hough transform.
For cases where the Hough transform provided inadequate candidate lines for the user's selection, the user interface could be adapted to allow drawing any line over theimage. This change would make the tool more usable in the case of very difficult datasets.
Labeling scores of easiness generally were the difficult metric for every dataset. A low score for labeling latitude and longitude in the Boston-1878 image brought its generalscore of ease down to 5.2%. While this score could be improved by finding the image's margins and providing a sub-image at that location to the user instead of the sub-image arounda line segment's end points, another improvement and line of possible research would be incorporating an optical character recognition (OCR) algorithm. The OCR algorithm could attemptto read the degrees, minutes, and seconds off the margin of the chart and simply ask for the user to verify those coordinates. To do this would be out of the scope of this project but would be an exciting next step.
Conclusion
The tool designed in this project would be an invaluable time-saver to the manual alternative of georeferencing very large nautical charts by hand. Historical nautical charts posemany obstacles to the ease at which they can be displayed in modern mapping toolkits like Google Maps. They are extremely high resolution (requiring so much memory that the entire filecannot be opened at once), they have no standard color encoding to differentiate land from graticule and margins, they have varied orientations, varied scales,and often exhibit artifacts from scanning (paper creases, tears, discolorations, etc.).
Here the Hough transform can be shown to serve as a valuable aid to the user's annotation of these charts. In just a few minutes accurate geo-references can be computedfor each dataset, which can then be used to project these charts into a geographic display. Doing this for historical charts is not only an interesting thing to do, but it is something thathas never been done for such imagery.
References
- (1) Duda, R. O. and P. E. Hart, "Use of the Hough Transformation to Detect Lines and Curves in Pictures," Comm. ACM, Vol. 15, pp. 11-15 (January, 1972)
- (2) Hough Transform [http://en.wikipedia.org/wiki/hough_transform - April 8, 2010]
- (3) NOAA Historical Map and Chart Project [http://historicalcharts.noaa.gov/historicals/historical_zoom.asp - April 8, 2010]
| Comments |
|