Goals
The primary objectives for this tutorial:
- How to use OpenCV imread to load satellite imagery.
- 使用imread函数读取卫星图像
- How to use OpenCV imread to load SRTM Digital Elevation Models
- 使用imread函数加载SRTM数字海拔模型
- Given the corner coordinates of both the image and DEM, correllate the elevation data to the image to find elevations for each pixel.
- 找到每一个坐标的海拔高度。
- Show a basic, easy-to-implement example of a terrain heat map.
- 显示一个地热图？
- Show a basic use of DEM data coupled with ortho-rectified imagery.
由于原图像无法找到，我的图像是8UC3，并不是16SC1，于是对代码进行了修改，首先进行了类型转换，然后使用split函数分离出一个通道。此外，这次并没有使用KDevelop开发环境，而是使用的qt creator进行的开发。其中重要的是，需要在pro后缀的文件中添加OpenCV的一些库。如下：
库的目录一般是：/usr/local/lib
opencv的路径一般是：/usr/local/include
这些应该是当初CMAKE确定的。
添加的内容如下：
INCLUDEPATH += /usr/local/include \
/usr/local/include/opencv \
/usr/local/include/opencv2
LIBS += /usr/local/lib/libopencv_calib3d.so \
/usr/local/lib/libopencv_core.so \
/usr/local/lib/libopencv_features2d.so \
/usr/local/lib/libopencv_flann.so \
/usr/local/lib/libopencv_highgui.so \
/usr/local/lib/libopencv_imgcodecs.so \
/usr/local/lib/libopencv_imgproc.so \
/usr/local/lib/libopencv_ml.so \
/usr/local/lib/libopencv_objdetect.so \
/usr/local/lib/libopencv_photo.so \
/usr/local/lib/libopencv_shape.so \
/usr/local/lib/libopencv_stitching.so \
/usr/local/lib/libopencv_superres.so \
/usr/local/lib/libopencv_videoio.so \
/usr/local/lib/libopencv_video.so \
/usr/local/lib/libopencv_videostab.so
Code
代码有点繁琐，其实本节内容就只需要掌握imread和imwrite即可，而这些之前都已经学过了。
/*
* gdal_image.cpp -- Load GIS data into OpenCV Containers using the Geospatial Data Abstraction Library
*/
// OpenCV Headers
#include "opencv2/core.hpp"
#include "opencv2/imgproc.hpp"
#include "opencv2/highgui.hpp"
// C++ Standard Libraries
#include <cmath>
#include <iostream>
#include <stdexcept>
#include <vector>
using namespace std;
// define the corner points
// Note that GDAL library can natively determine this
cv::Point2d tl( -122.441017, 37.815664 );
cv::Point2d tr( -122.370919, 37.815311 );
cv::Point2d bl( -122.441533, 37.747167 );
cv::Point2d br( -122.3715, 37.746814 );
// determine dem corners
cv::Point2d dem_bl( -122.0, 38);
cv::Point2d dem_tr( -123.0, 37);
// range of the heat map colors
std::vector<std::pair<cv::Vec3b,double> > color_range;
//关于pair的用法可以参考《C++ Primer》的379页
// List of all function prototypes
cv::Point2d lerp( const cv::Point2d&, const cv::Point2d&, const double& );
cv::Vec3b get_dem_color( const double& );
cv::Point2d world2dem( const cv::Point2d&, const cv::Size&);
cv::Point2d pixel2world( const int&, const int&, const cv::Size& );
void add_color( cv::Vec3b& pix, const uchar& b, const uchar& g, const uchar& r );
/*
* Linear Interpolation
* p1 - Point 1
* p2 - Point 2
* t - Ratio from Point 1 to Point 2
*/
cv::Point2d lerp( cv::Point2d const& p1, cv::Point2d const& p2, const double& t ){
return cv::Point2d( ((1-t)*p1.x) + (t*p2.x),
((1-t)*p1.y) + (t*p2.y));
}
/*
* Interpolate Colors
*/
//模板函数，可以参考《C++ primer》578页
template <typename DATATYPE, int N>
cv::Vec<DATATYPE,N> lerp( cv::Vec<DATATYPE,N> const& minColor,
cv::Vec<DATATYPE,N> const& maxColor,
double const& t ){
cv::Vec<DATATYPE,N> output;
for( int i=0; i<N; i++ ){
output[i] = (uchar)(((1-t)*minColor[i]) + (t * maxColor[i]));
}
return output;
}
/*
* Compute the dem color
*/
cv::Vec3b get_dem_color( const double& elevation ){
// if the elevation is below the minimum, return the minimum
if( elevation < color_range[0].second ){
return color_range[0].first;
}
// if the elevation is above the maximum, return the maximum
if( elevation > color_range.back().second ){
return color_range.back().first;
}
// otherwise, find the proper starting index
int idx=0;
double t = 0;
for( int x=0; x<(int)(color_range.size()-1); x++ ){
// if the current elevation is below the next item, then use the current
// two colors as our range
if( elevation < color_range[x+1].second ){
idx=x;
t = (color_range[x+1].second - elevation)/
(color_range[x+1].second - color_range[x].second);
break;
}
}
// interpolate the color
return lerp( color_range[idx].first, color_range[idx+1].first, t);
}
/*
* Given a pixel coordinate and the size of the input image, compute the pixel location
* on the DEM image.
*/
cv::Point2d world2dem( cv::Point2d const& coordinate, const cv::Size& dem_size ){
// relate this to the dem points
// ASSUMING THAT DEM DATA IS ORTHORECTIFIED
double demRatioX = ((dem_tr.x - coordinate.x)/(dem_tr.x - dem_bl.x));
double demRatioY = 1-((dem_tr.y - coordinate.y)/(dem_tr.y - dem_bl.y));
cv::Point2d output;
output.x = demRatioX * dem_size.width;
output.y = demRatioY * dem_size.height;
return output;
}
/*
* Convert a pixel coordinate to world coordinates
*/
cv::Point2d pixel2world( const int& x, const int& y, const cv::Size& size ){
// compute the ratio of the pixel location to its dimension
double rx = (double)x / size.width;
double ry = (double)y / size.height;
// compute LERP of each coordinate
cv::Point2d rightSide = lerp(tr, br, ry);
cv::Point2d leftSide = lerp(tl, bl, ry);
// compute the actual Lat/Lon coordinate of the interpolated coordinate
return lerp( leftSide, rightSide, rx );
}
/*
* Add color to a specific pixel color value
*/
void add_color( cv::Vec3b& pix, const uchar& b, const uchar& g, const uchar& r ){
if( pix[0] + b < 255 && pix[0] + b >= 0 ){ pix[0] += b; }
if( pix[1] + g < 255 && pix[1] + g >= 0 ){ pix[1] += g; }
if( pix[2] + r < 255 && pix[2] + r >= 0 ){ pix[2] += r; }
}
/*
* Main Function
*/
int main(){
/*
* Check input arguments
*/
// if( argc < 3 ){
// cout << "usage: " << argv[0] << " <image_name> <dem_model_name>" << endl;
// return -1;
// }
// load the image (note that we don't have the projection information. You will
// need to load that yourself or use the full GDAL driver. The values are pre-defined
// at the top of this file
cv::Mat image = cv::imread("lena.jpg", cv::IMREAD_LOAD_GDAL | cv::IMREAD_COLOR );
// load the dem model
cv::Mat dem = cv::imread("test1.jpg", cv::IMREAD_LOAD_GDAL | cv::IMREAD_ANYDEPTH );
dem.convertTo(dem,CV_16S);
cv::Mat temp[3];
cv::split(dem,temp);
// create our output products
cv::Mat output_dem(image.size(), CV_8UC3 );
cv::Mat output_dem_flood(image.size(), CV_8UC3 );
//cout << temp[0].type() << endl;
//cout << temp[1].channels() << endl;
//cout << temp[2].depth() << endl;
//for sanity sake, make sure GDAL Loads it as a signed short
if( temp[0].type() != CV_16SC1 ){ throw std::runtime_error("DEM image type must be CV_16SC1"); }
// define the color range to create our output DEM heat map
// Pair format ( Color, elevation ); Push from low to high
// Note: This would be perfect for a configuration file, but is here for a working demo.
color_range.push_back( std::pair<cv::Vec3b,double>(cv::Vec3b( 188, 154, 46), -1));
color_range.push_back( std::pair<cv::Vec3b,double>(cv::Vec3b( 110, 220, 110), 0.25));
color_range.push_back( std::pair<cv::Vec3b,double>(cv::Vec3b( 150, 250, 230), 20));
color_range.push_back( std::pair<cv::Vec3b,double>(cv::Vec3b( 160, 220, 200), 75));
color_range.push_back( std::pair<cv::Vec3b,double>(cv::Vec3b( 220, 190, 170), 100));
color_range.push_back( std::pair<cv::Vec3b,double>(cv::Vec3b( 250, 180, 140), 200));
// define a minimum elevation
double minElevation = -10;
// iterate over each pixel in the image, computing the dem point
for( int y=0; y<image.rows; y++ ){
for( int x=0; x<image.cols; x++ ){
// convert the pixel coordinate to lat/lon coordinates
cv::Point2d coordinate = pixel2world( x, y, image.size() );
// compute the dem image pixel coordinate from lat/lon
cv::Point2d dem_coordinate = world2dem( coordinate, temp[0].size() );
// extract the elevation
double dz;
if( dem_coordinate.x >= 0 && dem_coordinate.y >= 0 &&
dem_coordinate.x < temp[0].cols && dem_coordinate.y < temp[0].rows ){
dz = temp[0].at<short>(dem_coordinate);
}else{
dz = minElevation;
}
// write the pixel value to the file
output_dem_flood.at<cv::Vec3b>(y,x) = image.at<cv::Vec3b>(y,x);
// compute the color for the heat map output
cv::Vec3b actualColor = get_dem_color(dz);
output_dem.at<cv::Vec3b>(y,x) = actualColor;
// show effect of a 10 meter increase in ocean levels
if( dz < 10 ){
add_color( output_dem_flood.at<cv::Vec3b>(y,x), 90, 0, 0 );
}
// show effect of a 50 meter increase in ocean levels
else if( dz < 50 ){
add_color( output_dem_flood.at<cv::Vec3b>(y,x), 0, 90, 0 );
}
// show effect of a 100 meter increase in ocean levels
else if( dz < 100 ){
add_color( output_dem_flood.at<cv::Vec3b>(y,x), 0, 0, 90 );
}
}}
// print our heat map
cv::imwrite( "heat-map.jpg" , output_dem );
// print the flooding effect image
cv::imwrite( "flooded.jpg", output_dem_flood);
imshow("out", output_dem);
imshow("put", output_dem_flood);
cv::waitKey(0);
return 0;
}