WO2012169294A1 - Ｄｔｍ推定方法、ｄｔｍ推定プログラム、ｄｔｍ推定装置及び３次元建物モデルの作成方法、並びに、領域抽出方法、領域抽出プログラム及び領域抽出装置 - Google Patents

Ｄｔｍ推定方法、ｄｔｍ推定プログラム、ｄｔｍ推定装置及び３次元建物モデルの作成方法、並びに、領域抽出方法、領域抽出プログラム及び領域抽出装置 Download PDF

Info

Publication number
WO2012169294A1
WO2012169294A1 PCT/JP2012/061204 JP2012061204W WO2012169294A1 WO 2012169294 A1 WO2012169294 A1 WO 2012169294A1 JP 2012061204 W JP2012061204 W JP 2012061204W WO 2012169294 A1 WO2012169294 A1 WO 2012169294A1
Authority
WO
WIPO (PCT)
Prior art keywords
region
dtm
area
data
value
Prior art date
Application number
PCT/JP2012/061204
Other languages
English (en)
French (fr)
Inventor
純一 須▲崎▼
Original Assignee
国立大学法人京都大学
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Priority claimed from JP2011129316A external-priority patent/JP2012256230A/ja
Priority claimed from JP2011144266A external-priority patent/JP2013012034A/ja
Application filed by 国立大学法人京都大学 filed Critical 国立大学法人京都大学
Priority to US14/122,082 priority Critical patent/US20140081605A1/en
Publication of WO2012169294A1 publication Critical patent/WO2012169294A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/13—Architectural design, e.g. computer-aided architectural design [CAAD] related to design of buildings, bridges, landscapes, production plants or roads
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/05—Geographic models
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/50—Depth or shape recovery
    • G06T7/521—Depth or shape recovery from laser ranging, e.g. using interferometry; from the projection of structured light
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10028—Range image; Depth image; 3D point clouds
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10032—Satellite or aerial image; Remote sensing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30181—Earth observation

Definitions

  • the present invention relates to a technique for estimating a DTM (Digital Terrain Model) based on surface laser scanner data from an aircraft.
  • DTM Digital Terrain Model
  • the present invention also relates to a method, a program, and an apparatus for extracting a building area based on photograph data taken from an aircraft or an artificial satellite.
  • LiDAR Light Detection And Ranging
  • a DSM Digital Surface Model
  • the top surface of the building and the top end of the tree can be obtained for the desired region in addition to the ground surface (the ground surface such as a site or road).
  • the DTM representing only the ground surface can be estimated from the partially obtained ground surface data and the data obtained by interpolating the entire data.
  • the height of a building can be calculated as (DSM-DTM) (see Non-Patent Document 1, for example).
  • DTM the lowest elevation point in the predetermined range of LiDAR data is the ground surface. Then, if the ground surface is searched within the range and is interpolated, the DTM can be estimated.
  • Non-Patent Document 1 Feature extraction and area extraction in the case of bird's-eye view have been important themes from the past (for example, Non-Patent Document 1).
  • Non-Patent Documents 3 and 4 road extraction from aerial photographs or artificial satellite images has been actively studied (for example, Non-Patent Documents 3 and 4).
  • Remote sensing image classification methods can be broadly classified into pixel-based classification and area (object) -based classification.
  • Pixel-based classification techniques such as clustering (so-called unsupervised), maximum likelihood (supervised), and Support Vector Machines (SVM) (supervised) assign each pixel to a classification class based on statistical theory.
  • SVM Support Vector Machines
  • the region-by-region classification method information around the pixel (Context) is used.
  • classification using a mathematical Morphological approach is a typical example (for example, Non-Patent Documents 6 to 10).
  • 22 (a) and 22 (b) are examples of aerial photographs of an urban area (Higashiyama Ward, Kyoto City) where houses are densely packed. In such a dense urban area, the conventional region extraction method does not function effectively.
  • the factors first, (a) because the buildings are close to each other, if the brightness value and texture of the roof are similar, the boundary line of the building will not be clear.
  • roofs with corrugated cross-sectional shapes for example, slate roofs
  • DTM estimation method Since the above DTM is an estimation unlike LiDAR data, it is unavoidable that some errors and misrecognitions occur, but how to suppress the errors and misrecognitions is important. As a result of DTM estimation on various ground surfaces, it is difficult to recognize data points such as river beds, dikes, bridges, etc. as ground surfaces, especially when there are rivers on flat land as a whole. It has been found. In addition, the rooftop parking lot (the top of the building) that gradually rises from the road may be mistakenly recognized as the ground surface.
  • an object of the present invention is to improve the accuracy of DTM estimated based on LiDAR data.
  • the above factor (a) gives a problem of how to deal with the case where the edge of the region is not sufficiently obtained. It is also a problem to exclude unnecessary edges from edges that are excessively extracted due to the above factor (c).
  • the problem (c) is related to the problem (b).
  • Canny proposes an edge extraction operator that is resistant to noise. This is one of the typical edge extraction operators widely used so far.
  • Non-Patent Document 12 a method of removing noise using Wavelet (for example, Non-Patent Document 12) or utilizing an average value (for example, Non-Patent Document 13) or a multi-resolution image (for example, Non-Patent Document 14) has been reported. ing.
  • Non-Patent Documents 15 to 18 In addition, in connection with the above factor (b), research cases for restoring the luminance value of the shadow region have been reported (for example, Non-Patent Documents 15 to 18). However, as far as the restoration results in Non-Patent Document 18, for example, the boundary between the original shadow region and the non-shadow region appears clearly in the restored image, and the correction value in the shadow region Problems such as difficult decisions and sometimes overrecovering are still unsolved.
  • an object of the present invention is to provide a technique for more accurately extracting a building area based on data such as aerial photographs.
  • the present invention is a DTM estimation method for estimating, using a computer, DTM that is altitude data of only the ground surface for a predetermined range based on laser scanner data on the ground surface by an aircraft.
  • DTM that is altitude data of only the ground surface for a predetermined range based on laser scanner data on the ground surface by an aircraft.
  • a river region is extracted by connecting pixels in which no data exists in the unit grid, and a temporary DTM is estimated by setting a first maximum allowable slope value for data excluding the river region, and from the estimated DTM
  • a local inclination is calculated, and when the inclination exceeds a predetermined value, a second maximum allowable inclination value larger than the first maximum allowable inclination value is set and the DTM is estimated again.
  • the DTM estimation method as described above, it is possible to prevent the accuracy of the DTM from being deteriorated by searching the ground surface across the river by removing the river region in estimating the DTM. Also, when the estimated DTM locally has a slope exceeding a predetermined value, a second maximum allowable slope value that is larger than the first maximum allowable slope value is set and the DTM is estimated again, so that the first time The surface data that could not be searched can be added to improve the accuracy of DTM estimation.
  • the second maximum allowable slope value may be increased or decreased according to the features of the topography and the streets.
  • optimal DTM estimation can be performed according to the features of the topography and the streets.
  • the present invention is a method of creating a three-dimensional building model including the DTM estimation method of (1) above, wherein the laser scanner data is a DTM estimated by the DTM estimation method, a non-ground surface
  • This is a method of creating a three-dimensional building model that is separated into data, and for each region in the non-ground surface data, the presence or absence of a normal vector in a pairing direction is examined, and the shape of the roof is estimated based on the result.
  • accurate non-surface data is acquired based on the DTM estimated with high accuracy, and further, the presence of a normal vector in the direction of the pair is checked, so that in addition to a flat roof, gables, dormitories, etc.
  • a three-dimensional building model can be created by estimating the shape of the roof.
  • each region is extracted in advance by preferentially extracting at least a rectangular shape from a predetermined range of aerial photo data. May be.
  • a rectangle having a high probability as the shape of the roof it is possible to increase the accuracy of identifying the roof region and create a more accurate three-dimensional building model.
  • the present invention is a DTM estimation program (or DTM estimation program storage medium) that estimates DTM, which is altitude data only on the ground surface, for a predetermined range based on laser scanner data on the ground surface by an aircraft.
  • DTM which is altitude data only on the ground surface
  • This is a DTM estimation program (or DTM estimation program storage medium) for realizing the function of estimating the above by a computer.
  • DTM estimation program (or DTM estimation program storage medium) as described above, it is possible to prevent the DTM accuracy from deteriorating due to the search of the ground surface across the river by excluding the river area in estimating the DTM. Also, when the estimated DTM locally has a slope exceeding a predetermined value, a second maximum allowable slope value that is larger than the first maximum allowable slope value is set and the DTM is estimated again, so that the first time The surface data that could not be searched can be added to improve the accuracy of DTM estimation.
  • the present invention is a DTM estimation device for estimating DTM, which is altitude data only on the ground surface for a predetermined range, based on laser scanner data on the ground surface by an aircraft, and includes laser scanner data on the ground surface by an aircraft.
  • a provisional DTM is estimated by setting an inclination value, and a local inclination is calculated from the estimated DTM. If the inclination exceeds a predetermined value, a second maximum allowable value that is greater than the first maximum allowable inclination value.
  • a DTM estimation unit that sets a slope value and estimates DTM again, and a display unit that displays the estimated DTM are provided.
  • the river area is excluded in the estimation of the DTM, so that the accuracy of the DTM can be prevented from being deteriorated by searching the ground surface across the river. Also, when the estimated DTM locally has a slope exceeding a predetermined value, a second maximum allowable slope value that is larger than the first maximum allowable slope value is set and the DTM is estimated again, so that the first time The surface data that could not be searched can be added to improve the accuracy of DTM estimation.
  • the present invention is an area extraction method for extracting an area of a building using a computer based on data of a photograph taken from an aircraft or an artificial satellite, and sets a plurality of different discretization widths, Regarding the discretization width, the luminance value of the data is discretized into a plurality of values discretely set with the discretization width, and pixels having the same value in the discretized image obtained by discretization are connected, A region having a shape close to a rectangle is extracted as a candidate, and a region having a shape closer to a rectangle is adopted as a region of a building among the plurality of region groups extracted for each of the plurality of different discretization widths. To do.
  • a region having a texture with a large luminance value dispersion can be extracted as one region by discretization. Further, by connecting points having the same value and extracting a region having a shape close to a rectangle, it becomes easy to remove regions other than buildings. In addition, by adopting a region of a more rectangular shape as a building region among a plurality of region groups extracted for a plurality of different discretization widths, a locally optimal spatial smoothing parameter A function equivalent to that of applying can be exhibited.
  • a rectangular index defined by (area of area / area of rectangle surrounding area) can be used as an index representing a shape close to a rectangle.
  • the degree of “closeness” in a shape close to a rectangle can be expressed simply and accurately as an index.
  • the merger can be executed when the regions are closer to a rectangle. It may be. In this case, the area of the building can be captured more accurately.
  • a region that is estimated to be vegetation may be excluded from extraction targets based on RGB luminance values.
  • the region can be excluded from the extraction target in consideration of the color characteristics of the vegetation.
  • the present invention is an area extraction program (or area extraction program storage medium) that extracts a building area based on data of a photograph taken from an aircraft or an artificial satellite, and has a plurality of different discretization widths.
  • a function for discretizing the luminance value of the data into a plurality of values discretely set by the discretization width, and pixels having the same value in the discretized image obtained by discretization for each discretization width A function of connecting and extracting a region of a shape close to a rectangle as a candidate for a region of the building, and a region of a shape closer to a rectangle out of a plurality of region groups extracted for each of the plurality of different discretization widths,
  • the present invention is an area extraction device that extracts a building area based on data of a photograph taken from an aircraft or an artificial satellite, and sets a plurality of different discretization widths, and each discretization width is set.
  • a function of discretizing the luminance value of the data into a plurality of values discretely set with the discretization width, and connecting pixels having the same value in the discretized image obtained by discretization A function of extracting an area having a shape close to a rectangle as a candidate, and a function of adopting an area having a shape closer to a rectangle out of the plurality of area groups extracted for each of the plurality of different discretization widths, as a building area It is what has.
  • an area having a texture with a large luminance value distribution can be extracted as one area by discretization.
  • an area having a texture with a large luminance value distribution can be extracted as one area by discretization.
  • by connecting points having the same value and extracting a region having a shape close to a rectangle it becomes easy to remove regions other than buildings.
  • a region of a more rectangular shape as a building region among a plurality of region groups extracted for each of a plurality of different discretization widths, a locally optimal spatial smoothing parameter can be obtained. The function equivalent to what is applied can be exhibited.
  • (A) is an aerial photograph of Kiyomizu-dera, Hokkan-ji, and Kodai-ji around Higashiyama-ku, Kyoto. This area is a mixture of flat and dense hills.
  • (B) is the figure which showed the altitude distribution (DSM including the uppermost part of a building etc.) which the LiDAR data by an aircraft shows by color-coding for convenience (however, the color is not expressed in this figure).
  • (A) is the ground surface data extracted by the first process
  • (b) is the ground surface data extracted by the second process.
  • (A) shows DTM by the 1st process
  • (b) shows final DTM after performing the 2nd process.
  • (A) is an aerial photograph of the area around Gojo-dori, Nakagyo-ku, Kyoto. This area has a river in a flat area where buildings are densely packed.
  • (B) is the figure which showed the altitude distribution (DSM including the uppermost part of a building etc.) which the LiDAR data by an aircraft shows by color-coding for convenience (however, the color is not expressed in this figure).
  • (A) is the ground surface data extracted by the first process
  • (b) is the ground surface data extracted by the second process.
  • (A) shows DTM by the 1st process
  • (b) shows final DTM after performing the 2nd process.
  • (A) is an aerial photograph around Fushimi Momoyama, Fushimi-ku, Kyoto.
  • This area has a river with a large bank in a flat area where buildings are densely packed.
  • B is the figure which showed the altitude distribution (DSM including the uppermost part of a building etc.) which the LiDAR data by an aircraft shows by color-coding for convenience (however, the color is not expressed in this figure).
  • DTM altitude distribution
  • c shows the DTM by the first process
  • d shows the final DTM after the second process.
  • A) is the figure which expanded a part of LiDAR data shown to (b) of FIG. 7,
  • (b) is the ground photograph of the said area seen from the arrow direction in (a).
  • (c) shows a DTM (enlarged view) by the first process
  • (d) shows a final DTM (enlarged view) after the second process.
  • (A) shows the altitude distribution (DSM including the top of buildings, etc.) indicated by the LiDAR data by aircraft in the parking lot near Kiyomizu-dera, Higashiyama-ku, Kyoto, and its surroundings for convenience. (It is not expressed.) It is the figure shown, (b) is an aerial photograph. (A) shows DTM by the 1st process, (b) shows final DTM after performing the 2nd process. (A) is an aerial photograph of an urban area (75 m ⁇ 75 m), and (b) is a result of area division. (A), (b) is a three-dimensional building wire frame model, (a) has shown the planar shape, (b) has each shown the three-dimensional shape.
  • (A) is an aerial photograph of an urban area including a high-rise building (75 m ⁇ 75 m), and (b) is a result of area division.
  • (A), (b) is a three-dimensional building wire frame model, (a) has shown the planar shape, (b) has each shown the three-dimensional shape.
  • (A) is an aerial photograph of an urban area where a tree is adjacent to a building (75 m ⁇ 75 m), and (b) is a result of area division.
  • (A), (b) is a three-dimensional building wire frame model, (a) has shown the planar shape, (b) has each shown the three-dimensional shape.
  • (A) is an aerial photograph of an urban area with many dormitory roofs (75 m ⁇ 75 m), and (b) is a result of area division.
  • (A), (b) is a three-dimensional building wire frame model, (a) has shown the planar shape, (b) has each shown the three-dimensional shape.
  • It is a block diagram which shows an example of a structure of a DTM estimation apparatus. It is a flowchart which shows a DTM estimation method or an estimation program.
  • (A) is a figure which shows the normal vector of a gable roof
  • (b) is a figure which shows the normal vector of a dormitory roof.
  • FIG. 3 It is a figure which shows the result of the edge extraction in the area (area 3) where Takagi is adjacent to a building. It is a figure which shows the result of the edge extraction in the area (area 4) where many buildings with a dormitory roof exist. It is a figure which shows the result of the area
  • FIG. 2 It is a figure which shows the result of the area
  • FIG. It is a figure which shows the result of the area
  • FIG. It is a figure which shows the result of the area
  • FIG. It is a figure which shows the comparison result of the area
  • FIG. 1 It is a figure which shows the comparison result of the area
  • FIG. It is a figure which shows the comparison result of the area
  • FIG. It is a figure which shows the comparison result of the area
  • FIG. 1 It is a figure which shows the comparison result of the area
  • FIG. 19 is a block diagram showing an example of the configuration of a device that estimates DTM, which is altitude data only on the ground surface, for a predetermined range based on laser scanner data on the ground surface by an aircraft.
  • This DTM estimation apparatus includes a CPU 1, a memory 3 connected to the CPU 1 via a bus 2, an auxiliary storage device 4 such as a hard disk, interfaces 5 and 6, and a drive 7 and a display 8 connected to the interfaces 5 and 6, respectively.
  • this is a personal computer.
  • a storage medium 9 such as a CD or DVD containing LiDAR data (laser scanner data) to the drive 7 as an input device, the LiDAR data can be read.
  • the CPU1 is provided with the river extraction part 1a and the DTM estimation part 1b if it divides roughly as an internal function implement
  • the river extraction unit has a function of extracting a river region by connecting pixels in which no data exists in a unit grid in a predetermined range of LiDAR data (laser scanner data).
  • the DTM estimation unit 1b sets a first maximum allowable inclination value for data excluding the river region, estimates a provisional DTM, calculates a local inclination from the estimated DTM, and the inclination is a predetermined value. In the case of exceeding, the second maximum allowable inclination value larger than the first maximum allowable inclination value is set and the DTM is estimated again.
  • the estimated DTM can be displayed on the display 8 as a display unit.
  • the program is stored in the auxiliary storage device 4 (FIG. 19), which is a storage medium, and is executed by the CPU 1 using the memory 3.
  • the DTM estimation program is installed and executed from various storage media in the auxiliary storage device 4, or downloaded from a specific server storing the DTM estimation program via the network and installed in the auxiliary storage device 4. Executed. That is, the DTM estimation program can exist and be distributed as a storage (recording) medium.
  • the LiDAR data to be used is a commercially available product, and includes, for example, only last pulse data (reception data finally reflected on the ground surface) having an average point density of 1 point / 1 m 2 and three-dimensional coordinate data.
  • usable data is not limited to this.
  • continuous data from the first to the last called the full wave type can be used.
  • the CPU 1 extracts a river region from LiDAR data (step S1). Specifically, in an image set to a 1 m grid as a unit grid, pixels that do not have data (points) are connected, and an area occupying a certain area or more is extracted as a river area.
  • the CPU 1 extracts the local minimum value from the LiDAR data (step S2). Specifically, the local minimum value is extracted within a 50 m ⁇ 50 m window at the first time, and the local minimum value is searched within a 25 m ⁇ 25 m window at the second time.
  • the slope calculated from two points on the three-dimensional coordinates that are the first local minimum and the second local minimum are not significantly higher than the first local minimum (for example, within 50 cm). Is less than or equal to the first maximum allowable inclination value (for example, 3 degrees), the first local minimum value is updated to the second local minimum value. In other cases, the first local minimum value is adopted.
  • the CPU 1 estimates a plane from the LiDAR data by the least square method (step S3), and when the RMSE (Root ⁇ ⁇ ⁇ ⁇ Mean Square Errors) for the plane is small, those points are on a common plane. It shall be. Specifically, for example, a plane equation is calculated using all point groups included in a 4 m ⁇ 4 m region, and if the RMSE for a plane is 10 cm or less, the plane equations are assigned to all point groups. .
  • This plane is a candidate road surface, and may be a slope, but an extreme slope has a very low possibility of a road, so the vertical component of the normal of the plane is limited to 0.9 or more, for example.
  • the CPU 1 sets an initial value (the initial point group) of the ground surface data as a point having a small elevation difference from the local minimum value among the points on the plane and within the inclination viewed from the local minimum value. ) Is selected (step S4). Subsequently, the CPU 1 searches for a nearby point of the ground surface data, calculates the distance to the plane using the plane equation of the ground surface data, and if the distance is equal to or less than a threshold (for example, 10 cm), Candidates for additional ground surface data. If the slope formed with the ground surface data is equal to or less than the maximum allowable slope value (for example, 3 degrees), the ground surface data is added (step S5).
  • a threshold for example, 10 cm
  • step S5 The addition of the ground surface data (step S5) is performed until no points are added (step S6). When there are no more points to be added, the extraction of the ground surface data is complete. Subsequently, the CPU 1 searches for nearby ground surface data for a point where the ground surface data does not exist, and interpolates to estimate the DTM (step S7). If the river area is encountered during this search, the CPU 1 stops searching in that direction.
  • the CPU 1 calculates a local gradient based on the estimated DTM (step S8).
  • the average value of the ground surface height is calculated by using a 21 m ⁇ 21 m area as a unit. Search for the surrounding area (up / down / left / right / oblique) with respect to the ground surface height of a certain area, and calculate the slope from the difference between the distance between the center of the area and the average value of the ground surface height. Find the maximum slope value.
  • a threshold value for example, 4.5 degrees
  • the maximum allowable inclination value is set larger and updated to the second maximum allowable inclination value (for example, updated to 4.5 degrees).
  • the CPU 1 determines whether or not the processing of steps S5 to S8 is the first time (step S9).
  • the first is “Yes”, and the CPU 1 executes steps S5 and S6 again based on the second maximum allowable inclination value.
  • the second maximum allowable inclination value larger than the first maximum allowable inclination value, a point having a locally steep inclination is added as the ground surface data. That is, in the process of the first step S5, point cloud data such as river banks and hills that have been leaked from the additional data are extracted and added as ground surface data. As a result, the DTM is estimated with higher accuracy than the first time.
  • step S9 for the second time “No” is given, and the process is terminated.
  • FIG. 1 (a) is an aerial photograph of the area around Kiyomizu-dera, Hokanji, and Kodai-ji, Higashiyama-ku, Kyoto. This area is a mixture of flat lands with dense buildings and hilly terrain.
  • FIG. 1B is a diagram showing the altitude distribution (DSM including the uppermost part of a building or the like) indicated by the LiDAR data obtained by the aircraft by color-coding for convenience (however, the color is not expressed in this figure). It is.
  • the black part on the left is a part with a relatively low altitude
  • the black part on the right is a part with a relatively high altitude.
  • the intermediate whitish part is their intermediate elevation. Since DSM includes buildings, the ground surface is difficult to understand.
  • FIG. 2 is the ground surface data extracted in the first process
  • (b) is the ground surface data extracted in the second process.
  • the black part in the figure represents the ground surface (ground surface).
  • 3A shows a DTM by the first process
  • FIG. 3B shows a final DTM after the second process. 2 and 3, as is clear from the comparison between (a) and (b), the accuracy of the detailed data on the hills (right side) and flat land is greatly improved by performing the second processing. It is well obtained.
  • FIG. 4 is an aerial photograph around Gojo-dori, Nakagyo-ku, Kyoto. In this area, there are rivers in a flat area where buildings are densely packed.
  • FIG. 4B is a diagram showing the altitude distribution (DSM including the uppermost part of the building or the like) indicated by the LiDAR data by the aircraft in different colors (however, the color is not expressed in this figure). It is.
  • FIG. 5 is the ground surface data extracted in the first process
  • (b) is the ground surface data extracted in the second process
  • 6A shows the DTM by the first process
  • FIG. 6B shows the final DTM after the second process.
  • there is no dramatic difference between the first time and the second time In other words, in a flat terrain having no locally steep slope, the same first maximum allowable slope value as that in the first time is used for the second processing, so that excessive extraction is not performed.
  • FIG. 7A is an aerial photograph around Fushimi Momoyama in Fushimi-ku, Kyoto. In this area, there is a river with a large bank in a flat area with dense buildings.
  • FIG. 7B is a diagram showing the altitude distribution (DSM including the uppermost part of the building or the like) indicated by the LiDAR data by the aircraft in different colors (however, the color is not expressed in this figure). It is.
  • the ellipse in (a) is a highway and was not yet constructed at the time of LiDAR data acquisition of (b).
  • (c) of FIG. 7 shows DTM by the 1st process
  • (d) shows the last DTM after performing the 2nd process.
  • FIG. 8A is an enlarged view of a part of the LiDAR data shown in FIG. 7B, and FIG. 8B is a ground photograph of the area as seen from the arrow direction in FIG. It is.
  • FIG. 8C shows a DTM (enlarged view) by the first process
  • FIG. 8D shows a final DTM (enlarged view) after performing the second process.
  • the white line extending vertically in the center is a river, and the black parts on the left and right are banks. From the comparison between (c) and (d), it can be seen that the bank that could not be handled as the ground surface at the first time was recognized with high accuracy as the ground surface the second time.
  • FIG. 9A shows the altitude distribution (DSM including the top of the building, etc.) indicated by the LiDAR data by aircraft in the parking lot near Kiyomizu-dera, Higashiyama-ku, Kyoto, and its surroundings for convenience. (The color is not represented.)
  • B is an aerial photograph. 10A shows the DTM by the first process, and FIG. 10B shows the final DTM after the second process. From the comparison between (a) and (b), it can be seen that the parking lot (the black portion in the center) that could not be captured as the ground surface at the first time was accurately recognized as the ground surface at the second time.
  • Example 5 Although not shown, the first and second DTMs were similarly obtained for the rooftop parking lot of the building, but there was almost no change. In this case, the rooftop parking lot was not mistakenly recognized as the ground surface.
  • the river area is excluded in the estimation of the DTM, so that the accuracy of the DTM can be prevented from being deteriorated due to the search for the ground surface across the river.
  • a second maximum allowable slope value that is larger than the first maximum allowable slope value is set and the DTM is estimated again, so that the first time The surface data that could not be searched can be added to improve the accuracy of DTM estimation.
  • the second maximum allowable slope value may be increased or decreased according to the features of the topography and the townscape.
  • Kyoto City it was confirmed that 4.5 degrees is suitable. It is assumed that 4.5 degrees or more is suitable for land with more slope than Kyoto City. Conversely, it is assumed that less than 4.5 degrees (greater than 3 degrees) is suitable for land with less slope than Kyoto City.
  • region division As an outline of region division, first, a region having a texture having a large variance of luminance values is also extracted. Therefore, region division is attempted after discretizing luminance values into a small number of values. Plural widths of luminance values to be discretized are prepared to cope with different degrees of dispersion. In order to preferentially extract a region close to a rectangle, an index called a rectangle index is calculated. Considering the possibility of being divided by shadows, if the rectangle index is improved in combination with adjacent areas, both are integrated.
  • the luminance value of the original image is “0 to 39” and the same discretized value is given, and similarly “40 to 79”, “80 to 119”, “120 to 159”, “160 to” 199 ”,“ 200 to 239 ”, and“ 240 to 255 ”are discretized into seven sections.
  • the following parameters were used in this experiment.
  • a label number is assigned to a region closed by an edge, and a region that seems to be vegetated is removed based on the RGB luminance value.
  • a large area (6400 pixels) with a certain area or more is removed, and a small area with a certain area or less (30 pixels in this experiment) is integrated if there is an area with a certain area in the periphery, or removed if there is no area.
  • An index called a rectangular index for each region is calculated as follows.
  • the first axis and the second axis are determined from the two-dimensional coordinates of a set of edge of the region (referred to as edge group).
  • edge group The range of the area is expressed by the values of the first axis and the second axis, and (maximum value-minimum value + 1) of the first axis is multiplied by (maximum value-minimum value + 1) of the second axis. , Get the rectangular area surrounding the region.
  • C The area of the actual area / the area of the rectangle surrounding the area is defined as the rectangle index.
  • D If the rectangular index falls below a certain value (0.4), it is judged that there is a high possibility that it is not a building and is excluded.
  • the area obtained by the area division is not a roof unit, and there is a case where a pair of sloped roofs in the same building or a roof of another building is not separated.
  • the roof / building may not be extracted, and conversely, trees and roads may be extracted.
  • the following measures are taken to reduce modeling leakage.
  • the estimated DTM is removed from the LiDAR data. Most of the trees can be deleted by setting the RMSE tolerance when calculating the roof plane and normal vectors. If the roofs facing each other in the gable roof building are not separated, the distribution of normal vectors is examined and if there is a high possibility that they are mixed, the gable roof is generated. In the situation where the roof on one side is not separated in the adjacent gable roof building (the number of roofs is 1 to 2 across the ridge), two sets of two are separated by separating the unseparated area into two roofs Create a one-to-one roof pair and create two gable roof buildings.
  • a filtering process is performed on LiDAR data, which is point cloud data of three-dimensional coordinates, using a DTM to separate the ground surface and the non-ground surface.
  • the normal vector in the area of the roof or building is calculated from the non-surface data. If there are 6 or more points in a 9 m ⁇ 9 m area, a plane equation is calculated. If RMSE is 10 cm or less, normal vectors are assigned to all point groups used in the calculation. If the normal of the region is inclined, a paired region is searched.
  • ⁇ For non-ground surface data, classify each point group whose horizontal distance and vertical distance from the nearest point data are within a certain value. Of these, the main azimuth of the point cloud data is determined using only the point cloud corresponding to the area of the roof or building.
  • a wire frame model is created in the order of dormitory roof, gable roof, and other roofs. When creating each model, care should be taken that the ridges and contours are the same or perpendicular to the main azimuth of the point cloud data.
  • the dormitory roof is generated when there are two pairs of roofs to be paired.
  • the gable roof there is a region where a pair of roofs exist, and the gable roof is generated in both regions.
  • a region in which the normal vector distribution indicates a mixed state of a plurality of normals is divided into two in the region, and a gable roof is generated when the normal vectors are clearly separated. If not, create a flat roof.
  • FIG. 21 shows the normal vector of the gable roof
  • (b) shows the normal vector of the dormitory roof.
  • a flat roof is generated. A point that is not extracted as a building but has a certain height from the ground surface and is on a specific plane is selected and included in the final result.
  • FIGS. Example 1
  • a of FIG. 11 is an aerial photograph of a city area (75 m ⁇ 75 m).
  • B is a result of area division. For convenience, they are color-coded (colors are not displayed in this figure), but the colors themselves have no meaning.
  • the outline of the roof is generally captured. There are also places where the characteristics of gables and dormitories are captured.
  • 12A and 12B are three-dimensional building wire frame models, where FIG. 12A shows a planar shape and FIG. 12B shows a three-dimensional shape. It can be seen that the roof shapes such as flat, gable, and dormitory are clearly visible.
  • FIG. 13A is an aerial photograph of a city area including a high-rise building (75 m ⁇ 75 m).
  • B is a result of area division. For convenience, they are color-coded (colors are not displayed in this figure), but the colors themselves have no meaning.
  • the outline of the roof is generally captured. There are also places where the characteristics of gables and dormitories are captured.
  • 14A and 14B are three-dimensional building wire frame models, in which FIG. 14A shows a planar shape and FIG. 14B shows a three-dimensional shape. It can be seen that the roof shapes such as flat, gable, and dormitory are clearly visible.
  • FIG. 15 is an aerial photograph of the urban area where Takagi is adjacent to the building (75 m ⁇ 75 m).
  • (B) is a result of area division. For convenience, they are color-coded (colors are not displayed in this figure), but the colors themselves have no meaning.
  • the outline of the roof is generally captured. There are also places where the characteristics of gables and dormitories are captured.
  • 16A and 16B are three-dimensional building wire frame models, where FIG. 16A shows a planar shape and FIG. 16B shows a three-dimensional shape. It can be seen that the roof shapes such as flat, gable, and dormitory are clearly visible.
  • FIG. 17 is an aerial photograph of an urban area with many dormitory roofs (75 m ⁇ 75 m).
  • B is a result of area division. For convenience, they are color-coded (colors are not displayed in this figure), but the colors themselves have no meaning.
  • the outline of the roof is generally captured. There are also places where the characteristics of gables and dormitories are captured.
  • 18A and 18B are three-dimensional building wire frame models, where FIG. 18A shows a planar shape and FIG. 18B shows a three-dimensional shape. It can be seen that the roof shapes such as flat, gable, and dormitory are clearly visible.
  • a three-dimensional building model can be created by estimating the shape of a roof such as a building. Further, by extracting a region for creating a three-dimensional building model in advance, the accuracy of identifying the roof region can be improved, and a more accurate three-dimensional building model can be created.
  • FIG. 23 is a photograph of the area taken on the ground. Further, as the aerial photograph data of this area, data (commercially available product) having been subjected to ortho (orthogonal projection) processing at 25 cm resolution was used. This data is aligned with the planar rectangular coordinate system. In addition, even if it has been ortho-processed, the side of the building can be seen slightly, but there is no particular problem in practical use.
  • FIG. 24 is a block diagram showing an example of the computer device 10 for performing the region extraction method, and the region extraction program causes the computer device 10 to realize a function.
  • the computer apparatus 10 having the function of the area extraction program is an area extraction apparatus according to an embodiment of the present invention.
  • the computer 10 includes a CPU 1, a memory 3 connected to the CPU 1 via a bus 2, an auxiliary storage device 4 such as a hard disk, interfaces 5, 6 and drives 7 connected to the interfaces 5 and 6, respectively. And a display 8, typically this is a personal computer.
  • a storage medium 9 such as a CD or DVD containing aerial photographs to the drive 7 as an input device, the aerial photograph data can be read.
  • the area extraction program is installed and executed from various storage media in the auxiliary storage device 4 separately from the aerial photograph, or downloaded from a specific server storing the area extraction program via the network, It is installed in the storage device 4 and executed. That is, the area extraction program can exist and be distributed as a storage (recording) medium.
  • region extraction a main step of the method / program
  • region extraction in order to effectively extract a roof having a texture with a large variance of luminance values, the luminance values are discretized into a predetermined number of values, and regions having the same discretized values are labeled. And the area
  • the spatially smoothing parameters that were optimal locally were applied by sequentially adopting regions that were closer to a rectangle. The same processing is realized.
  • each step will be described in detail.
  • N disc types of discretization widths are set for luminance values of an arbitrary 1 band of 1 byte (luminance value range: 0 to 255) ⁇ 3-band (RGB) image.
  • FIG. 25 shows an example of searching for the edge group 1a in the upper left direction and the edge group 1b in the lower right direction. For example, (a) group 1a has d1 or more edge groups, and group 1b has d1 or more edge groups, (b) group 2a has d2 or less edge groups, and group 2b has only d2 or less edge groups. If both conditions (a) and (b) are satisfied, the non-edge center pixel is complemented as an edge.
  • the condition (b) prevents the pixels inside the area near the corners of the rectangular area from being erroneously extracted as edges.
  • a label number is assigned (labeled) to the region closed by the edge, and the region that seems to be vegetation is removed based on the RGB luminance value.
  • a large area for example, 6400 pixels or more
  • a small area having a certain area or less for example, 30 pixels or less
  • B red , B grn , and B blue represent the luminance values of the red, green, and blue bands
  • T veg and DN represent the luminance value threshold values
  • T veg, g2b_ratio , and T veg and r2b_ratio represent the index threshold values . If any of R grn_veg and R red_veg is greater than or equal to the threshold value T veg, ratio , it is removed as a vegetation region.
  • an index called a rectangular index for each region is calculated as follows.
  • A The first axis and the second axis are determined from the two-dimensional coordinates of a set of edges in the region (referred to as an edge group).
  • B Express the existence range of the area with the value of the first axis and the second axis, and multiply by (maximum value-minimum value + 1) of the first axis and (maximum value-minimum value + 1) of the second axis. , Get the rectangular area surrounding the region.
  • C (Area area / Area area surrounding the area) is defined as a rectangle index. Here, if the rectangular index falls below a certain value (for example, 0.4), it is determined that there is a high possibility that it is not a building, and is excluded from the extraction target.
  • FIG. 26 shows a conceptual diagram of rectangular index calculation.
  • the first and second axes are calculated from the edge group of a certain region, and the smallest rectangle as shown in the figure surrounding the target region is obtained from the rectangles whose sides are parallel to the first and second axes.
  • the slope of the straight line is voted by using an edge group in which the distance between two points is a constant value (between 5 pixels and 20 pixels), and the slope at which the maximum frequency occurs is determined.
  • the orientation was the first axis, and the second axis was determined to be orthogonal to the first axis.
  • i dx S actual / S rect (1) It is defined as Here, i dx represents a rectangular index, S actual represents the area of a certain region, and S rect represents the area of a rectangle surrounding the region.
  • the rectangular index takes a value in the range of 0 to 1, and the closer to 1, the closer to the rectangle. With such an index, the degree of “closeness” in a shape close to a rectangle can be expressed simply and accurately as an index.
  • step 7 using the rectangular index of equation (1), the roof / building divided by the shadow is integrated to estimate the original area. For example, search for a region near a certain region A, (A) The rectangular index when merged is improved over each rectangular index. (B) The specified threshold value (0.7) or more. (C) The difference in the average luminance value of the region is within a certain value (30). It is determined whether or not the condition is satisfied. When satisfy
  • a rectangular index when the regions ⁇ , ⁇ , and ⁇ adjacent to the target region ⁇ are integrated is calculated.
  • the first and second axes are calculated using an edge group excluding the edge groups close to each other among the edge groups of the two areas.
  • the axis is determined by voting the slope of a straight line using an edge group in which the distance between two points is a constant value (5 pixels or more and 20 pixels or less), and the slope with the highest frequency is taken as the direction of the first axis.
  • the second axis was determined to be orthogonal to the first axis.
  • FIG. 28 is a diagram illustrating a result of edge extraction in an urban area where low-rise buildings are densely gathered (hereinafter, referred to as area 1).
  • A is an aerial photograph in which an area having an actual size of about 75 m ⁇ 75 m is photographed with a photographing resolution of 300 ⁇ 300 pixels.
  • B is a result of using a known Canny filter for the data of (a).
  • C is an edge extraction result (first stage) when the discretization width 40 is used in the discretization described above. Further, (d) is an edge extraction result (second stage) when the discretization width 30 is used, and (e) is an edge extraction result (third stage) when the discretization width 20 is used.
  • FIG. 29 is a diagram showing a result of edge extraction in an urban area where high-rise buildings and low-rise buildings are mixed (hereinafter referred to as region 2).
  • region 2 is an aerial photograph
  • region 3 is a result of using a canny filter
  • region 4 is an edge extraction result when using a discretization width 40.
  • First stage is an edge extraction result (second stage) when the discretization width 30 is used
  • e is an edge extraction result (third stage) when the discretization width 20 is used. is there.
  • FIG. 30 is a diagram showing the result of edge extraction in an area where Takagi is adjacent to a building (hereinafter referred to as area 3).
  • area 3 is a diagram showing the result of edge extraction in an area where Takagi is adjacent to a building (hereinafter referred to as area 3).
  • (A) to (e) are the same as those in FIG. 28,
  • (a) is an aerial photograph
  • (b) is a result of using a canny filter
  • (c) is an edge extraction result when using a discretization width 40.
  • First stage is an edge extraction result (second stage) when the discretization width 30 is used, and
  • e) is an edge extraction result (third stage) when the discretization width 20 is used. is there.
  • FIG. 31 is a diagram showing the result of edge extraction in an area where there are many buildings with a dormitory roof (hereinafter referred to as area 4).
  • area 4 a diagram showing the result of edge extraction in an area where there are many buildings with a dormitory roof (hereinafter referred to as area 4).
  • (A) to (e) are the same as those in FIG. 28,
  • (a) is an aerial photograph
  • (b) is a result of using a canny filter
  • (c) is an edge extraction result when using a discretization width 40.
  • First stage (d) is an edge extraction result (second stage) when the discretization width 30 is used, and (e) is an edge extraction result (third stage) when the discretization width 20 is used. is there.
  • FIG. 32, 33, 34, and 35 show the results of region extraction in the above-described region 1 (FIG. 7), region 2 (FIG. 29), region 3 (FIG. 30), and region 4 (FIG. 31), respectively.
  • FIG. 32 to 35 (a) is an aerial photograph of each region, (b) is a region extraction result when using a discretization width 40 (first stage), and (c) is a discretization width 30. Region extraction result when used (second stage), (d) is an edge extraction result when using the discretization width 20 (third stage), and (e) is a result of three results (( It is a figure of the final result which integrated b)-(d)). In any of FIGS. 32 to 35, (e) is in a state where the best of the three results are locally taken, as is clear when compared with the results of (b) to (d). Regions (edges) are extracted more accurately.
  • FIG. 36 is a diagram for region 1.
  • the edge of the building is clearly displayed in the color image, but in this figure, the edge looks like the edge of the building.
  • FIG. 37 is a diagram for region 2.
  • FIG. 38 is a diagram for region 3.
  • FIG. 39 is a diagram for the area 4.
  • edges obtained by applying different discretization widths are integrated. This is essentially the same as conversion to different spatial resolutions and processing. However, the edge is preserve
  • the optimal spatial scale parameter can be selected locally according to the size of the roof or building through the integration process.
  • labeling using the discretized image takes time, and the processing time increases as the image size increases. For example, in the case of a computer using a 6 GB memory for a CPU operating at a clock of 3.2 GHz, the average processing time in each target area was about 90 seconds.
  • the triangular area seen on the dormitory roof is also well extracted.
  • the ideal triangle rectangle index remains at 0.5 and is not preferentially adopted.
  • a factor that has been successfully extracted this time is that the rectangular area around the triangular area has been successfully extracted.
  • the triangular area and the rectangular area may be integrated.
  • the rectangle index of the area where the triangle area and the rectangle area are integrated is lower than the rectangle index of the original rectangle alone area. It is hard to be adopted. Based on such a selection process, the rectangular area was first selected, and then the triangular area was extracted.
  • the result of the above method has a wider range where no region exists. This is partly due to the upper limit on the area of the area, but it is largely due to the fact that if the rectangular index of the area is less than a certain value (0.4), the possibility of roofs / buildings is low. is doing. Most but not all of the vegetation and the area including the cemetery in FIG. 36 are far from being rectangular in shape and can be removed as a result. In other words, it can be said that it functions sufficiently in terms of building area extraction.
  • the rectangular index will be described from the viewpoint of area.
  • the rectangle index is extracted sequentially, a relatively small area is often extracted preferentially, and an area that captures a rough outline is selected even if the rectangle index is slightly small. It may not be.
  • we examined the method. For example, by executing the correction using the following equation (2), a region corresponding to a plurality of roofs / buildings is obtained in a part of the region which is larger than the results shown in FIGS. ⁇ i dx C 1 ⁇ ln (S actual ⁇ C 2 + C 3 ) (2)
  • ⁇ i dx represents a correction value of the rectangular index
  • C 1 to C 3 represent empirically determined coefficients.
  • the number of buildings extracted decreased. This time, we are targeting the area where traditional Japanese buildings are lined up, but the functions and parameter values in this rectangular index correction should be adjusted with the building to be extracted in mind.
  • FIGS. 36 to 39 (c) show the result of region extraction when the axis is determined by principal component analysis.
  • the slate roof region was divided and partly extracted. If the principal component analysis is applied using the edge group unless it is a complete rectangle, the first principal component is not always parallel to the roof edge, especially in the roof that is divided by the shadow as the object of this time. No results were obtained.
  • 1-c four roofs are extracted as one large region in the region extraction result using principal component analysis.
  • the region extraction result lacks stability by using principal component analysis in this way, the principal component analysis is not applied, and the slope of the straight line using an edge group in which the distance between two points is a constant value.
  • the slope at which the maximum frequency occurs was defined as the direction of the first axis, and the second axis was determined to be orthogonal to the first axis.
  • This threshold value should be determined empirically according to the characteristics of the target area.
  • the purpose is to extract buildings in urban areas, and vegetation was a target to be removed.
  • a method of removing pixels with high vegetation by referring to luminance values in pixel units is also conceivable.
  • the original shape of the roof or building may not be obtained, and it may be divided or the area may be too small to be extracted.
  • a building or a roof can be effectively extracted as a region in a densely built-up area.
  • the luminance values are discretized into a small number of values, and regions having the same discretized values are labeled.
  • the rectangular index calculated from the edge of the area is used to preferentially extract the rectangular area seen on the roof or building.
  • a process of closing the edge when the possibility of a rectangular or triangular shape is high was added.
  • the triangular area seen on the dormitory roof has a low rectangular index and is not preferentially adopted, but as a result, it was successfully extracted.
  • the rectangular index of the area where the triangular area and the rectangular area are integrated is less likely to be adopted than the rectangular index of the original single rectangle area, and the triangular area is extracted by selecting the rectangular area first.
  • the aerial photograph is used as the original data, but it is also possible to use data of a photograph taken from an artificial satellite photographed with high accuracy instead of the aerial photograph.
  • the embodiment disclosed this time is illustrative and not restrictive in all respects.
  • the scope of the present invention is defined by the terms of the claims, and is intended to include any modifications within the scope and meaning equivalent to the terms of the claims.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Software Systems (AREA)
  • Computer Graphics (AREA)
  • Remote Sensing (AREA)
  • Computer Hardware Design (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Optics & Photonics (AREA)
  • Structural Engineering (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Civil Engineering (AREA)
  • Architecture (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Image Processing (AREA)
  • Processing Or Creating Images (AREA)

Abstract

　ＤＴＭ推定に関しては、所定範囲における単位グリッド内にデータが存在しない画素を連結して河川領域を抽出し、河川領域を除くデータについて、第１の最大許容傾斜値（例えば３度）を設定して暫定的なＤＴＭを推定し、推定したＤＴＭから局所的な傾斜を計算し、傾斜が所定値を超える場合は、第１の最大許容傾斜値より大きい第２の最大許容傾斜値（例えば４．５度）を設定して再度ＤＴＭを推定する。　領域抽出に関しては、航空写真等のデータに基づいて建物の領域を抽出するに際し、複数の異なる離散化幅を設定し、それぞれの離散化幅について、データの輝度値を当該離散化幅で離散設定された複数の値に離散化する。また、離散化して得られた離散化画像において同一値を持つ画素を連結し、建物の領域の候補として長方形に近い形状の領域を抽出する。そして、複数の異なる離散化幅ごとに抽出された複数の領域群のうち、より長方形に近い形状の領域を、建物の領域として採用する。

Description

ＤＴＭ推定方法、ＤＴＭ推定プログラム、ＤＴＭ推定装置及び３次元建物モデルの作成方法、並びに、領域抽出方法、領域抽出プログラム及び領域抽出装置
　［ＤＴＭ推定方法、ＤＴＭ推定プログラム、ＤＴＭ推定装置及び３次元建物モデルの作成方法について］
　本発明は、航空機による地表のレーザースキャナーデータに基づいてＤＴＭ（Digital Terrain Model）を推定する技術に関する。
　［領域抽出方法、領域抽出プログラム及び領域抽出装置について］
　また、本発明は、航空機又は人工衛星から撮影した写真のデータに基づいて建物の領域を抽出する方法、プログラム及び装置に関する。
　［ＤＴＭ推定方法、ＤＴＭ推定プログラム、ＤＴＭ推定装置及び３次元建物モデルの作成方法について］
　航空機から地表をレーザースキャンした反射波として得られるＬｉＤＡＲ（Light Detection And Ranging）データは、高精度な三次元座標の点群データとして、建物の高さや森林の樹高を測定することに利用可能である（例えば、特許文献１参照。）。このＬｉＤＡＲデータにより、所望の領域について、地表面（敷地や道路等の地表面）の他、建物の最上部や樹木の最上端を含むＤＳＭ（Digital Surface Model）を得ることができる。
　一方、部分的に得られた地表面のデータ及び、これらから全体を補間したデータによって、地表面のみを表すＤＴＭを推定することができる。これによって、例えば、建物の高さは、（ＤＳＭ－ＤＴＭ）として算出することができる（例えば、非特許文献１参照。）。ＤＴＭを推定するには、ＬｉＤＡＲデータの所定範囲について、その中の標高最低点が地表面であるとする。そして、当該範囲内で地表面を探索し、かつ、内挿（補間）すれば、ＤＴＭを推定することができる。
　［領域抽出方法、領域抽出プログラム及び領域抽出装置について］
　都市を俯瞰した場合の特徴抽出、領域抽出は昔から重要なテーマであった（例えば、非特許文献１）。例えば、航空写真又は人工衛星画像からの道路抽出は盛んに研究されてきた（例えば、非特許文献３，４）。リモートセンシング画像の分類手法は、画素単位の分類と領域（オブジェクト）単位の分類とに大別できる。クラスタリング（いわゆる教師なし）、最尤法（教師つき）、Support Vector Machines (SVM)（教師つき）に代表される画素単位の分類手法は、統計的理論に基づき、各画素を分類クラスに割り当てていく（例えば、非特許文献５）。一方、領域単位の分類手法では、画素の周辺の情報（Context）を利用する。例えば、数学的 Morphological アプローチを用いた分類などが代表例として挙げられる（例えば、非特許文献６～１０）。
　図２２の（ａ）及び（ｂ）は、家屋が密集した市街地（京都市東山区）の航空写真の一例である。このような密集市街地では、従来の領域抽出手法が効果的に機能しない。その要因を分析すると、まず、（イ）建物同士が接近しているため、屋根の輝度値、テクスチャが似ていると、建物の境界線が明瞭でなくなる点である。また、（ロ）隣接した建物の影が写り込みやすい、という密集市街地特有の要因がある。さらに、日本の伝統的建物では波型の横断面形状を持つ屋根（例えばスレート屋根）が使用されていることもあり、（ハ）輝度値にばらつきのあるテクスチャの屋根の輪郭抽出が困難である、という対象地域の建物特性による要因も確認できる。
特許第４０５８２９３号
「航空機ＬｉＤＡＲデータからの密集市街地における地表面データのフィルタリングアルゴリムの構築」、須▲崎▼純一、幸良淳志、児島利治、応用測量論文集Ｖｏｌ．２１、第７８～８９頁 Weng, Q. and Quattrochi, D.A., "Urban Remote Sensing", CRC Press, 2007 Hu, J., Razdan, A., Femiani, J.C., Cui, M. and Wonka, P., "Road network extraction and intersection detection from aerial images by tracking road footprints", IEEE Transactions on Geoscience and Remote Sensing, Vol. 45, No. 12, pp. 4144 - 4157, 2007 Movaghati, S., Moghaddamjoo, A. and Tavakoli, A., "Road extraction from satellite images using particle filtering and extended Kalman filtering", IEEE Transactions on Geoscience and Remote Sensing, Vol. 48, No. 7, pp. 2807 - 2817, 2010 Tso, B. and Mather, P.M., "Classification methods for remotely sensed data -- 2nd ed.", CRC Press, 2009 Soille, P. and Pesaresi, M. , " Advances in mathematical morphology applied to geoscience and remote sensing", IEEE Transactions on Geoscience and Remote Sensing, Vol. 40, No. 9, pp. 2042-2055, 2002 Benediktsson, J.A., Pesaresi, M. and Amason, K., "Classification and feature extraction for remote sensing images from urban areas based on morphological transformations", IEEE Transactions on Geoscience and Remote Sensing, Vol. 41, No. 9, pp. 1940 - 1949, 2003 Benediktsson, J.A., Palmason, J.A. and Sveinsson, J.R., "Classification of hyperspectral data from urban areas based on extended morphological profiles", IEEE Transactions on Geoscience and Remote Sensing, Vol. 43, No. 3, pp. 480 - 491, 2005 Bellens, R., Gautama, S., Martinez-Fonte, L., Philips, W., Chan, J.C.-W., Canters, F., "Improved Classification of VHR Images of Urban Areas Using Directional Morphological Profiles", IEEE Transactions on Geoscience and Remote Sensing, Vol. 46, No. 10, pp. 2803 - 2813, 2008 Tuia, D., Pacifici, F., Kanevski, M. and Emery, W.J., "Classification of Very High Spatial Resolution Imagery Using Mathematical Morphology and Support Vector Machines", IEEE Transactions on Geoscience and Remote Sensing, Vol. 47, No. 11, pp. 3866 - 3879, 2009 Canny, J., "A Computational Approach to Edge Detection", IEEE Transactions on Pattern Analysis and Machine Intelligence, Vol. PAMI-8, No. 6, pp. 679 - 698, 1986 Zhong, J. and Ning, R., "Image denoising based on wavelets and multifractals for singularity detection", IEEE Transactions on Image Processing, Vol. 14, No. 10, pp. 1435 - 1447, 2005 Tonazzini, A., Bedini, L. and Salerno, E., "A Markov model for blind image separation by a mean-field EM algorithm", IEEE Transactions on Image Processing, Vol. 15, No. 2, pp. 473 - 482, 2006 Berlemont, S. and Olivo-Marin, J.-C., "Combining Local Filtering and Multiscale Analysis for Edge, Ridge, and Curvilinear Objects Detection", IEEE Transactions on Image Processing, Vol. 19, No. 1, pp. 74 - 84, 2010 Ding, J., Ma, R. and Chen, S., "A Scale-Based Connected Coherence Tree Algorithm for Image Segmentation", IEEE Transactions on Image Processing, Vol. 17, No. 2, pp. 204 - 216, 2008 Chien, S .Y., Ma, S.Y. and Chen, L.G., "Efficient moving object segmentation algorithm using background registration technique", IEEE Transactions on Circuits and Systems for Video Technology, Vol. 12 , No. 7, pp. 577 - 586, 2002 Tsai, V.J.D., "A comparative study on shadow compensation of color aerial images in invariant color models", IEEE Transactions on Geoscience and Remote Sensing, Vol. 44, No. 6, pp. 1661 - 1671, 2006 Ma, H., Qin, Q. and Shen, X., "Shadow Segmentation and Compensation in High Resolution Satellite Images", Proceedings of IEEE International Geoscience and Remote Sensing Symposium, 2008, Vol. 2, pp. II-1036 - II-1039, 2008
　［ＤＴＭ推定方法、ＤＴＭ推定プログラム、ＤＴＭ推定装置及び３次元建物モデルの作成方法について］
　上記のＤＴＭは、ＬｉＤＡＲデータと違って推定であるため、多少の誤差・誤認識が生じることは避けられないが、いかに当該誤差・誤認識を抑制するかが重要である。種々の地表面に対してＤＴＭ推定を行ってみた結果、特に、全体として平坦な土地に河川がある場合に、河床、堤防、橋等のデータの点群を、地表面として認識するのが難しいことが判明した。また、道路から徐々に上がっていく屋上駐車場（建物の最上部）が、誤って地表面と認識されることもあった。
　かかる従来の問題点に鑑み、本発明は、ＬｉＤＡＲデータに基づいて推定されるＤＴＭの精度を高めることを目的とする。
　［領域抽出方法、領域抽出プログラム及び領域抽出装置について］
　上記（イ）の要因は、言い換えれば、領域のエッジが十分に得られていない場合に、どう対処するか、という課題を与えている。また、上記（ハ）の要因により過剰に抽出されるエッジから不要なエッジを除外することも課題である。ここで、（ハ）の課題は（イ）の課題とも関連していると考えられる。例えば、非特許文献１１において、Cannyはノイズに強いエッジ抽出オペレータを提唱している。これは、これまでに幅広く使われている代表的なエッジ抽出オペレータの一つといえる。またWaveletを使ってノイズを除去したり（例えば、非特許文献１２）、平均値（例えば、非特許文献１３）や多重解像度画像（例えば、非特許文献１４）を活用したりする手法が報告されている。
　また上記（ロ）の要因に関連して、影領域の輝度値を回復する研究事例が報告されている（例えば、非特許文献１５～１８）。しかしながら例えば非特許文献１８における回復結果を見る限りでは、回復後の画像において元々の影領域と非影領域との境界線が明瞭に出てしまうという問題点、及び、影領域内の補正値の決定が難しく、過剰に回復することもある、といった問題点が、依然として解決できていない。
　かかる従来の問題点に鑑み、本発明は、航空写真等のデータに基づいて、建物の領域をより正確に抽出する技術を提供することを目的とする。
　［ＤＴＭ推定方法、ＤＴＭ推定プログラム、ＤＴＭ推定装置及び３次元建物モデルの作成方法について］
　（１）本発明は、航空機による地表のレーザースキャナーデータに基づいて、その所定範囲について地表面のみの標高データであるＤＴＭを、コンピュータを用いて推定するＤＴＭ推定方法であって、前記所定範囲における単位グリッド内にデータが存在しない画素を連結して河川領域を抽出し、前記河川領域を除くデータについて、第１の最大許容傾斜値を設定して暫定的なＤＴＭを推定し、推定したＤＴＭから局所的な傾斜を計算し、傾斜が所定値を超える場合は、前記第１の最大許容傾斜値より大きい第２の最大許容傾斜値を設定して再度ＤＴＭを推定することを特徴とする。
　上記のようなＤＴＭ推定方法では、ＤＴＭの推定にあたって河川領域を除くことで、河川をまたぐ地表面の探索によってＤＴＭの精度が悪くなることを防止できる。また、推定したＤＴＭにおいて局所的に傾斜が所定値を超えるときは、第１の最大許容傾斜値より大きい第２の最大許容傾斜値を設定して再度ＤＴＭを推定することにより、１回目には探索できなかった地表面データを追加し、ＤＴＭ推定の精度を向上させることができる。
　（２）また、上記（１）のＤＴＭ推定方法において、第２の最大許容傾斜値は、地形及び町並みの特徴に応じて増減されるものであってもよい。
　この場合、地形及び町並みの特徴に応じて最適なＤＴＭ推定を行うことができる。
　（３）また、本発明は、上記（１）のＤＴＭ推定方法を含む３次元建物モデルの作成方法であって、前記レーザースキャナーデータを、前記ＤＴＭ推定方法によって推定されたＤＴＭと、非地表面データとに分離し、前記非地表面データにおける各領域について、対となる向きの法線ベクトルの有無を調べ、その結果に基づいて屋根の形を推定する３次元建物モデルの作成方法である。
　この場合、精度良く推定されたＤＴＭに基づいて正確な非地表面データを取得し、さらに、対となる向きの法線ベクトルの有無を調べることにより、平坦な屋根の他、切妻、寄棟等の屋根の形を推定して３次元建物モデルを作成することができる。
　（４）また、上記（３）の３次元建物モデルの作成方法において、各領域は、所定範囲の航空写真のデータから少なくとも長方形の形状を優先的に抽出することにより予め抽出されるものであってもよい。
　この場合、屋根の形状として確率の高い長方形により、屋根の領域の識別精度を高め、より正確な３次元建物モデルを作成することができる。
　（５）一方、本発明は、航空機による地表のレーザースキャナーデータに基づいて、その所定範囲について地表面のみの標高データであるＤＴＭを推定するＤＴＭ推定プログラム（又はＤＴＭ推定プログラム記憶媒体）であって、前記所定範囲における単位グリッド内にデータが存在しない画素を連結して河川領域を抽出する機能、前記河川領域を除くデータについて、第１の最大許容傾斜値を設定して暫定的なＤＴＭを推定する機能、推定したＤＴＭから局所的な傾斜を計算する機能、及び、傾斜が所定値を超える場合は、前記第１の最大許容傾斜値より大きい第２の最大許容傾斜値を設定して再度ＤＴＭを推定する機能を、コンピュータによって実現させるためのＤＴＭ推定プログラム（又はＤＴＭ推定プログラム記憶媒体）である。
　上記のようなＤＴＭ推定プログラム（又はＤＴＭ推定プログラム記憶媒体）では、ＤＴＭの推定にあたって河川領域を除くことで、河川をまたぐ地表面の探索によってＤＴＭの精度が悪くなることを防止できる。また、推定したＤＴＭにおいて局所的に傾斜が所定値を超えるときは、第１の最大許容傾斜値より大きい第２の最大許容傾斜値を設定して再度ＤＴＭを推定することにより、１回目には探索できなかった地表面データを追加し、ＤＴＭ推定の精度を向上させることができる。
　（６）また、本発明は、航空機による地表のレーザースキャナーデータに基づいて、その所定範囲について地表面のみの標高データであるＤＴＭを推定するＤＴＭ推定装置であって、航空機による地表のレーザースキャナーデータを読み込む入力装置と、前記レーザースキャナーデータの所定範囲における単位グリッド内にデータが存在しない画素を連結して河川領域を抽出する河川抽出部と、前記河川領域を除くデータについて、第１の最大許容傾斜値を設定して暫定的なＤＴＭを推定し、推定したＤＴＭから局所的な傾斜を計算し、傾斜が所定値を超える場合は、前記第１の最大許容傾斜値より大きい第２の最大許容傾斜値を設定して再度ＤＴＭを推定するＤＴＭ推定部と、推定されたＤＴＭを表示する表示部とを備えたものである。
　上記のようなＤＴＭ推定装置では、ＤＴＭの推定にあたって河川領域を除くことで、河川をまたぐ地表面の探索によってＤＴＭの精度が悪くなることを防止できる。また、推定したＤＴＭにおいて局所的に傾斜が所定値を超えるときは、第１の最大許容傾斜値より大きい第２の最大許容傾斜値を設定して再度ＤＴＭを推定することにより、１回目には探索できなかった地表面データを追加し、ＤＴＭ推定の精度を向上させることができる。
　［領域抽出方法、領域抽出プログラム及び領域抽出装置について］
　（７）本発明は、航空機又は人工衛星から撮影した写真のデータに基づいて、コンピュータを用いて建物の領域を抽出する領域抽出方法であって、複数の異なる離散化幅を設定し、それぞれの離散化幅について、前記データの輝度値を当該離散化幅で離散設定された複数の値に離散化し、離散化して得られた離散化画像において同一値を持つ画素を連結し、建物の領域の候補として長方形に近い形状の領域を抽出し、前記複数の異なる離散化幅ごとに抽出された複数の領域群のうち、より長方形に近い形状の領域を、建物の領域として採用することを特徴とする。
　上記のような領域抽出方法では、離散化によって、輝度値の分散が大きいテクスチャを持つ領域も、１つの領域として抽出可能である。また、同一値を持つ点を連結し、長方形に近い形状の領域を抽出することで、建物以外の領域を除去し易くなる。また、複数の異なる離散化幅ごとに抽出された複数の領域群のうち、より長方形に近い形状の領域を、建物の領域として採用することにより、局所的に最適な、空間的な平滑化パラメータを適用していることに等しい機能を発揮させることができる。
　（８）また、上記（７）の領域抽出方法において、長方形に近い形状を表す指数として、（領域の面積／領域を取り囲む長方形の面積）によって定義される長方形指数を用いることができる。
　この場合、長方形に近い形状における「近さ」の度合いを、指数として簡単かつ的確に表すことができる。
　（９）また、上記（７）又は（８）の領域抽出方法において、上記抽出において、互いに隣接する領域同士が併合された場合にさらに長方形に近くなる場合は、当該併合を実行可能とするようにしてもよい。
　この場合、より正確に建物の領域を捉えることができる。
　（１０）また、上記（８）の領域抽出方法において、長方形指数が所定値より小さい場合は、建物の領域として採用しないことが好ましい。
　この場合、建物でない可能性が高い領域を抽出の対象から除外することができる。
　（１１）また、上記（７）の領域抽出方法において、ＲＧＢ輝度値に基づいて、植生と推測される領域を抽出対象から除外してもよい。
　この場合、植生の色の特徴を考慮して当該領域を抽出対象から除外することができる。
　（１２）一方、本発明は、航空機又は人工衛星から撮影した写真のデータに基づいて建物の領域を抽出する領域抽出プログラム（又は領域抽出プログラム記憶媒体）であって、複数の異なる離散化幅を設定し、それぞれの離散化幅について、前記データの輝度値を当該離散化幅で離散設定された複数の値に離散化する機能、離散化して得られた離散化画像において同一値を持つ画素を連結し、建物の領域の候補として長方形に近い形状の領域を抽出する機能、及び、前記複数の異なる離散化幅ごとに抽出された複数の領域群のうち、より長方形に近い形状の領域を、建物の領域として採用する機能を、コンピュータによって実現させるための領域抽出プログラム（又は領域抽出プログラム記憶媒体）である。
　（１３）また、本発明は、航空機又は人工衛星から撮影した写真のデータに基づいて建物の領域を抽出する領域抽出装置であって、複数の異なる離散化幅を設定し、それぞれの離散化幅について、前記データの輝度値を当該離散化幅で離散設定された複数の値に離散化する機能と、離散化して得られた離散化画像において同一値を持つ画素を連結し、建物の領域の候補として長方形に近い形状の領域を抽出する機能と、前記複数の異なる離散化幅ごとに抽出された複数の領域群のうち、より長方形に近い形状の領域を、建物の領域として採用する機能とを有するものである。
上記（１２）、（１３）の領域抽出プログラム（又は領域抽出プログラム記憶媒体）／領域抽出装置では、離散化によって、輝度値の分散が大きいテクスチャを持つ領域も、１つの領域として抽出可能である。また、同一値を持つ点を連結し、長方形に近い形状の領域を抽出することで、建物以外の領域を除去し易くなる。また、複数の異なる離散化幅ごとに抽出された複数の領域群のうち、より長方形に近い形状の領域を、建物の領域として採用することにより、局所的に最適な空間的な平滑化パラメータを適用していることに等しい機能を発揮させることができる。
　［ＤＴＭ推定方法、ＤＴＭ推定プログラム、ＤＴＭ推定装置及び３次元建物モデルの作成方法について］
　本発明によれば、ＤＴＭ推定の精度を向上させることができる。
　［領域抽出方法、領域抽出プログラム及び領域抽出装置について］
　本発明によれば、航空写真等の写真のデータに基づいて、建物の領域をより正確に抽出することができる。
（ａ）は、京都市東山区清水寺、法観寺、高台寺の周辺の航空写真であり、この地域は、建物が密集している平地と、丘陵地とが混在している。（ｂ）は、航空機によるＬｉＤＡＲデータが示す標高分布（建物等の最上部を含むＤＳＭ）を便宜的に色分けして（但し、この図ではカラーは表現されていない。）示した図である。 （ａ）は、１回目の処理で抽出された地表面データであり、（ｂ）は２回目の処理で抽出された地表面データである。 （ａ）は、１回目の処理によるＤＴＭを示し、（ｂ）は、２回目の処理を行った後の最終のＤＴＭを示す。 （ａ）は、京都市中京区五条通り周辺の航空写真であり、この地域は、建物が密集している平地の中に、河川がある。（ｂ）は、航空機によるＬｉＤＡＲデータが示す標高分布（建物等の最上部を含むＤＳＭ）を便宜的に色分けして（但し、この図ではカラーは表現されていない。）示した図である。 （ａ）は、１回目の処理で抽出された地表面データであり、（ｂ）は２回目の処理で抽出された地表面データである。 （ａ）は、１回目の処理によるＤＴＭを示し、（ｂ）は、２回目の処理を行った後の最終のＤＴＭを示す。 （ａ）は、京都市伏見区伏見桃山周辺の航空写真であり、この地域は、建物が密集している平地の中に、大きな土手を備えた河川がある。（ｂ）は、航空機によるＬｉＤＡＲデータが示す標高分布（建物等の最上部を含むＤＳＭ）を便宜的に色分けして（但し、この図ではカラーは表現されていない。）示した図である。さらに、（ｃ）は、１回目の処理によるＤＴＭを示し、（ｄ）は、２回目の処理を行った後の最終のＤＴＭを示す。 （ａ）は、図７の（ｂ）に示すＬｉＤＡＲデータの一部を拡大した図であり、（ｂ）は、（ａ）における矢印方向から見た当該地域の地上写真である。さらに、（ｃ）は、１回目の処理によるＤＴＭ（拡大図）を示し、（ｄ）は、２回目の処理を行った後の最終のＤＴＭ（拡大図）を示す。 （ａ）は、京都市東山区清水寺近くの駐車場及びその周辺の、航空機によるＬｉＤＡＲデータが示す標高分布（建物等の最上部を含むＤＳＭ）を便宜的に色分けして（但し、図ではカラーは表現されていない。）示した図であり、（ｂ）は航空写真である。 （ａ）は、１回目の処理によるＤＴＭを示し、（ｂ）は、２回目の処理を行った後の最終のＤＴＭを示す。 （ａ）は市街地の航空写真であり（７５ｍ×７５ｍ）、（ｂ）は領域分割の結果である。 （ａ）、（ｂ）は、３次元建物ワイヤーフレームモデルで、（ａ）は平面的な形状を、（ｂ）は立体的な形状を、それぞれ示している。 （ａ）は高層ビルを含む市街地の航空写真であり（７５ｍ×７５ｍ）、（ｂ）は領域分割の結果である。 （ａ）、（ｂ）は、３次元建物ワイヤーフレームモデルで、（ａ）は平面的な形状を、（ｂ）は立体的な形状を、それぞれ示している。 （ａ）は高木が建物に隣接している市街地の航空写真であり（７５ｍ×７５ｍ）、（ｂ）は領域分割の結果である。 （ａ）、（ｂ）は、３次元建物ワイヤーフレームモデルで、（ａ）は平面的な形状を、（ｂ）は立体的な形状を、それぞれ示している。 （ａ）は寄棟屋根の多い市街地の航空写真であり（７５ｍ×７５ｍ）、（ｂ）は領域分割の結果である。 （ａ）、（ｂ）は、３次元建物ワイヤーフレームモデルで、（ａ）は平面的な形状を、（ｂ）は立体的な形状を、それぞれ示している。 ＤＴＭ推定装置の構成の一例を示すブロック図である。 ＤＴＭ推定方法又は推定プログラムを示すフローチャートである。 （ａ）が切妻屋根の法線ベクトル、（ｂ）が寄棟屋根の法線ベクトルを示す図である。 家屋が密集した市街地（京都市東山区）の航空写真の一例である。 京都市東山区法観寺から東大路通り周辺を地上で撮影した写真である。 領域抽出方法を実施するためのコンピュータ装置の一例を示すブロック図である。 画素の連結の概念図である。 長方形指数算出の概念図である。 対象領域と、それに隣接する領域との統合の概念図である。 低層建物が密集する市街地（地域１）におけるエッジ抽出の結果を示す図である。 高層ビルと低層建物とが混在する市街地（地域２）におけるエッジ抽出の結果を示す図である。 高木が建物に隣接している地域（地域３）におけるエッジ抽出の結果を示す図である。 寄棟屋根作りの建物が多数存在する地域（地域４）におけるエッジ抽出の結果を示す図である。 地域１における領域抽出の結果を示す図である。 地域２における領域抽出の結果を示す図である。 地域３における領域抽出の結果を示す図である。 地域４における領域抽出の結果を示す図である。 地域１についての領域抽出の比較結果を示す図である。 地域２についての領域抽出の比較結果を示す図である。 地域３についての領域抽出の比較結果を示す図である。 地域４についての領域抽出の比較結果を示す図である。
　［ＤＴＭ推定方法、ＤＴＭ推定プログラム、ＤＴＭ推定装置及び３次元建物モデルの作成方法について］
　《ＤＴＭの推定》
　図１９は、航空機による地表のレーザースキャナーデータに基づいて、その所定範囲について地表面のみの標高データであるＤＴＭを推定する装置の構成の一例を示すブロック図である。このＤＴＭ推定装置は、ＣＰＵ１、バス２を介してＣＰＵ１と接続されたメモリ３、ハードディスク等の補助記憶装置４、インターフェース５，６、並びに、インターフェース５及び６にそれぞれ接続されたドライブ７及びディスプレイ８を含むものであり、典型的には、これは、パーソナルコンピュータである。入力装置としてのドライブ７に、ＬｉＤＡＲデータ（レーザースキャナーデータ）を収めたＣＤ，ＤＶＤ等の記憶媒体９を装着することにより、ＬｉＤＡＲデータが読込み可能となる。
　ＣＰＵ１は、ソフトウェアによって実現される内部機能として、大別すれば、河川抽出部１ａと、ＤＴＭ推定部１ｂとを備えている。これらの機能についての詳細は後述するが、河川抽出部は、ＬｉＤＡＲデータ（レーザースキャナーデータ）の所定範囲における単位グリッド内にデータが存在しない画素を連結して河川領域を抽出する機能を有する。また、ＤＴＭ推定部１ｂは、河川領域を除くデータについて、第１の最大許容傾斜値を設定して暫定的なＤＴＭを推定し、推定したＤＴＭから局所的な傾斜を計算し、傾斜が所定値を超える場合は、第１の最大許容傾斜値より大きい第２の最大許容傾斜値を設定して再度ＤＴＭを推定するという機能を有する。推定されたＤＴＭは、表示部としてのディスプレイ８に表示可能である。
　次にＣＰＵ１における上記機能すなわちＤＴＭ推定方法又は推定プログラムについて、図２０のフローチャートを参照して説明する。なお、プログラムは記憶媒体である補助記憶装置４（図１９）に記憶され、メモリ３を用いて、ＣＰＵ１により実行される。なお、ＤＴＭ推定プログラムは、各種の記憶媒体から補助記憶装置４にインストールして実行されるか、又は、ＤＴＭ推定プログラムを格納した特定のサーバからネットワーク経由でダウンロードし、補助記憶装置４にインストールして実行される。すなわち、ＤＴＭ推定プログラムは、記憶（記録）媒体としても存在・流通し得る。
　使用するＬｉＤＡＲデータは、市販品であり、例えば、平均点密度１点／１ｍ2で、３次元座標データを持つラストパルスデータ（最後に地表面で反射される受信データ）のみを含むものである。なお、使用可能なデータは、これに限られるものではない。例えば、フルウェーブタイプと呼ばれる、ファーストからラストまで連続したデータも使用可能である。
　図２０において、まず、ＣＰＵ１は、ＬｉＤＡＲデータから河川領域を抽出する（ステップＳ１）。具体的には、単位グリッドとしての１ｍグリッドに設定した画像において、データ（点）が存在しない画素を連結し、一定面積以上を占める領域を、河川領域として抽出する。
　次に、ＣＰＵ１は、ＬｉＤＡＲデータから局所最低値を抽出する（ステップＳ２）。具体的には、１回目は５０ｍ×５０ｍの窓内で局所最低値を抽出し、２回目は２５ｍ×２５ｍの窓内で局所最低値を探索する。２回目の局所最低値が１回目の局所最低値より大幅に高くならず（例えば５０ｃｍ以内）、かつ、１回目、２回目の局所最低値である３次元座標上の２点から計算される傾斜が第１の最大許容傾斜値（例えば３度）以下のときは、１回目の局所最低値を２回目の局所最低値に更新する。それ以外のときは、１回目の局所最低値を採用する。
　次に、ＣＰＵ１は、ＬｉＤＡＲデータから最小二乗法により平面を推定し（ステップＳ３）、平面に対するＲＭＳＥ（Root Mean Square Errors：平均二乗誤差）が小さい場合の、それらの点は、共通の平面上にあるものとする。具体的には、例えば、４ｍ×４ｍの領域に含まれる全ての点群を用いて平面の方程式を計算し、平面に対するＲＭＳＥが１０ｃｍ以下であれば、全ての点群に平面の方程式を付与する。この平面とは道路面の候補であり、斜面でもよいが、極端な斜面は道路の可能性が非常に低いので、例えば、平面の法線の鉛直成分が０．９以上に限定する。
　そして、ＣＰＵ１は、上記平面上の点のうち、局所最低値の点との標高差が小さく、かつ、局所最低値から見た傾斜以内である点を地表面データの初期値（初期の点群）として選定する（ステップＳ４）。
　続いて、ＣＰＵ１は、当該地表面データの近隣の点を探索し、地表面データが有する平面の方程式を使って平面までの距離を計算し、その距離が閾値（例えば１０ｃｍ）以下であれば、追加する地表面データの候補となる。そして、地表面データと成す傾斜が最大許容傾斜値（例えば３度）以下であれば、地表面データとして追加する（ステップＳ５）。
　地表面データの追加（ステップＳ５）は、追加される点が無くなるまで行われる（ステップＳ６）。追加される点が無くなると、地表面データの抽出は完了である。
　続いて、ＣＰＵ１は、地表面データが存在しない地点について近隣の地表面データを探索し、内挿してＤＴＭを推定する（ステップＳ７）。この探索時に河川領域に遭遇すれば、ＣＰＵ１は、その方向への探索を中止する。
　次に、ＣＰＵ１は、推定したＤＴＭに基づいて局所的な傾斜を計算する（ステップＳ８）。傾斜は２１ｍ×２１ｍの領域を単位として、まず地表面高さの平均値を計算する。ある領域の地表面高さに対し、周囲の８つ（上下／左右／斜め）の近傍の領域を探索し、領域の中心間の距離と地表面高さの平均値の差から傾斜を計算し、最大の傾斜値を求める。傾斜が閾値（例えば４．５度）以上の場合には、最大許容傾斜値を、より大きく設定して、第２の最大許容傾斜値に更新する（例えば４．５度に更新する。）。
　続いて、ＣＰＵ１は、ステップＳ５～Ｓ８の処理が１回目か否かを判定する（ステップＳ９）。最初は「Ｙｅｓ」であり、ＣＰＵ１は、第２の最大許容傾斜値に基づいて再びステップＳ５，Ｓ６を実行する。ここで、第１の最大許容傾斜値より大きい第２の最大許容傾斜値を設定していることによって、局所的に急な傾斜を持つ地点が、地表面データとして追加される。すなわち、最初のステップＳ５の処理では追加データから漏れていた河川の土手や丘陵地等の点群データが地表面データとして抽出され、追加される。その結果、ＤＴＭの推定が１回目よりも精度良く行われることになる。なお、２回目のステップＳ９では「Ｎｏ」となって、処理終了となる。
　上記のような処理によるＤＴＭの精度向上を、図１～１０の実例で示す。
　（実例１）
　図１の（ａ）は、京都市東山区清水寺、法観寺、高台寺の周辺の航空写真である。この地域は、建物が密集している平地と、丘陵地とが混在している。図１の（ｂ）は、航空機によるＬｉＤＡＲデータが示す標高分布（建物等の最上部を含むＤＳＭ）を便宜的に色分けして（但し、この図ではカラーは表現されていない。）示した図である。（ｂ）において、左側の黒っぽい部分が相対的に標高の低い部分であり、右側の黒っぽい部分は相対的に標高の高い部分である。中間の白っぽい部分は、それらの中間的な標高である。ＤＳＭは建物等を含んでいるため、地表面がわかりにくい。
　図２の（ａ）は、１回目の処理で抽出された地表面データであり、（ｂ）は２回目の処理で抽出された地表面データである。図中の黒い部分が地表面（地盤面）を表している。また、図３の（ａ）は、１回目の処理によるＤＴＭを示し、（ｂ）は、２回目の処理を行った後の最終のＤＴＭを示す。図２，図３のそれぞれにおいて、（ａ）、（ｂ）の比較により明らかなように、２回目の処理を行ったことによって丘陵地（向かって右側）や平地の細部のデータが格段に精度良く得られている。
　（実例２）
　図４の（ａ）は、京都市中京区五条通り周辺の航空写真である。この地域は、建物が密集している平地の中に、河川がある。図４の（ｂ）は、航空機によるＬｉＤＡＲデータが示す標高分布（建物等の最上部を含むＤＳＭ）を便宜的に色分けして（但し、この図ではカラーは表現されていない。）示した図である。
　図５の（ａ）は、１回目の処理で抽出された地表面データであり、（ｂ）は２回目の処理で抽出された地表面データである。また、図６の（ａ）は、１回目の処理によるＤＴＭを示し、（ｂ）は、２回目の処理を行った後の最終のＤＴＭを示す。この例では、１回目と、２回目とで、劇的な差は現れていない。言い換えれば、局部的に急な傾斜が無い平坦な地形では２回目の処理に１回目と同じ第１の最大許容傾斜値を用いることになるので、過剰な抽出が行われることはない。
　（実例３）
　図７の（ａ）は、京都市伏見区伏見桃山周辺の航空写真である。この地域は、建物が密集している平地の中に、大きな土手を備えた河川がある。図７の（ｂ）は、航空機によるＬｉＤＡＲデータが示す標高分布（建物等の最上部を含むＤＳＭ）を便宜的に色分けして（但し、この図ではカラーは表現されていない。）示した図である。なお、（ａ）における楕円は、高速道路であり、（ｂ）のＬｉＤＡＲデータ取得時にはまだ建設されていなかった。また、図７の（ｃ）は、１回目の処理によるＤＴＭを示し、（ｄ）は、２回目の処理を行った後の最終のＤＴＭを示す。
　図８の（ａ）は、図７の（ｂ）に示すＬｉＤＡＲデータの一部を拡大した図であり、図８の（ｂ）は、（ａ）における矢印方向から見た当該地域の地上写真である。また、図８の（ｃ）は、１回目の処理によるＤＴＭ（拡大図）を示し、（ｄ）は、２回目の処理を行った後の最終のＤＴＭ（拡大図）を示す。（ｄ）において、中央に縦に延びている白い線は河川であり、その左右の黒っぽい部分は土手である。（ｃ）、（ｄ）の比較により、１回目は地表面として扱えなかった土手が、２回目は地表面として精度良く認識されていることがわかる。
　（実例４）
　図９の（ａ）は、京都市東山区清水寺近くの駐車場及びその周辺の、航空機によるＬｉＤＡＲデータが示す標高分布（建物等の最上部を含むＤＳＭ）を便宜的に色分けして（但し、図ではカラーは表現されていない。）示した図である。（ｂ）は航空写真である。また、図１０の（ａ）は、１回目の処理によるＤＴＭを示し、（ｂ）は、２回目の処理を行った後の最終のＤＴＭを示す。（ａ）、（ｂ）の比較により、１回目は地表面としてその形状を捉えきれなかった駐車場（中央の黒い部分）が、２回目は地表面として精度良く認識されていることがわかる。
　（実例５）
　図示しないが、ビルの屋上駐車場についても同様に１回目、２回目のＤＴＭを得たが、ほとんど変化がなかった。この場合、屋上駐車場を地表面と誤って認識することはなかった。
　《まとめ》
　上記のようなＤＴＭ推定方法／プログラム（又はプログラム記憶媒体）／装置では、ＤＴＭの推定にあたって河川領域を除くことで、河川をまたぐ地表面の探索によってＤＴＭの精度が悪くなることを防止できる。また、推定したＤＴＭにおいて局所的に傾斜が所定値を超えるときは、第１の最大許容傾斜値より大きい第２の最大許容傾斜値を設定して再度ＤＴＭを推定することにより、１回目には探索できなかった地表面データを追加し、ＤＴＭ推定の精度を向上させることができる。
　なお、第２の最大許容傾斜値は、地形及び町並みの特徴に応じて増減すればよい。京都市の場合は４．５度が好適であることが確認された。京都市よりもさらに傾斜の多い土地では４．５度以上も適すると想定される。逆に、京都市よりも傾斜の少ない土地では４．５度未満（３度よりは大きい）が適すると想定される。実験的に最適な数値を選ぶことで、地形及び町並みの特徴に応じて最適なＤＴＭ推定を行うことができる。
　《３次元建物モデルの作成》
　ＬｉＤＡＲデータ及び、上述の処理により精度良く推定されたＤＴＭを用いて、（ＤＳＭ－ＤＴＭ）の演算を行うことにより、データ上で地表面を除去し、建物や木等の、地上に存在する物の高さのデータのみを取得することができる。以下、このようなデータを用いて３次元建物モデルを作成する方法について説明する。まず、航空写真のデータを用いた建物等の領域分割について説明する。
　（領域分割）
　領域分割の概要は、まず、輝度値の分散が大きいテクスチャを持つ領域も抽出するため、輝度値を少数個の値に離散化した上で領域分割を試みる。離散化する輝度値の幅は、異なる分散の程度に対応するために複数用意する。長方形に近い領域を優先的に抽出するため、長方形指数と呼ぶ指数を計算する。影で分断されている可能性を考慮し、隣接する領域と組み合わせて、長方形指数が向上すれば両者を統合する。
　以下、細部に及ぶが、領域分割の具体例について参考のため説明する。
（１）１バイト（輝度値の範囲：０～２５５）×３バンド（ＲＧＢ）画像の任意の１バンドの輝度値に対し、Ｎｄｉｓｃ種類の離散化幅を設定する。各離散化幅の下で、Ｎｏｆｆ種類の異なるオフセット値を適用し、画像輝度値を離散化する。例えば、「輝度値の離散化幅＝４０」「オフセット数＝５」の時に「オフセット値幅＝８」となり、「オフセット値＝｛０，８，１６，２４，３２｝」のＮｏｆｆ種類の離散化画像が得られる。「オフセット値＝０」の下では、原画像の輝度値が「０～３９」が同一離散化値を与えられ、同様に「４０～７９」「８０～１１９」「１２０～１５９」「１６０～１９９」「２００～２３９」「２４０～２５５」を含めた７つの区分に離散化される。本実験では次のパラメータを使用した。
　使用バンド：赤（Ｒ）バンド
　輝度値の離散化幅Δｄ＝｛４０，３０，２０｝
　オフセット数Ｎｏｆｆ＝５
　オフセット値幅Δｏｆｆ＝Δｄ／Ｎｏｆｆ＝｛８，６，４｝
（２）各離散化画像において、４方向を探索し、同じ値を持つ画素を連結し領域として抽出する。一定面積以上の大領域（本実験では６４００画素：以下、実験で採用した値）は除去する。また一定面積未満（８０画素）の小領域は周辺に一定面積以上の領域があれば統合し、なければ除去する。その後、各領域のエッジを抽出する。
（３）ある離散化幅におけるＮｏｆｆ種類の離散化画像から得られたエッジを全て１枚に重ね合わせる。
（４）一定強度以上のエッジを残し、それに連結するエッジも抽出する。この時点で本来なら閉じるはずの領域がノイズなどの影響で閉じていないことも多いため、エッジが存在しなくても周辺に直線状のエッジが存在することが確認できれば、エッジを連結する。
（５）エッジで閉じた領域にラベル番号を付与し、ＲＧＢ輝度値から判断し、植生らしい領域を除去する。（２）と同様に一定面積以上の大領域（６４００画素）を除去し、一定面積以下（本実験では３０画素）の小領域は周辺に一定面積以上の領域があれば統合し、なければ除去する。
（６）各領域の長方形指数と呼ぶ指数を次のように計算する。
（ア）領域のエッジの集合（エッジ群と呼ぶ）の２次元座標から、第１軸と第２軸を決定する。
（イ）領域の存在範囲を第１軸と第２軸の値で表現し、第１軸の（最大値－最小値＋１）、第２軸の（最大値－最小値＋１）を掛け合わせて、領域を取り囲む矩形の面積を得る。
（ウ）実際の領域の面積／領域を取り囲む矩形の面積を長方形指数と定義する。
（エ）長方形指数が一定値（０．４）を下回ると、建物でない可能性が高いと判断して除外される。
（７）ある領域Ａの近隣の領域を探索し、
（ア）併合した場合の長方形指数が各々の長方形指数よりも向上する。
（イ）指定した閾値（０．７）以上である。
（ウ）領域の平均輝度値の差が一定値（３０）以内である。
という条件を満たす近隣の領域のうち、併合時の長方形指数が最大となる領域を併合の候補とし、仮に領域Ｂとする。領域Ｂでも同様に探索し、領域Ａが（ア）～（ウ）の全ての条件を満たし、かつ併合時の長方形指数が最大となる領域である場合に、領域ＡとＢを併合する。
（８）Ｎｄｉｓｃ種類の離散化幅の下で得られた一定面積以上の領域を、長方形指数が高い順に選んでいく。ただし、当該領域の一部でも既に選ばれた領域に重なっていれば、その領域は選ばない。
（９）選ばれなかった領域に対し、再度選定していく。今度は既に選ばれた領域との重なりが一定値未満（本実験では、３０％）で、かつ一定面積未満（３０画素未満）であれば、重なっていない部分だけを新たな領域として追加選出する。ただし重なっていない部分に対しても長方形指数を算出し、閾値（０．４５）以上である場合に限る。
（１０）領域内の穴を埋める。
　（３次元建物モデル）
　領域分割で得られた領域は屋根単位ではなく、同一建物で対となる傾斜のある屋根、あるいは別の建物の屋根が分離されていないこともある。また領域分割では屋根・建物の抽出漏れ、逆に樹木や道路も抽出していることも起こりうる。しかしながら、下記の対策を取ることで、モデリングの漏れを少なくしている。
　推定したＤＴＭは、ＬｉＤＡＲデータから除去する。屋根の平面、法線ベクトルを計算する際にＲＭＳＥの許容値を設けることで樹木の大半は削除できる。切妻屋根の建物において向かい合う屋根が分離されていない場合、法線ベクトルの分布状況を調べて混合している可能性が高い場合には分離して、切妻屋根を生成する。隣接する切妻屋根の建物で片側の屋根が分離されていない状況において（屋根の数が棟を挟んで１対２の関係）、分離されていない領域を２つの屋根に分離することで２組の１対１の屋根のペアを生成し、２棟の切妻屋根の建物を生成する。
　密集市街地における建物の配置に着目し、屋根面がもつ傾斜や法線ベクトルといった情報を利用することで、切妻屋根や寄棟屋根、平屋根の判定を行い、より本来の形に近い建物３次元モデルを生成する。具体的には、近隣の点群間の水平・鉛直距離が一定値以内の点群から主たる方位角を決めて、それに基づいて建物の棟や輪郭を決定することで、現実に近い建物の向きを決定できる。
　具体的には、３次元座標の点群データであるＬｉＤＡＲデータに対し、ＤＴＭを用いて地表面と非地表面を分離するフィルタリング処理を行う。非地表面データから屋根や建物の領域における法線ベクトルを計算する。９ｍ×９ｍの領域内で、６点以上点が存在する場合に平面の方程式を計算し、ＲＭＳＥが１０ｃｍ以下であれば、計算に使用した全ての点群に法線ベクトルを付与する。領域の法線が傾いていれば対となる領域を探索する。
　非地表面データに対し、最近隣の点データからの水平距離、鉛直距離が一定値以内である点群ごとに分類する。このうち屋根や建物の領域に該当する点群だけを用いて、点群データの主たる方位角を決定する。寄棟屋根、切妻屋根、それ以外の屋根の順にワイヤーフレームモデルを作成する。各モデルを作成する際には棟や輪郭が、点群データの主たる方位角と同一あるいは直交するように配慮する。
　寄棟屋根については、対となる屋根の組が２組存在する時に生成される。切妻屋根については、対となる屋根が存在する領域があり、両者の領域において切妻屋根を生成する。法線ベクトル分布が複数の法線の混合状態を示す領域については、領域内において２つに分割し、法線ベクトルが明瞭に分離される場合には、切妻屋根を生成する。そうでない場合には平屋根を生成する。なお、図２１は、（ａ）が切妻屋根の法線ベクトル、（ｂ）が寄棟屋根の法線ベクトルを示す。
　それ以外の屋根については、平屋根を生成する。建物として抽出されなかったものの、地表面から一定以上の高さを持ち、特定の平面上にある点を選定し、最終結果に含める。
　上記のような処理による領域分割、３次元建物モデルの作成を、図１１～１８の実例で示す。
　（実例１）
　図１１の（ａ）は市街地の航空写真である（７５ｍ×７５ｍ）。（ｂ）は領域分割の結果である。便宜上色分けしている（本図では色は表示されていない。）が、色自体に意味はない。（ｂ）に示すように屋根の輪郭が概ね捉えられている。また、切妻、寄棟の特徴が捉えられている個所もある。図１２の（ａ）、（ｂ）は、３次元建物ワイヤーフレームモデルで、（ａ）は平面的な形状を、（ｂ）は立体的な形状を、それぞれ示している。平坦、切妻、寄棟等の屋根の形状が明瞭に現れていることが分かる。
　（実例２）
　図１３の（ａ）は高層ビルを含む市街地の航空写真である（７５ｍ×７５ｍ）。（ｂ）は領域分割の結果である。便宜上色分けしている（本図では色は表示されていない。）が、色自体に意味はない。（ｂ）に示すように屋根の輪郭が概ね捉えられている。また、切妻、寄棟の特徴が捉えられている個所もある。図１４の（ａ）、（ｂ）は、３次元建物ワイヤーフレームモデルで、（ａ）は平面的な形状を、（ｂ）は立体的な形状を、それぞれ示している。平坦、切妻、寄棟等の屋根の形状が明瞭に現れていることが分かる。
　（実例３）
　図１５の（ａ）は高木が建物に隣接している市街地の航空写真である（７５ｍ×７５ｍ）。（ｂ）は領域分割の結果である。便宜上色分けしている（本図では色は表示されていない。）が、色自体に意味はない。（ｂ）に示すように屋根の輪郭が概ね捉えられている。また、切妻、寄棟の特徴が捉えられている個所もある。図１６の（ａ）、（ｂ）は、３次元建物ワイヤーフレームモデルで、（ａ）は平面的な形状を、（ｂ）は立体的な形状を、それぞれ示している。平坦、切妻、寄棟等の屋根の形状が明瞭に現れていることが分かる。
　（実例４）
　図１７の（ａ）は寄棟屋根の多い市街地の航空写真である（７５ｍ×７５ｍ）。（ｂ）は領域分割の結果である。便宜上色分けしている（本図では色は表示されていない。）が、色自体に意味はない。（ｂ）に示すように屋根の輪郭が概ね捉えられている。また、切妻、寄棟の特徴が捉えられている個所もある。図１８の（ａ）、（ｂ）は、３次元建物ワイヤーフレームモデルで、（ａ）は平面的な形状を、（ｂ）は立体的な形状を、それぞれ示している。平坦、切妻、寄棟等の屋根の形状が明瞭に現れていることが分かる。
　《まとめ》
　以上のように、精度良く推定されたＤＴＭに基づいて正確な非地表面データを取得し、さらに、対となる向きの法線ベクトルの有無を調べることにより、平坦な屋根の他、切妻、寄棟等の屋根の形を推定して３次元建物モデルを作成することができる。
　また、３次元建物モデルを作成するための領域を予め抽出することにより、屋根の領域の識別精度を高め、より正確な３次元建物モデルを作成することができる。
　なお、今回開示された実施の形態はすべての点で例示であって制限的なものではないと考えられるべきである。本発明の範囲は特許請求の範囲によって示され、特許請求の範囲と均等の意味及び範囲内での全ての変更が含まれることが意図される。
　［領域抽出方法、領域抽出プログラム及び領域抽出装置について］
　以下、本発明の一実施形態に係る領域抽出方法・領域抽出プログラムについて、説明する。領域抽出のモデル地域としては、京都市東山区法観寺から東大路通り周辺を対象とした。この地域では木造建築が敷地間口いっぱいに立ち並んでいる。図２３は、当該地域を地上で撮影した写真である。また、この地域の航空写真データとして、２５ｃｍ解像度、オルソ（正射投影）処理済みのデータ（市販品）を使用した。このデータは、平面直角座標系に合わせてある。なお、オルソ処理されていても、若干は建物の側面が見えるが、実用上は特に問題ない。
　次に、建物の領域（屋根の輪郭形状）を抽出する領域抽出方法・領域抽出プログラムについて説明する。図２４は、領域抽出方法を実施するためのコンピュータ装置１０の一例を示すブロック図であり、また、領域抽出プログラムは、このコンピュータ装置１０に機能を実現させるものである。また、領域抽出プログラムの機能が搭載されたコンピュータ装置１０は、本発明の一実施形態に係る領域抽出装置である。
　このコンピュータ装置１０は、ＣＰＵ１と、バス２を介してＣＰＵ１と接続されたメモリ３と、ハードディスク等の補助記憶装置４と、インターフェース５，６と、インターフェース５及び６にそれぞれ接続されたドライブ７及びディスプレイ８とを含むものであり、典型的には、これは、パーソナルコンピュータである。入力装置としてのドライブ７に、航空写真を収めたＣＤ，ＤＶＤ等の記憶媒体９を装着することにより、航空写真のデータが読込み可能となる。なお、領域抽出プログラムは、航空写真とは別に、各種の記憶媒体から補助記憶装置４にインストールして実行されるか、又は、領域抽出プログラムを格納した特定のサーバからネットワーク経由でダウンロードし、補助記憶装置４にインストールして実行される。すなわち、領域抽出プログラムは、記憶（記録）媒体としても存在・流通し得る。
　《領域抽出の手順》
　以下、領域抽出（方法・プログラムの主要ステップ）について詳細に説明する。
　領域抽出の概要としては、輝度値の分散が大きいテクスチャを持つ屋根も効果的に抽出すべく、輝度値を所定数の値に離散化し、同一の離散化値を持つ領域をラベリングする。そして、平面視した建物に見られる長方形に近い形状の領域を優先的に抽出する。また、複数の異なる離散化幅を適用して得られた領域群に対し、より長方形に近い領域から順次採用していくことで、局所的に最適な、空間的な平滑化パラメータを適用したのと同様な処理を実現する。以下、ステップごとに詳細に説明する。
　（第１ステップ）
　まず、１バイト（輝度値の範囲:０～２５５）×３バンド（ＲＧＢ）画像の任意の１バンドの輝度値に対し、Ｎdisc種類の離散化幅を設定する。各離散化幅の下で、Ｎoff種類の異なるオフセット値を適用し、画像輝度値を離散化する。例えば、「輝度値の離散化幅＝４０」「オフセット数＝５」のときに「オフセット値幅＝８」となり、「オフセット値＝｛０，８，１６，２４，３２｝」のＮoff種類の離散化画像が得られる。「オフセット値＝０」の下では、原画像の輝度値「０～３９」が同一離散化値を与えられ、同様に「４０～７９」「８０～１１９」「１２０～１５９」「１６０～１９９」「２００～２３９」「２４０～２５５」を含めた７つの区分に離散化される。本実験では一例として次のパラメータを使用した。
　使用バンド：赤（Ｒ）バンド
　輝度値の離散化幅Δｄ＝｛４０，３０，２０｝
　オフセット数Ｎoff＝５
　オフセット値幅Δoff＝Δｄ／Ｎoff＝｛８，６，４｝
　（第２ステップ）
　次に、各離散化画像において、４方向（上下左右）を探索し、同じ値を持つ画素を連結し領域として抽出する。一定面積以上の大領域（例えば６４００画素以上）は除去する。また一定面積未満（例えば８０画素未満）の小領域は周辺に一定面積以上の領域があれば統合し、なければ除去する。その後、各領域のエッジ（輪郭）を抽出する。
　（第３ステップ）
　次に、ある離散化幅におけるＮoff種類の離散化画像から得られたエッジを全て１枚に重ね合わせる。
　（第４ステップ）
　次に、一定強度以上のエッジを残し、それに連結するエッジも抽出する。この時点で隣接する屋根や建物がほぼ同様の輝度値を持つ場合にそれらのエッジが抽出されていないことも多いため、エッジが存在しなくても周辺に直線状のエッジが存在することが確認できれば、エッジを連結する。
　このエッジ連結について、補足説明する。まず、エッジではない画素の周辺８方向（上下左右斜め）を探索し、一定画素数ｄの中に存在するエッジの数を調べる。
　図２５は、左上の方向にあるエッジ群グループ１ａと、右下の方向にあるエッジ群グループ１ｂとを探索する例を示している。例えば、（ａ）グループ１ａにｄ１以上、かつ、グループ１ｂにもｄ１以上にエッジ群が存在し、（ｂ）グループ２ａにｄ２以下、かつ、グループ２ｂにもｄ２以下のエッジ群しか存在しない、という条件を用意し、条件（ａ）、（ｂ）を共に満足する場合に、非エッジである中心画素をエッジとして補完する。（ｂ）の条件があることで、矩形領域の隅に近い領域内部の画素を誤ってエッジと抽出することを防いでいる。探索は、上下方向、左右方向、左上から右下、右上から左下の４回にわたって調べる。今回は、ｄ＝７、ｄ１＝５、ｄ２＝１と設定した。
　（第５ステップ）
　次に、エッジで閉じた領域にラベル番号を付与（ラベリング）し、ＲＧＢ輝度値から判断して、植生らしい領域を除去する。ステップ２と同様に、一定面積以上の大領域（例えば６４００画素以上）を除去し、一定面積以下（例えば３０画素以下）の小領域は周辺に一定面積以上の領域があれば統合し、なければ除去する。
　植生の除去について補足説明する。対象地域では、植生には緑系統と赤系統の植生が確認できた。そのため、
（１）Bgrn ≧ Tveg, DN, Bblue ≧ Tveg, DN, and Bgrn / Bblue ≧ Tveg, g2b_ratio
（２）Bred ≧ Tveg, DN, Bblue ≧ Tveg, DN, and Bred / Bblue ≧ Tveg, r2b_ratio
のいずれかを満足する画素が占有する割合Rgrn_veg, Rred_vegを計算する。ここでBred、 Bgrn 、Bblueは赤、緑、青バンドの輝度値を、Tveg, DNは輝度値の閾値を、Tveg, g2b_ratio、Tveg, r2b_ratioは指標の閾値を意味する。Rgrn_veg, Rred_vegのいずれかが閾値Tveg, ratio以上であれば、植生領域として除去する。ここでは、Tveg, DN ＝２０、Tveg, g2b_ratio＝１．２５、Tveg, r2b_ratio＝２．０、Tveg, ratio＝０．３の値を採用した。
　（第６ステップ）
　次に、各領域の長方形指数と呼ぶ指数を、次のように計算する。
　(ア)領域のエッジの集合（エッジ群と呼ぶ。）の２次元座標から、第１軸及び第２軸を決定する。
　(イ)領域の存在範囲を第１軸と第２軸の値で表現し、第１軸の(最大値－最小値＋１)、第２軸の(最大値－最小値＋１)を掛け合わせて、領域を取り囲む矩形の面積を得る。
　(ウ)（実際の領域の面積／領域を取り囲む矩形の面積）を長方形指数と定義する。
　ここで、長方形指数が一定値（例えば０．４）を下回ると、建物でない可能性が高いと判断して、抽出の対象から除外する。
　長方形指数について補足説明する。図２６に長方形指数算出の概念図を示す。ある領域のエッジ群から第１，第２軸を計算し、各々の辺が第１，第２軸に平行な長方形のうち、対象領域を取り囲む、図示のような最小の長方形を求める。この第１，第２軸の決定には、２点間の距離が一定値（５画素以上２０画素以下）にあるエッジ群を用いて直線の傾きを投票し、最大の頻度が発生する傾きを第１軸の向きとし、第２軸は第１軸に直交するように定めた。
　ｉdx　＝　Ｓactual／Ｓrect　　　・・・（１）
と定義する。ここでｉdxは長方形指数を、Ｓactualはある領域の面積を、Ｓrectはその領域を取り囲む長方形の面積を表す。長方形指数は０から１の範囲の値を取り、１に近いほど長方形に近い形状を有する。このような指数により、長方形に近い形状における「近さ」の度合いを、指数として簡単かつ的確に表すことができる。
　（第７ステップ）
　式（１）の長方形指数を用いて、ステップ７では、影で分断されている屋根・建物を統合し、本来の領域を推定する。
　例えば、ある領域Ａの近隣の領域を探索し、
　(ア)併合した場合の長方形指数が各々の長方形指数よりも向上する。
　(イ)指定した閾値（０．７）以上である。
　(ウ)領域の平均輝度値の差が一定値（３０）以内である。
という条件を満たすか否かを判定する。満たす場合には、近隣の領域のうち、併合時の長方形指数が最大となる領域を併合の候補とし、これを、仮に領域Ｂとする。領域Ｂでも同様に探索し、領域Ａが（ア）～（ウ）の全ての条件を満たし、かつ併合時の長方形指数が最大となる領域である場合に、領域ＡとＢとを互いに併合する。
　例えば図２７に示すように、対象領域αに隣接する領域β、γ、δの全てを対象に、仮に統合した場合の長方形指数を計算する。この際の第１，第２軸は、2つの領域のエッジ群のうち、互いの領域に近いエッジ群は除外したエッジ群を用いて計算される。軸の決定には、２点間の距離が一定値（５画素以上２０画素以下）にあるエッジ群を用いて直線の傾きを投票し、最大の頻度が発生する傾きを第１軸の向きとし、第２軸は第１軸に直交するように定めた。
　（第８ステップ）
　次に、Ｎdisc種類の離散化幅の下で得られた一定面積以上の領域を、長方形指数が高い順に選んでいく。ただし、当該領域の一部でも既に選ばれた領域に重なっていれば、その領域は選ばない。
　（第９ステップ）
　選ばれなかった領域に対し、再度選定していく。今度は既に選ばれた領域との重なりが一定値未満（３０％）で、かつ一定面積（８０画素）未満であれば、重なっていない部分だけを新たな領域として追加選出する。ただし重なっていない部分に対しても長方形指数を算出し、閾値（０．４５）以上である場合に限る。
　（第１０ステップ）
　最後に、領域内の穴を埋める。
　《実例》
　以下、実例を示す。図２８は、低層建物が密集する市街地（以下、地域１という。）におけるエッジ抽出の結果を示す図である。（ａ）は撮影の解像度３００×３００画素で、実寸法約７５ｍ×７５ｍのエリアを撮影した航空写真である。（ｂ）は、（ａ）のデータに対して既知のキャニー（Canny）フィルタを用いた結果である。（ｃ）は、上述の離散化において離散化幅４０を用いた際のエッジ抽出結果（第１段階）である。さらに、（ｄ）は、離散化幅３０を用いた際のエッジ抽出結果（第２段階）、（ｅ）は、離散化幅２０を用いた際のエッジ抽出結果（第３段階）である。
　図２９は、高層ビルと低層建物とが混在する市街地（以下、地域２という。）におけるエッジ抽出の結果を示す図である。（ａ）～（ｅ）は、図２８と同様であり、（ａ）は航空写真、（ｂ）はキャニーフィルタを用いた結果、（ｃ）は離散化幅４０を用いた際のエッジ抽出結果（第１段階）、（ｄ）は離散化幅３０を用いた際のエッジ抽出結果（第２段階）、（ｅ）は離散化幅２０を用いた際のエッジ抽出結果（第３段階）である。
　図３０は、高木が建物に隣接している地域（以下、地域３という。）におけるエッジ抽出の結果を示す図である。（ａ）～（ｅ）は、図２８と同様であり、（ａ）は航空写真、（ｂ）はキャニーフィルタを用いた結果、（ｃ）は離散化幅４０を用いた際のエッジ抽出結果（第１段階）、（ｄ）は離散化幅３０を用いた際のエッジ抽出結果（第２段階）、（ｅ）は離散化幅２０を用いた際のエッジ抽出結果（第３段階）である。
　図３１は、寄棟屋根作りの建物が多数存在する地域（以下、地域４という。）におけるエッジ抽出の結果を示す図である。（ａ）～（ｅ）は、図２８と同様であり、（ａ）は航空写真、（ｂ）はキャニーフィルタを用いた結果、（ｃ）は離散化幅４０を用いた際のエッジ抽出結果（第１段階）、（ｄ）は離散化幅３０を用いた際のエッジ抽出結果（第２段階）、（ｅ）は離散化幅２０を用いた際のエッジ抽出結果（第３段階）である。
　図３２，図３３，図３４，図３５はそれぞれ、上述の地域１（図７），地域２（図２９），地域３（図３０），地域４（図３１）における領域抽出の結果を示す図である。図３２～３５の各図において、（ａ）は各地域の航空写真、（ｂ）は離散化幅４０を用いた際の領域抽出結果（第１段階）、（ｃ）は離散化幅３０を用いた際の領域抽出結果（第２段階）、（ｄ）は離散化幅２０を用いた際のエッジ抽出結果（第３段階）、（ｅ）は、長方形指数を用いて３つの結果（（ｂ）～（ｄ））を統合した最終結果の図である。図３２～３５のいずれにおいても、（ｅ）は、（ｂ）～（ｄ）の各結果と比較すると明らかなように、局所的に３つの結果の最も良いところを取った状態となり、建物の領域（エッジ）がより正確に抽出されている。
　図３６は、地域１についての図であり、（ａ）～（ｆ）はそれぞれ、（ａ）航空写真、（ｂ）領域抽出結果を比較するための参照点、（ｃ）領域抽出結果（主たる方向を計算する際に主成分分析を利用）、（ｄ）領域抽出結果（主たる方向を計算する際に２点間を通る直線の傾きを投票して決定）、（ｅ）ITT Visual Information Solutions の画像処理ソフトENVI EXでScale=50, Merge=50に設定した時の領域抽出結果、（ｆ）ENVI EXでScale=50, Merge=80に設定した時の領域抽出結果、である。なお、（ｃ）～（ｆ）はカラー画像では建物のエッジが明瞭に表示されるが、この図では建物の縁取りのように見えているのがエッジである。
　なお、ENVI EXの「Feature Extraction」という機能で設定を要求されるパラメータは、「Scale Level」と「Merge Level」である。本例の対象地域では分散値の大きいテクスチャを持つ屋根が多い。「Merge Level」を変更することで最終的に残される領域の大きさや数が変わるため、「Scale Level = 50.0」は共通とし、「Merge Level = 50.0」「Merge Level = 80.0」の2種類の値で実行した。
　図３７は、地域２についての図であり、（ａ）～（ｆ）はそれぞれ、（ａ）航空写真、（ｂ）領域抽出結果を比較するための参照点、（ｃ）領域抽出結果（主たる方向を計算する際に主成分分析を利用）、（ｄ）領域抽出結果（主たる方向を計算する際に２点間を通る直線の傾きを投票して決定）、（ｅ）ENVI EXでScale=50, Merge=50に設定した時の領域抽出結果、（ｆ）ENVI EXでScale=50, Merge=80に設定した時の領域抽出結果、である。
　図３８は、地域３についての図であり、（ａ）～（ｆ）はそれぞれ、（ａ）航空写真、（ｂ）領域抽出結果を比較するための参照点、（ｃ）領域抽出結果（主たる方向を計算する際に主成分分析を利用）、（ｄ）領域抽出結果（主たる方向を計算する際に２点間を通る直線の傾きを投票して決定）、（ｅ）ENVI EXでScale=50, Merge=50に設定した時の領域抽出結果、（ｆ）ENVI EXでScale=50, Merge=80に設定した時の領域抽出結果、である。
　図３９は、地域４についての図であり、（ａ）～（ｆ）はそれぞれ、（ａ）航空写真、（ｂ）領域抽出結果を比較するための参照点、（ｃ）領域抽出結果（主たる方向を計算する際に主成分分析を利用）、（ｄ）領域抽出結果（主たる方向を計算する際に２点間を通る直線の傾きを投票して決定）、（ｅ）ENVI EXでScale=50, Merge=50に設定した時の領域抽出結果、（ｆ）ENVI EXでScale=50, Merge=80に設定した時の領域抽出結果、である。
　図３６～３９のそれぞれにおいて、中段の（ｃ），（ｄ）と下段の（ｅ），（ｆ）との比較により明らかなように、本実施形態の領域抽出による（ｃ），（ｄ）では、建物の長方形の領域が、より正確に抽出されていることがわかる。ENVI EXによる（ｅ），（ｆ）では不要なエッジや、細かすぎるエッジが非常に多く、（ｃ），（ｄ）に比べると領域抽出の出来映えが良くない。
　《考察》
　上記実施形態に係る領域抽出の手法では異なる離散化幅を適用して得られたエッジを統合しているが、このことは異なる空間解像度へ変換し処理することと本質的には同じである。ただし、オフセット値を変えて適用し統合することで、単なる平滑化フィルタと異なりエッジは保存されている。何より重要なのは、複数の異なる離散化幅を適用し得られた領域群に対し、長方形指数の大きい領域から順次採用していくことは局所的に最適な、空間的な平滑化パラメータを適用していることに等しい、ということである。
　図３２～３５からは、統合する処理を通じて、屋根や建物の大きさに応じて局所的に最適な空間スケールパラメータを選択できていることが示されている。ただし、離散化後の画像を用いたラベリングには時間を要し、画像サイズが大きくなるほど処理時間が増大する。例えば、クロック３．２ＧＨｚで動作するＣＰＵに、６ＧＢのメモリを使用するコンピュータの場合、各対象地域での平均処理時間は約９０秒であった。
　なお、領域のエッジから算出される形状に関する指標を活用して、屋根や建物に見られる長方形の形状の領域を優先的に抽出することが好ましい。その際にエッジを抽出しても完全に閉じない屋根や建物が多く、このことが領域抽出の精度を下げる要因であるが、完全に閉じていないエッジ群に対しては、長方形や三角形の形状の可能性が高い場合にエッジを閉じる処理を加えることが好ましい。
　また上述の、例えば３つの異なる離散化幅を用いることで、特にΔｄ＝２０の結果（図２８～３１における（ｅ））から図３６～３９の参照点１－ａ，１－ｂ，３－ａ，４－ａのような影がかかっている屋根のエッジを明瞭に抽出できていることがわかる。この場合でもエッジを閉じる処理がなければ分離されないことも確認できており、エッジ補完（連結）処理の有効性が確認できた。すなわち、エッジ補完処理と、異なる離散化幅の結果の統合という２つの要因が効果を発揮して、領域抽出精度が向上したと考えられる。ただし、ステップ７に示した「領域の平均輝度値の差が一定値以内」という条件がないと、明瞭に分かれている屋根も過剰に統合する結果が確認された。
　寄棟屋根に見られる三角形の領域も、図３６～３９の結果からは良好に抽出されていることが分かる。理想的な三角形の長方形指数は０．５に留まり、優先的に採用される訳ではない。今回良好に抽出できた要因として、三角形領域周辺の長方形領域が良好に抽出された点が挙げられる。離散化幅やオフセット値によっては、三角形領域と長方形領域が統合されることもある。しかしながら、３種類の離散化幅で得られた領域抽出結果を長方形指数という尺度で評価すると、三角形領域と長方形領域が統合される領域の長方形指数は、本来の長方形単独の領域の長方形指数より低く採用されにくい。このような選考過程に基づいて、長方形領域が先に選抜され、続いて三角形領域が抽出された。
　既存の領域抽出手法の結果に比べて、上述の手法の結果では領域が存在しない範囲が広い。これは領域面積に上限を設けたことも一因であるが、領域の長方形指数が一定値（０．４）未満であると屋根・建物の可能性が低いとして除去していることに大きく起因している。全てではないが大半の植生や、図３６の墓地を含む領域は形状が長方形には程遠く、結果的に除去できている。すなわち、建物領域抽出という観点では十分に機能していると言える。
　長方形指数を領域面積の観点から述べる。本手法では長方形指数が高い領域から順次抽出していくため、比較的小さな領域が優先的に抽出されることが多くなり、多少長方形指数が小さくても大まかな外形を捉えている領域が選ばれないこともありえる。そのような大まかな外形を捉えたい応用に対応するため、その方法を検討した。例えば、下記の式（２）を用いた補正を実行することで、一部の領域では図３６～３９に示す結果よりも大きく、複数の屋根・建物に対応した領域が得られた。
　Δｉdx＝Ｃ1・ｌｎ（Ｓactual・Ｃ2＋Ｃ3）　　　・・・（２）
　ここでΔｉdxは長方形指数の補正値、Ｃ1～Ｃ3は経験的に定める係数を表す。前述の対象地域に適用して経験的に確認した結果、Ｃ1＝０．０５，Ｃ2＝１００．０，Ｃ3＝１．０とした。式（２）の係数の中でもＣ1が特に最終抽出結果に大きな影響を与えることが判明した。Ｃ1＝１．０と設定すると、補正値が大きすぎて道路やそれに繋がる大きな面積を持つ領域の長方形指数が高くなり、それに応じて優先的に選ばれるようになり、本来目的としている屋根・建物の抽出数が減少する結果となった。今回は伝統的な日本建物が建ち並ぶ地域を対象にしているが、この長方形指数補正における関数や、パラメータ値は、抽出したい建物を念頭に調整するべきである。
　一方、長方形指数を算出する場合、図３６～３９の（ｃ）には主成分分析で軸を決めた場合の領域抽出の結果が示されている。参照点２－ａ、３－ｂ、３－ｃ、４－ｂにおいては、主成分分析を適用するとスレート屋根の領域が分割されて、一部が欠けて抽出された。完全な長方形でない限り、エッジ群を用いて主成分分析を適用すると、特に今回対象にしているような影で分断されている屋根では、第１主成分は屋根の辺に平行になるとは限られない結果が得られた。その反面、１－ｃでは主成分分析を用いた領域抽出結果では、４つの屋根が一つの大きな領域として抽出されている。
　このように主成分分析を用いることで領域抽出結果の安定性を欠くことが確認できたため、主成分分析を適用せず、２点間の距離が一定値にあるエッジ群を用いて直線の傾きを投票し、最大の頻度が発生する傾きを第１軸の向きとし、第２軸は第１軸に直交するように定めた。この２点間の距離の上限、下限を変更すると一部領域抽出の結果にも変化が生じたが、主成分分析を用いた結果ほど大きな変化は見られなかった。この閾値は対象地域の特性に応じて経験的に決めることが求められる。
　最後に、植生領域の除去について述べる。本実施形態では市街地の建物抽出が目的であり、植生は除去すべき対象であった。領域抽出の前処理として、画素単位の輝度値を参照して植生らしさが高い画素を除去する方法も考えられる。しかしながら、屋根や道路に部分的にかかる植生を除去することで本来の屋根や建物の形状が得られず、分断されたり、あるいはそもそも領域が小さくなりすぎて抽出されなくなったりすることがある。
　そこで、検討の結果、植生領域も領域抽出の対象に含めながら処理を進め、最終段階で適用することにした。この方針は現実的に高精度の抽出に貢献しているといえる。例えば赤色の植生を除去しようとすると、赤色の屋根が誤って除去されることがあった。そのため、赤色植生を慎重に除去すべく、閾値をTveg, r2b_ratio＝２．０と高めに設定した。同様の発想で影も一つの領域として抽出し、長方形指数が向上する場合には近隣の領域と統合するように対処した。影らしい領域には道路が多数含まれており、輝度値を用いて除去することも検討したが、誤って屋根や建物も除去される事例が確認できたため、最終的に除去しなかった。
　以上のように、本実施形態の領域抽出によれば、密集市街地で効果的に建物あるいは屋根を領域として抽出できる。当該手法では、輝度値の分散が大きいテクスチャを持つ屋根も効果的に抽出するため、輝度値を少数個の値に離散化し、同一の離散化値を持つ領域をラベリングする。そして、領域のエッジから算出される長方形指数を活用して、屋根や建物に見られる長方形の形状の領域を優先的に抽出するようにした。また完全に閉じていないエッジ群に対し、長方形や三角形の形状の可能性が高い場合にエッジを閉じる処理を加えた。
　京都市の伝統的な建物が密集する地域が撮影された２５ｃｍ解像度の航空写真を用いた実験では、３つの異なる離散化幅を適用したが、エッジ補完処理と異なる離散化幅の結果の統合という２つの要因が効果を発揮して、既存手法では十分に抽出できにくい影領域も明瞭に抽出されるようになった。当該手法の最も重要な特徴としては、複数の異なる離散化幅を適用し得られた領域群に対し、長方形指数の大きい領域から順次採用していくことで、局所的に最適な、空間的な平滑化パラメータを適用した処理を実現している点である。
　また寄棟屋根に見られる三角形の領域は、長方形指数は低く優先的に採用される訳ではないが、結果的には良好に抽出できた。三角形領域と長方形領域が統合される領域の長方形指数は、本来の長方形単独の領域の長方形指数より低く採用されにくく、長方形領域が先に選抜されることで、続いて三角形領域が抽出された。最終的に、低層建物が密集する市街地、高層ビルと低層建物が混在する市街地、高木が建物に隣接している地域、寄棟屋根作りの建物が多数存在する地域の全ての対象地域において、従来の領域抽出手法より良好な結果を得ることができた。よって、輝度値の離散化と長方形指数の活用に特徴を有する手法は有用であることが確認できた。
　なお、上記実施形態では元データとして航空写真を利用したが、航空写真に代えて、高精度に撮影された人工衛星から撮影した写真のデータを用いることも可能である。
　また、今回開示された実施の形態はすべての点で例示であって制限的なものではないと考えられるべきである。本発明の範囲は特許請求の範囲によって示され、特許請求の範囲と均等の意味及び範囲内での全ての変更が含まれることが意図される。
　［ＤＴＭ推定方法、ＤＴＭ推定プログラム、ＤＴＭ推定装置及び３次元建物モデルの作成方法について］
１ａ：河川抽出部
１ｂ：ＤＴＭ推定部
７：ドライブ（入力装置）
８：ディスプレイ（表示部）
　［領域抽出方法、領域抽出プログラム及び領域抽出装置について］
１０：コンピュータ装置（領域抽出装置）

Claims (13)

  1. 　航空機による地表のレーザースキャナーデータに基づいて、その所定範囲について地表面のみの標高データであるＤＴＭを、コンピュータを用いて推定するＤＴＭ推定方法であって、
    　前記所定範囲における単位グリッド内にデータが存在しない画素を連結して河川領域を抽出し、
    　前記河川領域を除くデータについて、第１の最大許容傾斜値を設定して暫定的なＤＴＭを推定し、
    　推定したＤＴＭから局所的な傾斜を計算し、
    　傾斜が所定値を超える場合は、前記第１の最大許容傾斜値より大きい第２の最大許容傾斜値を設定して再度ＤＴＭを推定する
    　ことを特徴とするＤＴＭ推定方法。
  2. 　前記第２の最大許容傾斜値は、地形及び町並みの特徴に応じて増減される請求項１記載のＤＴＭ推定方法。
  3. 　請求項１のＤＴＭ推定方法を含む３次元建物モデルの作成方法であって、
    　前記レーザースキャナーデータを、前記ＤＴＭ推定方法によって推定されたＤＴＭと、非地表面データとに分離し、
    　前記非地表面データにおける各領域について、対となる向きの法線ベクトルの有無を調べ、その結果に基づいて屋根の形を推定する３次元建物モデルの作成方法。
  4. 　前記各領域は、前記所定範囲の航空写真のデータから少なくとも長方形の形状を優先的に抽出することにより予め抽出されるものである請求項３記載の３次元建物モデルの作成方法。
  5. 　航空機による地表のレーザースキャナーデータに基づいて、その所定範囲について地表面のみの標高データであるＤＴＭを推定するＤＴＭ推定プログラムであって、
    　前記所定範囲における単位グリッド内にデータが存在しない画素を連結して河川領域を抽出する機能、
    　前記河川領域を除くデータについて、第１の最大許容傾斜値を設定して暫定的なＤＴＭを推定する機能、
    　推定したＤＴＭから局所的な傾斜を計算する機能、及び、
    　傾斜が所定値を超える場合は、前記第１の最大許容傾斜値より大きい第２の最大許容傾斜値を設定して再度ＤＴＭを推定する機能を、コンピュータによって実現させるためのＤＴＭ推定プログラム。
  6. 　航空機による地表のレーザースキャナーデータに基づいて、その所定範囲について地表面のみの標高データであるＤＴＭを推定するＤＴＭ推定装置であって、
    　航空機による地表のレーザースキャナーデータを読み込む入力装置と、
    　前記レーザースキャナーデータの所定範囲における単位グリッド内にデータが存在しない画素を連結して河川領域を抽出する河川抽出部と、
    　前記河川領域を除くデータについて、第１の最大許容傾斜値を設定して暫定的なＤＴＭを推定し、推定したＤＴＭから局所的な傾斜を計算し、傾斜が所定値を超える場合は、前記第１の最大許容傾斜値より大きい第２の最大許容傾斜値を設定して再度ＤＴＭを推定するＤＴＭ推定部と、
    　推定されたＤＴＭを表示する表示部と
    　を備えていることを特徴とするＤＴＭ推定装置。
  7. 　航空機又は人工衛星から撮影した写真のデータに基づいて、コンピュータを用いて建物の領域を抽出する領域抽出方法であって、
    　複数の異なる離散化幅を設定し、それぞれの離散化幅について、前記データの輝度値を当該離散化幅で離散設定された複数の値に離散化し、
    　離散化して得られた離散化画像において同一値を持つ画素を連結し、建物の領域の候補として長方形に近い形状の領域を抽出し、
    　前記複数の異なる離散化幅ごとに抽出された複数の領域群のうち、より長方形に近い形状の領域を、建物の領域として採用する
    　ことを特徴とする領域抽出方法。
  8. 　長方形に近い形状を表す指数として、（領域の面積／領域を取り囲む長方形の面積）によって定義される長方形指数を用いる請求項７記載の領域抽出方法。
  9. 　前記抽出において、互いに隣接する領域同士が併合された場合にさらに長方形に近くなる場合は、当該併合を実行可能とする請求項７又は８に記載の領域抽出方法。
  10. 　長方形指数が所定値より小さい場合は、建物の領域として採用しない請求項８記載の領域抽出方法。
  11. 　ＲＧＢ輝度値に基づいて、植生と推測される領域を抽出対象から除外する請求項７記載の領域抽出方法。
  12. 　航空機又は人工衛星から撮影した写真のデータに基づいて建物の領域を抽出する領域抽出プログラムであって、
    　複数の異なる離散化幅を設定し、それぞれの離散化幅について、前記データの輝度値を当該離散化幅で離散設定された複数の値に離散化する機能、
    　離散化して得られた離散化画像において同一値を持つ画素を連結し、建物の領域の候補として長方形に近い形状の領域を抽出する機能、及び、
    　前記複数の異なる離散化幅ごとに抽出された複数の領域群のうち、より長方形に近い形状の領域を、建物の領域として採用する機能を、コンピュータによって実現させるための領域抽出プログラム。
  13. 　航空機又は人工衛星から撮影した写真のデータに基づいて建物の領域を抽出する領域抽出装置であって、
    　複数の異なる離散化幅を設定し、それぞれの離散化幅について、前記データの輝度値を当該離散化幅で離散設定された複数の値に離散化する機能と、
    　離散化して得られた離散化画像において同一値を持つ画素を連結し、建物の領域の候補として長方形に近い形状の領域を抽出する機能と、
    　前記複数の異なる離散化幅ごとに抽出された複数の領域群のうち、より長方形に近い形状の領域を、建物の領域として採用する機能と
    　を有することを特徴とする領域抽出装置。
PCT/JP2012/061204 2011-06-09 2012-04-26 Ｄｔｍ推定方法、ｄｔｍ推定プログラム、ｄｔｍ推定装置及び３次元建物モデルの作成方法、並びに、領域抽出方法、領域抽出プログラム及び領域抽出装置 WO2012169294A1 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
US14/122,082 US20140081605A1 (en) 2011-06-09 2012-04-26 Dtm estimation method, dtm estimation program, dtm estimation device, and method for creating 3-dimensional building model, and region extraction method, region extraction program, and region extraction device

Applications Claiming Priority (4)

Application Number Priority Date Filing Date Title
JP2011-129316 2011-06-09
JP2011129316A JP2012256230A (ja) 2011-06-09 2011-06-09 Ｄｔｍ推定方法、ｄｔｍ推定プログラム及びｄｔｍ推定装置、並びに、３次元建物モデルの作成方法
JP2011144266A JP2013012034A (ja) 2011-06-29 2011-06-29 領域抽出方法、領域抽出プログラム及び領域抽出装置
JP2011-144266 2011-06-29

Publications (1)

Publication Number Publication Date
WO2012169294A1 true WO2012169294A1 (ja) 2012-12-13

Family

ID=47295865

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2012/061204 WO2012169294A1 (ja) 2011-06-09 2012-04-26 Ｄｔｍ推定方法、ｄｔｍ推定プログラム、ｄｔｍ推定装置及び３次元建物モデルの作成方法、並びに、領域抽出方法、領域抽出プログラム及び領域抽出装置

Country Status (2)

Country Link
US (1) US20140081605A1 (ja)
WO (1) WO2012169294A1 (ja)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2019111956A1 (ja) * 2017-12-08 2019-06-13 アジア航測株式会社 地物地盤高別色付画像生成装置及び地物高別色付画像生成プログラム
CN109887085A (zh) * 2019-02-22 2019-06-14 中国科学院地理科学与资源研究所 一种河流主支流分级方法及河流主支流分级装置
CN111788602A (zh) * 2017-12-29 2020-10-16 菲力尔系统公司 点云去噪系统和方法

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104751479A (zh) * 2015-04-20 2015-07-01 中测新图(北京)遥感技术有限责任公司 基于tin数据的建筑物提取方法和装置
EP3403050A4 (en) 2016-02-15 2019-08-21 Pictometry International Corp. AUTOMATED SYSTEM AND METHOD OF CHARACTERIZATION EXTRACTION
CN106204611B (zh) * 2016-07-19 2018-12-28 中国科学院地理科学与资源研究所 一种基于HASM模型的LiDAR点云数据处理方法及装置
CN109085560B (zh) * 2018-08-13 2022-07-15 陕西科技大学 一种提取激光雷达特征点的方法
US11144758B2 (en) 2018-11-15 2021-10-12 Geox Gis Innovations Ltd. System and method for object detection and classification in aerial imagery
CN111582194B (zh) * 2020-05-12 2022-03-29 吉林大学 基于多特征lstm网络的多时相高分辨率遥感影像的建筑物提取方法
CN111814715B (zh) * 2020-07-16 2023-07-21 武汉大势智慧科技有限公司 一种地物分类方法及装置

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH1083452A (ja) * 1996-09-09 1998-03-31 Kokusai Gijutsu Kaihatsu Kk パターン欠陥検出装置
JP2002157576A (ja) * 2000-11-22 2002-05-31 Nec Corp ステレオ画像処理装置及びステレオ画像処理方法並びにステレオ画像処理用プログラムを記録した記録媒体
JP2004341422A (ja) * 2003-05-19 2004-12-02 Hitachi Ltd 地図作成装置、地図配信方法及び地図作成プログラム
JP2007188117A (ja) * 2006-01-11 2007-07-26 Asahi Koyo Kk 崩壊地抽出方法、装置及びプログラム
JP2008146314A (ja) * 2006-12-08 2008-06-26 Hitachi Software Eng Co Ltd 濃淡画像のラベリング方法、及びそれを実行するためのプログラム
JP2009301510A (ja) * 2008-06-17 2009-12-24 Asahi Koyo Kk 傾斜分析装置、方法及びプログラム

Family Cites Families (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4906940A (en) * 1987-08-24 1990-03-06 Science Applications International Corporation Process and apparatus for the automatic detection and extraction of features in images and displays
US6018587A (en) * 1991-02-21 2000-01-25 Applied Spectral Imaging Ltd. Method for remote sensing analysis be decorrelation statistical analysis and hardware therefor
US5265173A (en) * 1991-03-20 1993-11-23 Hughes Aircraft Company Rectilinear object image matcher
US5640468A (en) * 1994-04-28 1997-06-17 Hsu; Shin-Yi Method for identifying objects and features in an image
US6404920B1 (en) * 1996-09-09 2002-06-11 Hsu Shin-Yi System for generalizing objects and features in an image
JP2000293696A (ja) * 1999-04-07 2000-10-20 Matsushita Electric Ind Co Ltd 画像認識装置
US8224078B2 (en) * 2000-11-06 2012-07-17 Nant Holdings Ip, Llc Image capture and identification system and process
JP4172941B2 (ja) * 2002-03-01 2008-10-29 日立ソフトウエアエンジニアリング株式会社 土地区画データ作成方法および装置
US7400770B2 (en) * 2002-11-06 2008-07-15 Hrl Laboratories Method and apparatus for automatically extracting geospatial features from multispectral imagery suitable for fast and robust extraction of landmarks
US8249346B2 (en) * 2008-01-28 2012-08-21 The United States Of America As Represented By The Secretary Of The Army Three dimensional imaging method and apparatus
US8224097B2 (en) * 2008-06-12 2012-07-17 Sri International Building segmentation for densely built urban regions using aerial LIDAR data
US8073279B2 (en) * 2008-07-08 2011-12-06 Harris Corporation Automated atmospheric characterization of remotely sensed multi-spectral imagery
US20130300740A1 (en) * 2010-09-13 2013-11-14 Alt Software (Us) Llc System and Method for Displaying Data Having Spatial Coordinates

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH1083452A (ja) * 1996-09-09 1998-03-31 Kokusai Gijutsu Kaihatsu Kk パターン欠陥検出装置
JP2002157576A (ja) * 2000-11-22 2002-05-31 Nec Corp ステレオ画像処理装置及びステレオ画像処理方法並びにステレオ画像処理用プログラムを記録した記録媒体
JP2004341422A (ja) * 2003-05-19 2004-12-02 Hitachi Ltd 地図作成装置、地図配信方法及び地図作成プログラム
JP2007188117A (ja) * 2006-01-11 2007-07-26 Asahi Koyo Kk 崩壊地抽出方法、装置及びプログラム
JP2008146314A (ja) * 2006-12-08 2008-06-26 Hitachi Software Eng Co Ltd 濃淡画像のラベリング方法、及びそれを実行するためのプログラム
JP2009301510A (ja) * 2008-06-17 2009-12-24 Asahi Koyo Kk 傾斜分析装置、方法及びプログラム

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
KATSUHIRO KITAGAWA: "Fundamental studies on the terrain analysis in mountainous forest areas", BULLETIN OF THE NAGOYA UNIVERSITY FORESTS, vol. 11, December 1991 (1991-12-01), pages 39 - 192 *

Cited By (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2019111956A1 (ja) * 2017-12-08 2019-06-13 アジア航測株式会社 地物地盤高別色付画像生成装置及び地物高別色付画像生成プログラム
CN111465971A (zh) * 2017-12-08 2020-07-28 亚洲航测株式会社 按地物地基高度着色图像生成装置以及按地物高度着色图像生成程序
KR20200096951A (ko) * 2017-12-08 2020-08-14 아시아코소쿠 가부시키가이샤 지물 지반 높이별 컬러 화상 생성 장치 및 지물 높이별 컬러 화상 생성 프로그램
AU2018381378B2 (en) * 2017-12-08 2021-10-07 Asia Air Survey Co., Ltd. Feature/Ground Height-Based Colored Image Generating Apparatus and Feature Height-Based Colored Image Generating Program
US11181367B2 (en) 2017-12-08 2021-11-23 Asia Air Survey Co., Ltd. Feature/ground height-based colored image generating apparatus and feature height-based colored image generating program
CN111465971B (zh) * 2017-12-08 2022-03-22 亚洲航测株式会社 按地物地基高度着色图像生成装置以及按地物高度着色图像生成方法
KR102415768B1 (ko) 2017-12-08 2022-06-30 아시아코소쿠 가부시키가이샤 지물 지반 높이별 컬러 화상 생성 장치 및 지물 높이별 컬러 화상 생성 프로그램
TWI787410B (zh) * 2017-12-08 2022-12-21 日商亞洲航測股份有限公司 依地物地盤高著色之畫像的產生裝置及依地物高著色之畫像的產生程式
CN111788602A (zh) * 2017-12-29 2020-10-16 菲力尔系统公司 点云去噪系统和方法
CN111788602B (zh) * 2017-12-29 2024-05-28 泰立戴恩菲力尔有限责任公司 点云去噪系统和方法
CN109887085A (zh) * 2019-02-22 2019-06-14 中国科学院地理科学与资源研究所 一种河流主支流分级方法及河流主支流分级装置
CN109887085B (zh) * 2019-02-22 2019-11-19 中国科学院地理科学与资源研究所 一种河流主支流分级方法及河流主支流分级装置

Also Published As

Publication number Publication date
US20140081605A1 (en) 2014-03-20

Similar Documents

Publication Publication Date Title
WO2012169294A1 (ja) Ｄｔｍ推定方法、ｄｔｍ推定プログラム、ｄｔｍ推定装置及び３次元建物モデルの作成方法、並びに、領域抽出方法、領域抽出プログラム及び領域抽出装置
CN107735794B (zh) 使用图像处理的状况检测
Vosselman Point cloud segmentation for urban scene classification
CN106780524B (zh) 一种三维点云道路边界自动提取方法
Zhou et al. Object-based land cover classification of shaded areas in high spatial resolution imagery of urban areas: A comparison study
CN107835997B (zh) 使用计算机视觉的用于电力线走廊监测的植被管理
Awrangjeb et al. Automatic extraction of building roofs using LIDAR data and multispectral imagery
Vosselman et al. Recognising structure in laser scanner point clouds
Xiao et al. Image-based street-side city modeling
Cheng et al. Integration of LiDAR data and optical multi-view images for 3D reconstruction of building roofs
Kedzierski et al. Methods of laser scanning point clouds integration in precise 3D building modelling
US20210201570A1 (en) Method and apparatus for generating digital surface model using satellite imagery
Tarsha Kurdi et al. Automatic filtering and 2D modeling of airborne laser scanning building point cloud
Wendel et al. Unsupervised facade segmentation using repetitive patterns
CN111652241A (zh) 融合影像特征与密集匹配点云特征的建筑物轮廓提取方法
JP7418281B2 (ja) 地物の分類システム、分類方法及びそのプログラム
CN111476723B (zh) 一种Landsat-7扫描线纠正器失效的遥感图像丢失像素恢复方法
Zakharov et al. Automatic building detection from satellite images using spectral graph theory
JP2013012034A (ja) 領域抽出方法、領域抽出プログラム及び領域抽出装置
CN111582156A (zh) 一种基于倾斜摄影城市三维模型的高大植被提取方法
JP2012256230A (ja) Ｄｔｍ推定方法、ｄｔｍ推定プログラム及びｄｔｍ推定装置、並びに、３次元建物モデルの作成方法
Novacheva Building roof reconstruction from LiDAR data and aerial images through plane extraction and colour edge detection
Bulatov et al. On Applications of Sequential Multi-view Dense Reconstruction from Aerial Images.
RU2488881C2 (ru) Способ определения линии поверхности земли
Cord et al. Automatic extraction and modelling of urban buildings from high resolution aerial images

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 12796085

Country of ref document: EP

Kind code of ref document: A1

WWE Wipo information: entry into national phase

Ref document number: 14122082

Country of ref document: US

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 12796085

Country of ref document: EP

Kind code of ref document: A1