On Saturday, 15. November 2008 22:13:35 [EMAIL PROTECTED] wrote:
Oliver Eichler wrote:
Just a small add on. Fabrice tried your patch and now he can grow rice in
the hills. Thus RasterIO will make staircases for interpolation.
so I don't commit it to svn;).
yes :) You have to do the interpolation with my little 2x2 filter.
I added this filter, but don't see other comments yet.
Now it work faster and look the same as before changes.
I don't commit this patch, because want think more about the first comment.
Oliver
On Saturday, 15. November 2008 20:53:42 [EMAIL PROTECTED] wrote:
Hi, Oliver
Another topic is to optimize the copy action of the elevation data.
It's less a question of how to do it, more one of who is doing it.
Either way, me or you, is ok for me. It's just important to assign
it to one, to avoid double work.
I want investigate this task. I have some ideas and need make some
experiments.
Ok, it's yours. I will ammuse myself with stacking raster, vector and
DEM layers this weekend. Need to find a good concept for that.
I'm create first version of the optimization patch. It's draft only,
but reflect the main idea: "get elevation of all points and make one
request of RasterIO".
Hi Andrew,
just my 2 cents from looking over the patch:
I would have done the API a bit different by passing a buffer to
IMap::getElevation(). Of course you need some additional data like the
buffer dimensions, topLeft, bottomRight and step size in meter for x
and y. By that you can get the elevation matrix with one single call.
And you can allocate the buffer on the stack instead on the heap with
new/delete. I don't like new and delete on simple objects as this is
prone to memory leaks.
To let RasterIO do the interpolation might be a good idea. It will do it
fast. But I don't think it does a good one. If you do an overzoom into a
raster map you will notice that the pixels will get squares and these
squares just get larger. That won't fit for elevation data. A 90m SRTM
data will become like stairs on a 1m resolution map. The interpolation
you disabled (float ele = w.c1 * e[0] + w.c2 * e[1] + w.c3 * e[2] + w.c4
* e[3];) will do better.
But maybe RasterIO recognize the special character of the DEM data and
applies a better interpolation. It's worth a try.
Just for the records: Contour shading does quite the same. It just does
not look as staircasey because the array is smoothed by a lowpass
filter. But it would have been too bad to derive the elevation value
from it. That is why I made this 4 point interpolation. Once we can
interpolate a whole region fast, I will use it for shading. That should
make much better results.
I can't apply the patch locally as my source tree is on hiatus for
Garmin typ file implementation. But I will try after that.
Oliver
------------------------------------------------------------------------
- This SF.Net email is sponsored by the Moblin Your Move Developer's
challenge Build the coolest Linux based applications with Moblin SDK &
win great prizes Grand prize is a trip for two to an Open Source event
anywhere in the world
http://moblin-contest.org/redirect.php?banner_id=100&url=/
_______________________________________________
QLandkarte-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/qlandkarte-users
-------------------------------------------------------------------------
This SF.Net email is sponsored by the Moblin Your Move Developer's
challenge Build the coolest Linux based applications with Moblin SDK &
win great prizes Grand prize is a trip for two to an Open Source event
anywhere in the world
http://moblin-contest.org/redirect.php?banner_id=100&url=/
_______________________________________________
QLandkarte-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/qlandkarte-users
>From 4d971841c823407223edd59e1bc288f3528c94ab Mon Sep 17 00:00:00 2001
From: Adnrew <[EMAIL PROTECTED]>
Date: Sun, 16 Nov 2008 01:35:27 +0300
Subject: [PATCH] add determination "Region".
it's needed for getting values for matrix of points.
this patch add funtionality for get elevations and
draw3DMap will work quicker.
It's only draft. We should add smooth out for elevations.
Update: add the interpolation with Oliver's little 2x2 filter
if DEM data have resolution less than required.
---
src/CMap3DWidget.cpp | 56 +++++++++++++++++-----------
src/CMap3DWidget.h | 3 +
src/CMapDEM.cpp | 100 ++++++++++++++++++++++++++++++++++++++++++++++++++
src/CMapDEM.h | 7 +++
src/IMap.h | 4 ++
5 files changed, 148 insertions(+), 22 deletions(-)
diff --git a/src/CMap3DWidget.cpp b/src/CMap3DWidget.cpp
index 305f2b3..ea90021 100644
--- a/src/CMap3DWidget.cpp
+++ b/src/CMap3DWidget.cpp
@@ -42,6 +42,7 @@ CMap3DWidget::CMap3DWidget(QWidget * parent)
zRot = 0;
xRotSens = 0.3;
zRotSens = 0.3;
+ step = 5;
xShift = 0;
yShift = 0;
@@ -215,6 +216,28 @@ void CMap3DWidget::drawFlatMap()
glDisable(GL_TEXTURE_2D);
}
+void CMap3DWidget::setEleRegion()
+{
+ QSize s = map->getSize();
+ double w = s.width();
+ double h = s.height();
+
+ IMap& dem = CMapDB::self().getDEM();
+ XY p1, p2;
+ p1.u = 0;
+ p1.v = 0;
+ p2.u = w;
+ p2.v = h;
+ map->convertPt2Rad(p1.u, p1.v);
+ map->convertPt2Rad(p2.u, p2.v);
+ dem.setRegion(p1, p2, w/step + 1, h/step + 1);
+}
+
+void CMap3DWidget::deleteEleRegion()
+{
+ IMap& dem = CMapDB::self().getDEM();
+ dem.deleteRegion();
+}
void CMap3DWidget::draw3DMap()
{
@@ -222,13 +245,12 @@ void CMap3DWidget::draw3DMap()
double w = s.width();
double h = s.height();
- int i, iv, it, j, k, cur, end;
- double step = 5;
+ int iv, it, j, k, end;
double x, y, u, v;
GLdouble *vertices;
GLdouble *texCoords;
GLuint idx[4];
- int num = (w / step + 2);
+ int num = (w / step + 1);
glEnableClientState(GL_VERTEX_ARRAY);
glEnableClientState(GL_TEXTURE_COORD_ARRAY);
glEnable(GL_TEXTURE_2D);
@@ -250,6 +272,7 @@ void CMap3DWidget::draw3DMap()
glColor3f(1.0, 0.0, 0.0);
glPointSize(2.0);
+ setEleRegion();
for (y = 0; y < h - step; y += step) {
it = it % (num * 4);
end = it + num * 2;
@@ -260,12 +283,7 @@ void CMap3DWidget::draw3DMap()
v = y;
texCoords[it + 0] = u / w;
texCoords[it + 1] = 1 - v / h;
- map->convertPt2Rad(u, v);
- // next step of optimization will be get elevation for the set of points
- // read line from DEM file is more effectively
- vertices[iv + 2] = dem.getElevation(u, v);
- if (vertices[iv + 2] > 5000 || vertices[iv + 2] < 0)
- qDebug() << vertices[iv + 0] << " " << vertices[iv + 1] << " " << vertices[iv + 2] << end;
+ vertices[iv + 2] = dem.getRegionValue(x/step, y/step);
convertPt23D(vertices[iv + 0], vertices[iv + 1], vertices[iv + 2]);
}
@@ -283,6 +301,7 @@ void CMap3DWidget::draw3DMap()
for (j = 0; j < 4; j++)
idx[j]++;
}
+ deleteEleRegion();
delete [] vertices;
delete [] texCoords;
glDisable(GL_TEXTURE_2D);
@@ -350,32 +369,25 @@ void CMap3DWidget::drawTrack()
void CMap3DWidget::updateElevationLimits()
{
- double x, y, u, v, ele;
+ double x, y, ele;
QSize s = map->getSize();
double w = s.width();
double h = s.height();
- double step = 5;
IMap& dem = CMapDB::self().getDEM();
- u = 0;
- v = 0;
- map->convertPt2Rad(u, v);
- minElevation = maxElevation = dem.getElevation(u,v);
+ setEleRegion();
+ minElevation = maxElevation = dem.getRegionValue(0, 0);
for (y = 0; y < h - step; y += step)
for (x = 0; x < w; x += step) {
- u = x;
- v = y;
- map->convertPt2Rad(u, v);
- // FIXME can't use map instead of dem. need investigation.
- ele = dem.getElevation(u,v);
+ ele = dem.getRegionValue(x / step, y /step);
if (ele > maxElevation)
maxElevation = ele;
if (ele < minElevation)
minElevation = ele;
}
-
+ deleteEleRegion();
if (! track.isNull() && (maxElevation - minElevation < 1)) {
/*selected track exist and dem isn't present for this map*/
QList<CTrack::pt_t>& trkpts = track->getTrackPoints();
@@ -556,7 +568,7 @@ void CMap3DWidget::convertDsp2Z0(QPoint &a)
{
GLdouble projection[16];
GLdouble modelview[16];
- GLdouble k1, z1, x1, y1, x0, xk, y0, yk, z0, zk;
+ GLdouble k1, z1, x0, xk, y0, yk, z0, zk;
GLint viewport[4];
GLsizei vx, vy;
diff --git a/src/CMap3DWidget.h b/src/CMap3DWidget.h
index b52d020..3a87f0d 100644
--- a/src/CMap3DWidget.h
+++ b/src/CMap3DWidget.h
@@ -54,6 +54,8 @@ class CMap3DWidget: public QGLWidget
void focusOutEvent ( QFocusEvent * event );
void createActions();
void updateElevationLimits();
+ void setEleRegion();
+ void deleteEleRegion();
QAction *eleZoomInAct;
QAction *eleZoomOutAct;
@@ -64,6 +66,7 @@ class CMap3DWidget: public QGLWidget
QSet<int> pressedKeys;
private:
+ double step;
QPointer<IMap> map;
void loadMap();
/// expand map relative to the center
diff --git a/src/CMapDEM.cpp b/src/CMapDEM.cpp
index 6d7c1e3..7cbd724 100644
--- a/src/CMapDEM.cpp
+++ b/src/CMapDEM.cpp
@@ -16,6 +16,7 @@
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111 USA
**********************************************************************************************/
+#include <assert.h>
#include "CMapDEM.h"
#include "CMapDB.h"
@@ -188,6 +189,105 @@ void CMapDEM::dimensions(double& lon1, double& lat1, double& lon2, double& lat2)
{
}
+void CMapDEM::setRegion(XY p1, XY p2, int w, int h)
+{
+
+ if(pjsrc == 0) return;
+
+ /*
+ Calculate area of DEM data to be read.
+ */
+
+ XY _p1 = p1;
+ XY _p2 = p2;
+
+ // 1. convert top left and bottom right point into the projection system used by the DEM data
+ pj_transform(pjtar, pjsrc, 1, 0, &_p1.u, &_p1.v, 0);
+ pj_transform(pjtar, pjsrc, 1, 0, &_p2.u, &_p2.v, 0);
+
+ // 2. get floating point offset of topleft corner
+ double xoff1_f = (_p1.u - xref1) / xscale;
+ double yoff1_f = (_p1.v - yref1) / yscale;
+
+ // 3. truncate floating point offset into integer offset
+ int xoff1 = xoff1_f; //qDebug() << "xoff1:" << xoff1 << xoff1_f;
+ int yoff1 = yoff1_f; //qDebug() << "yoff1:" << yoff1 << yoff1_f;
+
+ // 4. get floating point offset of bottom right corner
+ double xoff2_f = (_p2.u - xref1) / xscale;
+ double yoff2_f = (_p2.v - yref1) / yscale;
+
+ // 5. round up (!) floating point offset into integer offset
+ int xoff2 = ceil(xoff2_f); //qDebug() << "xoff2:" << xoff2 << xoff2_f;
+ int yoff2 = ceil(yoff2_f); //qDebug() << "yoff2:" << yoff2 << yoff2_f;
+ int w1 = xoff2 - xoff1; //while((w1 & 0x03) != 0) ++w1;
+ int h1 = yoff2 - yoff1; //qDebug() << "w1:" << w1 << "h1:" << h1;
+
+ // bail out if this is getting too big
+ if(w1 > 10000 || h1 > 10000) return;
+
+ int w2 = w, h2 = h;
+
+ if (w > w1) {
+ w1++; w2= w1;
+ }
+ if (h > h1) {
+ h1++; h2 = h1;
+ }
+
+ region_width = w2;
+ region_height = h2;
+
+ // read 16bit elevation data from GeoTiff
+ region_data = new qint16[w2 * h2];
+ GDALRasterBand * pBand;
+ pBand = dataset->GetRasterBand(1);
+ CPLErr err = pBand->RasterIO(GF_Read, xoff1, yoff1, w1, h1, region_data, w2, h2, GDT_Int16, 0, 0);
+ if(err == CE_Failure) {
+ qDebug() << "faillure" << endl;
+ //delete [] region_data;
+ //FIXME add handle error
+ return;
+ }
+ if ((w2 != w) || (h2 != h)) {
+ // do interpolation if DEM data resolution less than required.
+ qint16 *region_data_tmp = new qint16[w * h];
+ int i, j;
+ double x, y, c, r;
+ w2--; h2--;
+ for (i = 0; i < w; i++) {
+ x = (i / (float) w ) * w2;
+ c = x - (int) x;
+ c = c * abs(xscale);
+
+ for (j = 0; j < h; j++) {
+ y = (j / (float) h ) * h2;
+ r = y - (int) y;
+ r = r * abs(yscale);
+ const weight_t& wt = weights[((int) r * (int)abs(xscale)) + (int) c];
+ region_data_tmp[i + j * w] = wt.c1 * getRegionValue(x, y) + \
+ wt.c2 * getRegionValue(x + 1, y) + \
+ wt.c3 * getRegionValue(x, y + 1) + \
+ wt.c4 * getRegionValue(x + 1, y + 1);
+ }
+ }
+ delete [] region_data;
+ region_data = region_data_tmp;
+ region_width = w;
+ region_height = h;
+ }
+}
+
+float CMapDEM::getRegionValue(int x, int y)
+{
+ float ele = region_data[region_width * y + x];
+ return ele;
+}
+
+void CMapDEM::deleteRegion()
+{
+ delete [] region_data;
+}
float CMapDEM::getElevation(float lon, float lat)
{
diff --git a/src/CMapDEM.h b/src/CMapDEM.h
index 7b91e06..7eaa098 100644
--- a/src/CMapDEM.h
+++ b/src/CMapDEM.h
@@ -49,8 +49,15 @@ class CMapDEM : public IMap
void select(const QRect& rect);
void dimensions(double& lon1, double& lat1, double& lon2, double& lat2);
float getElevation(float lon, float lat);
+ void setRegion(XY p1, XY p2, int w, int h);
+ float getRegionValue(int x, int y);
+ void deleteRegion();
+
private:
+ qint16 * region_data;
+ int region_width;
+ int region_height;
void shading(QImage& img, qint16 * data);
void contour(QImage& img, qint16 * data);
diff --git a/src/IMap.h b/src/IMap.h
index a920cf4..d5a8374 100644
--- a/src/IMap.h
+++ b/src/IMap.h
@@ -199,6 +199,10 @@ class IMap : public QObject
/// select map area for export or further processing
virtual void select(IMapSelection& ms, const QRect& rect){};
+ virtual void setRegion(XY p1, XY p2, int w, int h) {};
+ virtual float getRegionValue(int x, int y) {return 0.0;};
+ virtual void deleteRegion() {};
+
const maptype_e maptype;
signals:
void sigChanged();
--
1.5.6.5
-------------------------------------------------------------------------
This SF.Net email is sponsored by the Moblin Your Move Developer's challenge
Build the coolest Linux based applications with Moblin SDK & win great prizes
Grand prize is a trip for two to an Open Source event anywhere in the world
http://moblin-contest.org/redirect.php?banner_id=100&url=/
_______________________________________________
QLandkarte-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/qlandkarte-users
