Patent Description:
The present invention relates to the generation of shadows in images, and, more particularly, to methods and apparatus for performing shadowing on the basis of very large <NUM>-dimensional (<NUM>. 5D) elevation data sets.

Shadow mapping over an illuminated terrain is an inherent part of topographic relief shading techniques. Unless blocked from the light source, a given grid point on a surface is subject to a direct beam radiance defined by a local normal to the surface at the given grid point and local slope and aspect angles, which can be estimated using as small a local area as a <NUM>×<NUM> pixel kernel around the given grid point.

Preliminary Definitions: As used herein, and in any appended claims, the following terms shall have the meanings indicated, unless the context dictates otherwise:.

In areas of complex topography, geometry mandates that one must know the sky view in all directions for a point in order to determine if a given sun vector (specified, for example, in coordinates of (azimuth, elevation), referred to, herein, also as (AZ, EL)) is in view from that point, or not. Techniques to compute viewshed for a topographic point as a polar plot of radial coordinate for elevation and tangential coordinate for azimuth are designed in the paper by <NPL>). A schematic process to generate point-based viewshed maps designed by the authors is shown in <FIG> shows a digital elevation model (DEM) of a scene centered about a specified grid point <NUM>. <FIG> is an elevation cross section of the skyview relative to grid point <NUM>, while <FIG> is an azimuth-elevation (AZ-EL) representation of the viewshed relative to grid point <NUM>.

In addition to the absence of full 3D geometric data, another daunting feature of GIS data sets is the very large quantity of data that must be processed, often in a near real-time mode. The initial process to generate a viewshed for each point is computationally very expensive; however, the following shadow map generation for any given time of day (or intervals of time) is very fast. The accuracy limitation of the method is defined by a discrete number of azimuth directions selected to generate a viewshed profile in polar coordinates: interpolating between just <NUM> directions is obviously of lower fidelity than, say, <NUM> directions. The paper "<NPL>et al. discloses a real-time soft shadowing method applicable for height maps.

In accordance with embodiments of the present invention a computer-implemented method is provided for creating a shadow mask representing irradiation of a physical scene based upon a plurality of elevation grid points of an array of elevation data associated with the physical scene, as defined in claim <NUM>. The method has steps of:.

In accordance with other embodiments of the present invention, there may an additional step of calculating an irradiance contribution to a total insolation incident upon a specified grid point based on the irradiance of the light source at the elevation grid point as modulated by the binary shadow mask.

In further embodiments of the present invention, the light source may be the Sun.

In yet further embodiments of the present invention, propagating a shadow horizon may have a further step of applying a transformation to the direction of the light source to reflect curvature of an underlying surface of the physical scene. The underlying surface of the physical scene may be substantially spherical, and may, more particularly, be the Earth.

In accordance with another aspect of the present invention, a computer program product is provided for use of a computer system for generating a shadow mask representing irradiation of a physical scene based upon a plurality of elevation grid points of an array of elevation data associated with the physical scene, as defined in claim The computer program product has a non-transitory computer-readable medium on which are stored computer instructions such that, when executed by a processor, the instructions cause the processor to:.

In other embodiments, the instructions may also cause the processor to calculate an irradiance contribution to a total insolation incident upon a specified grid point based on the irradiance of the light source at the elevation grid point as modulated by the binary shadow mask. The instructions may further cause the processor to apply a transformation to the direction of the light source to reflect curvature of an underlying surface of the physical scene.

The invention will be more fully understood by referring to the following Detailed Description of Specific Embodiments in conjunction with the Drawings, of which:.

Further Definitions: As used herein, and in any appended claims, the following terms shall have the meanings indicated, unless the context dictates otherwise:.

Three components contribute to the illumination at each elevation grid point: direct beam irradiance, diffuse radiation (skylight), and ground albedo (neighbourhood terrain illuminates any non-horizontal surface).

In accordance with embodiments of the present invention, methods are provided for generating a shadow mask in a manner that may be, advantageously, both memory-efficient and computationally fast. Many methods for creating shadow maps may be practiced in accordance with the prior art. The present invention pertains to various novel methods that are now described and claimed herein.

A method of shadow estimation, in accordance with one embodiment of the present invention, is now described with reference to <FIG>, where it is designated generally by numeral <NUM>. It is to be understood that steps in accordance with the method need not be performed in the sequence in which they are shown, or described. A step that is first - as a matter of enumeration, is not necessarily temporally first, for example.

Elevation point data are acquired using aircraft or satellites or any other means, and it is assumed here, for heuristic convenience, that the data are received by a computer system <NUM> (shown in <FIG>) and stored in a memory device <NUM> in an array designated generally by numeral <NUM> (shown in <FIG> in which rows correspond to the East-West direction (<NUM>° and <NUM>° azimuth) and columns correspond to the North-South direction (<NUM>° and <NUM>° azimuth). It is to be understood that the arbitrary assignment of an orientation to the data array is a heuristic convenience, and that the rotation of the array about the anchor point and corresponding modification of the algorithm described herein is within the competence of a reader of ordinary skill. Array <NUM> may also be referred to herein as "grid of points" <NUM>, or as "grid" <NUM>, and it is to be understood that, within the scope of the present invention, various points might be missing from array <NUM> and that the methods described and claimed herein may still be performed.

In order to begin a process of shadow mask generation, designated generally by numeral <NUM>, elevation data are acquired, received and stored (<NUM>) in matrix format, which is to say that data are characterized by two indices indicating orthogonal cardinal directions, which, for heuristic purposes alone, are described herein as being North-South and East-West directions. Using the same orientation convention, a quadrant containing a specified light vector azimuth is identified (<NUM>). A "quadrant," otherwise referred to herein as a "sector," refers to a division of the azimuth circle into four parts, where azimuth is conventionally defined in a clockwise sense relative to <NUM>° at a reference North direction. The azimuthal circle is a circle of <NUM>° surrounding a specified anchor point. A "light vector azimuth" is the value of the azimuth component of a light vector, as defined above. The four quadrants (or sectors), respectively centered about the four cardinal directions, are the following:.

Based on the sector of the light vector azimuth, one of four shadow frontline propagation modes is selected (<NUM>):.

As used herein and in any appended claims, the term "shadow frontline" shall refer to a row (or column, as the case may be) that is currently subject to computational evaluation as to whether a given point blocks a light vector from a specified source of illumination to a specified anchor point.

The "first line" corresponds to the extremal row or column (as the case may be) of the array <NUM> of elevation data toward the source of illumination. Thus, in the case of the source of illumination lying in the NW-NE quadrant, the first line is the top row <NUM> (referred to, hereinafter, as the first line, for heuristic convenience) of array <NUM>. The shadow frontline is, initially, the first line. In step <NUM> (shown in <FIG>) each point <NUM> of the first line is flagged as "illuminated" or "lit", and the elevation associated with each point <NUM> is marked as the "light-blocking height" in the shadow frontline.

Points <NUM> of the shadow frontline are projected (step <NUM>) along light vector <NUM> onto the successive line <NUM> of array <NUM> and recorded. "Projected," as the term is used herein, means the connection of an elevation of the point <NUM> with the elevation of a point <NUM> in the successive line <NUM> by a line that is parallel to the light vector <NUM>. In the successive line <NUM>, light-blocking heights are recorded (<NUM>) as the greater of the elevation of each grid point in the elevation data or the projection of the elevation of a point in the preceding line <NUM>. Grid points may be connected by a spline fit, such as a Hermitian spline, generating a curve <NUM> (shown in <FIG>), referred to herein as a "shadow frontline" (defined above), separating light from dark, looking toward the direction of the illuminating source. The shadow frontline is stored (<NUM>) in a frontline buffer of computer memory, and recursively updated in successive lines (rows or columns). Propagation of a shadow to successive lines (rows or columns, as the case may be, depending on the quadrant in which the illuminating source appears) results in a shadow frontline, described in detail below with reference to <FIG>. Recording the maximum heights in successive lines may be referred to as "collecting a new shadow frontline. " Propagation of blocking heights is then repeated to successive lines until a final line <NUM> of array <NUM> is reached. Each grid point is thus marked by an elevation below which light from the source of illumination is blocked.

Frontline propagation, as heretofore described, may advantageously provide for memory-efficient and computationally fast generation of a shadow mask since input images need to be read line by line only (row-by-row or column-by-column) and the output mask image is written line-by-line without requiring an area cache. Once a given row or column is processed, the input line and output line can be removed from active memory and never accessed again; only a single shadow frontline <NUM> is retained, as stored in the frontline buffer.

The contour of a shadow frontline <NUM> is now described with reference to <FIG>, where solid ringed dots <NUM> represent sunlit points at the second column <NUM>, solid blue dots <NUM> represent shadow points, and heavy open circles <NUM> represent the light-blocking heights collected into the shadow frontline <NUM> at the current step. Light open circles <NUM> represent the previous shadow frontline projected to the current line along the direction of light vector <NUM> (shown in <FIG>), namely (-<NUM>°, -<NUM>°) of horizontal and vertical slope, in the current example. Spline curve <NUM> represents the "horizon line" separating light from dark at the current column looking west (to the Left, in <FIG>). The terms "horizon line" and "shadow horizon" are entirely synonymous with "shadow frontline," as the term has been defined above. All elevated observers above the horizon line see the light of the illuminating source, whereas observers at positions <NUM>, <NUM>, and <NUM> in <FIG> are in the shadow, by virtue of the fact that their elevation is beneath the shadow frontline.

In accordance with other aspects of the present invention, shadow casting employs approximations with polynomial models of the rigorous 3D transforms between geodetic coordinates (latitude, longitude, height) or map-projected elevation grids (easting, northing, height) and anchored topocentric Cartesian grids, but, instead of applying the transform to the entire scene, the transform is applied to the light vector <NUM>.

While a rigorous 3D transformation relating geodetic or map-projected data to topocentric data requires numerous nonlinear transformational terms, polynomial approximations may be used, whether continuous or piecewise smooth, within the scope of the present invention, and are notoriously well-known and commonly used to speed up XYZ transform of data consisting of massive numbers of points in three dimensions. This approach is used, in accordance with embodiments of the present invention, to model angular corrections to the direction of the local light vector <NUM>.

As an example, if grid of points <NUM> (shown in <FIG>) spans <NUM> degrees in latitude and <NUM> degrees in longitude at some mid-latitude anchor point (such as mid-Europe, or such as along the <NUM>th parallel North that forms the US-Canada border), then the overall shape of array <NUM> in a topocentric Cartesian frame is as shown in <FIG> in horizontal <NUM> and in <FIG> in vertical <NUM>.

Both rows and columns are not straight lines any more, with rows more curved than columns. Nevertheless, one can assume that a line segment between two immediately neighboring grid nodes is indeed a straight line, and there is only a need to correct the light direction (i.e., the direction of light vector <NUM>, shown in <FIG>) emanating from a grid node of the previous row or column with respect to horizontal and vertical orientation of the line segment connecting two points on both sides of the light ray as it intersects the next row or column.

As far as a single cell-rectangle is concerned, the process of light ray emanating from a corner and bisecting the opposite edge with a given ratio needs only a single correction to that ratio. All the local shape deformations such as differential scaling, rotation, shear, etc. cumulatively contribute to a single change in ratio at which the opposite edge is bisected - one correction to a horizontal ratio and one correction to a vertical elevation angle due to earth curvature.

These corrections as functions of the grid node (i ,j) position in the matrix can be modelled by special bi-variate polynomials that retain only those terms that significantly contribute to the correction shape. The special polynomials that are used and their parameterization, based on L<NUM> norm minimization or otherwise, are very well - known in the art, and a rigorous description of their use need not be described here, but is left to the design choice of a reader of ordinary skill, and may be found in "<NPL>, Any of a multitude of specific polynomial functions, whether continuous or piecewise smooth, may be selected, as a matter of design choice depending on the vertical /horizontal datum for a given elevation grid, and all are within the scope of the present invention. Geodetic grids obviously would have different correction forms than the map-projected grids, and the transverse Mercator projections may have a shape different from the conic Lambert projections, etc., though all are within the scope of the present invention as claimed. In accordance with embodiments of the present invention, the transformation that would be applied to the coordinate frame is applied to the incident light vector <NUM> as the shadow casting process described above is performed.

One example of a method for correcting the direction of light vector <NUM> is now described. For each row (assuming that frontline propagation is proceeding row-wise, and, else, mutatis mutandis, each column), a single polynomial is used for all points (x, y, z) relative to a fiducial anchor point (X, Y, Z). Discrepancies are collected, <MAT> where (Fx, Fy) is the projection from geodetic (or Universal Transverse Mercator (UTM), etc.) onto local topocentric Cartesian coordinate system, otherwise referred to herein as "mapping axes.

A single polynomial is used to model each transverse dimension, such as <MAT> While the foregoing <NUM>-term trivariate polynomial is provided by way of example, any model function may be used within the scope of the present invention. The polynomial that is employed serves to model xtopo = xUTM + dx at any point in three-dimensional UTM space (inside the grid <NUM>, on the grid nodes, as well as outside), although the approximation accuracy degrades outside the grid cube used to fit the model. Other polynomial terms may be retained as a matter of design choice, based on the source datum - geodetic, UTM, projected conic Lambert, etc..

The parameter vector of coefficients px = [co,. c<NUM>o] is defined as a solution to a least squares minimization problem in which {B* px - Fx(X,Y,Z)- x} is minimized using standard numerical techniques. Thus, for each pixel in the original raster elevation grid <NUM>, the correction to the bisection ratio is found for the succeeding line in the frontline propagation described above.

Referring to <FIG>, a "bisection ratio" is defined for a conformal mapping (a holomorphic mapping that preserves angles locally, such as the mapping from geodetic to topocentric coordinates) as follows: For an angle φ that bisects an opposing line segment BC in the untransformed coordinate system (i. e, in the domain of the conformal mapping), the bisection ratio is the ratio of the length of the segment B'D' to the length of the segment B'C' in the transformed frame, where the ratio is denoted as α, as shown in <FIG>.

While, in the original rectangular grid, given in, say, UTM coordinates, the bisection ratio is constant - it is α versus (<NUM>-α) - for any two adjacent points on any column, in the transformed grid, becoming a quadrilateral mesh, the bisection ratio varies with pixel position.

As far as the elevation EL of the light vector <NUM> is concerned, it is corrected by applying a polynomial correction that is second-degree in the distance from the anchor point, where the coefficients have been determined by minimizing residuals to the discrepancy between geodetic and topocentric Cartesian coordinates relative to the anchor point.

Various aspects of the invention may also be implemented as specialized software executing in a general-purpose computer system <NUM> such as that shown in <FIG>. Alternatively, aspects of the present invention are advantageously implemented on mobile phones, tablets, and other devices. The computer system <NUM> may include a database server <NUM> connected to one or more memory devices <NUM>, such as a disk drive, memory, or other device for storing data. Database server <NUM> stores point elevation data comprising elevations at a plurality of points, or other data to which the present invention may be applied. A processor <NUM> (also referred to herein as a central processing unit, or processing server) <NUM> contains computer-executable software configured to cast a shadow as taught in the foregoing description. Alternatively, a computer program product may have a non-transitory computer-readable medium on which are stored computer instructions such that, when executed by processor <NUM>, the instructions cause the processor to cast a shadow as taught in the foregoing description. Memory <NUM> is typically used for storing programs and data during operation of the computer system <NUM>. Components of computer system <NUM> may be coupled by an interconnection mechanism <NUM>, which may include one or more buses (e.g., between components that are integrated within a same machine) and/or a network (e.g., between components that reside on separate discrete machines). The interconnection mechanism <NUM> enables communications (e.g., data, instructions) to be exchanged between system components of system <NUM>. Computer system <NUM> also includes one or more input devices <NUM>, for example, a keyboard, mouse, trackball, microphone, touch screen, and one or more output devices <NUM>, for example, a printing device, display screen, speaker. In addition, computer system <NUM> may contain one or more interfaces (not shown) that connect computer system <NUM> to a communication network (in addition or as an alternative to the interconnection mechanism).

The computer system may include specially-programmed, special-purpose hardware, for example, an application-specific integrated circuit (ASIC). Aspects of the invention may be implemented in software, hardware or firmware, or any combination thereof. Further, such methods, acts, systems, system elements and components thereof may be implemented as part of the computer system described above or as an independent component.

Although computer system <NUM> is shown by way of example as one type of computer system upon which various aspects of the invention may be practiced, it should be appreciated that aspects of the invention are not limited to being implemented on the computer system as shown in <FIG>. Various aspects of the invention may be practiced on one or more computers having a different architecture or components than that shown in <FIG>.

Processors <NUM> and operating systems employed in conjunction with servers <NUM> and <NUM> define a computer platform for which application programs in high-level programming languages are written. It should be understood that the invention is not limited to a particular computer system platform, processor, operating system, or network. Also, it should be apparent to those skilled in the art that the present invention is not limited to a specific programming language or computer system. Further, it should be appreciated that other appropriate programming languages and other appropriate computer systems could also be used.

One or more portions of the computer system may be distributed across one or more computer systems (not shown) coupled to a communications network. These computer systems also may be general-purpose computer systems. For example, various aspects of the invention may be distributed among one or more computer systems configured to provide a service (e.g., servers) to one or more client computers, or to perform an overall task as part of a distributed system. For example, various aspects of the invention may be performed on a client-server system that includes components distributed among one or more server systems that perform various functions according to various embodiments of the invention. These components may be executable, intermediate, or interpreted code which communicate over a communication network (e.g., the Internet) using a communication protocol (e.g., TCP/IP).

It should be appreciated that the invention is not limited to executing on any particular system or group of systems. Also, it should be appreciated that the invention is not limited to any particular distributed architecture, network, or communication protocol.

Having now described some illustrative embodiments of the invention, it should be apparent to those skilled in the art that the foregoing is merely illustrative and not limiting, having been presented by way of example only. Numerous modifications and other illustrative embodiments are within the scope of one of ordinary skill in the art and are contemplated as falling within the scope of the invention, that is defined by the accompanying claims. In particular, while descriptions have been provided in terms of shadow casting, they are not limited to this context. The procedures are applicable to a wide variety of illumination rendering areas.

Moreover, where examples presented herein involve specific combinations of method acts or system elements, it should be understood that those acts and those elements may be combined in other ways to accomplish the same objective of shadow mapping. Acts, elements and features discussed only in connection with one embodiment are not intended to be excluded from a similar role in other embodiments.

Use of ordinal terms such as "first", "second", "third", etc., in the claims to modify a claim element does not by itself connote any priority, precedence, or order of one claim element over another or the temporal order in which acts of a method are performed, but are used merely as labels to distinguish one claim element having a certain name from another element having a same name (but for use of the ordinal term) to distinguish the claim elements. Additionally, single device features may fulfill the requirements of separately recited elements of a claim.

Claim 1:
A computer-implemented method for creating a shadow mask representing irradiation of a physical scene based upon a plurality of elevation grid points of an array (<NUM>) of elevation data associated with the physical scene, the method comprising:
a. receiving elevation point data in which an elevation value is associated with each of a plurality of spatial coordinates;
b. determining, for a specified anchor point, a direction to a light source characterized by an irradiance value at each elevation grid point, wherein the anchor point is a given point with respect to which directions are determined;
c. associating with the light source a quadrant of an azimuthal circle;
d. selecting a shadow frontline propagation mode based on the quadrant associated with the light source;
e. propagating a shadow horizon along the array of elevation data by projection onto successive lines (<NUM>) of elevation data;
f. recursively updating the shadow horizon stored in a frontline buffer of computer memory; and
g. generating a binary shadow mask based on whether a specified elevation point lies above or below the shadow horizon;
wherein a single shadow horizon is stored in the frontline buffer.