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".

Pls, see my patch and may be you will have any comments...
Oliver


>From db1c354e807235036b5c62d5aaebf71a71e449d5 Mon Sep 17 00:00:00 2001
From: Adnrew <[EMAIL PROTECTED]>
Date: Sat, 15 Nov 2008 22:22:46 +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.
---
 src/CMap3DWidget.cpp |   56 ++++++++++++++++++++-------------
 src/CMap3DWidget.h   |    3 ++
 src/CMapDEM.cpp      |   86 +++++++++++++++++++++++++++++++++++++++++++++++++-
 src/CMapDEM.h        |    7 ++++
 src/IMap.h           |    4 ++
 5 files changed, 133 insertions(+), 23 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..7ffae53 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,88 @@ void CMapDEM::dimensions(double& lon1, double& lat1, double& lon2, double& lat2)
 {
 }
 
+void CMapDEM::setRegion(XY p1, XY p2, int w, int h)
+{
+    region_width = w;
+    region_height = 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;
+
+    /*
+        The defined area into DEM data (xoff1,yoff1,xoff2,yoff2) covers a larger or equal
+        area in world coordinates [m] than the current viewport.
+
+        Next the width and height of the area in the DEM data and the corresponding sceen
+        size is calculated.
+    */
+
+    /*
+        Calculate the width and the height of the area. Extend width until it's a multiple of 4.
+        This will be of advantag as QImage will process 32bit alligned bitmaps much faster.
+    */
+    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;
+
+/*    // as the first point off the DEM data will not match exactly the given top left corner
+    // the bitmap has to be drawn with an offset.
+    int pxx = (xoff1_f - xoff1) * xscale / my_xscale;
+    //qDebug() << "pxx:" << pxx << "pxy:" << pxy;
+    int pxy = (yoff1_f - yoff1) * yscale / my_yscale;
+*/
+    // read 16bit elevation data from GeoTiff
+    region_data = new qint16[w * h];
+    GDALRasterBand * pBand;
+    pBand = dataset->GetRasterBand(1);
+    CPLErr err = pBand->RasterIO(GF_Read, xoff1, yoff1, w1, h1, region_data, w, h, GDT_Int16, 0, 0);
+    if(err == CE_Failure) {
+        qDebug() << "faillure" << endl;
+        //delete [] region_data;
+        //FIXME add handle error
+        return;
+    }
+}
+
+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)
 {
@@ -214,7 +297,8 @@ float CMapDEM::getElevation(float lon, float lat)
 
     const weight_t& w = weights[(r * (int)abs(xscale)) + c];
 
-    float ele = w.c1 * e[0] + w.c2 * e[1] + w.c3 * e[2] + w.c4 * e[3];
+//    float ele = w.c1 * e[0] + w.c2 * e[1] + w.c3 * e[2] + w.c4 * e[3];
+    float ele = e[0];
 
     //     qDebug() << c << r << "\t" << w.c1 << e[0] << w.c2 << e[1];
     //     qDebug() << "\t"           << w.c3 << e[2] << w.c4 << e[3];
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

Reply via email to