source: subversion/applications/editors/osm-editor/qt3/SRTMDataGrid.cpp @ 16590

Last change on this file since 16590 was 1158, checked in by nick, 13 years ago

moved current version to qt3 directory in preparation for qt4

File size: 9.1 KB
Line 
1/*
2    Copyright (C) 2005 Nick Whitelegg, Hogweed Software, nick@hogweed.org
3
4    This program is free software; you can redistribute it and/or modify
5    it under the terms of the GNU Lesser General Public License as published by
6    the Free Software Foundation; either version 2 of the License, or
7    (at your option) any later version.
8
9    This program is distributed in the hope that it will be useful,
10    but WITHOUT ANY WARRANTY; without even the implied warranty of
11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12    GNU Lesser General Public License for more details.
13
14    You should have received a copy of the GNU Lesser General Public License
15    along with this program; if not, write to the Free Software
16    Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111 USA
17
18 */
19#include "SRTMDataGrid.h"
20#include "functions.h"
21#include <cmath>
22
23namespace OpenStreetMap
24{
25
26// DataGrid constructor
27// Loads heights from one or more .hgt files.
28// Each rectangle in the input array "rects" defines the rectangle (in sampling
29// point indices) of a grid square. There's normally only one rectangle but
30// may be up to 4 if two latitude/longitude lines pass through the visible area.
31//
32// sampling_pts is a definition of the area in terms of nationally-indexed
33// sampling points; this function fills it in.
34//
35// Returns an array of all the heights indexed nationally.
36
37SRTMDataGrid::SRTMDataGrid(LATLON_TILE **rects,int w, int h, Map& map, int f)
38{
39        // Do each input rectangle
40        int index_w=0, index_h=0,pts_w,pts_h;
41        samplewidth=0; 
42        sampleheight=0;
43
44        this->f=f;
45
46        for(int hcount=0;hcount<h; hcount++)
47        {
48                sampleheight+=
49                                (rects[hcount][0].bottom-rects[hcount][0].top)+1;
50        }
51
52        for(int wcount=0; wcount<w; wcount++)
53        {
54                samplewidth+=
55                                (rects[0][wcount].right-rects[0][wcount].left)+1;
56        }
57
58        points=new SRTM_SAMPLE_POINT [sampleheight*samplewidth];
59
60        for(int hcount=0; hcount<h; hcount++)
61        {
62                index_w = index_h;
63                pts_h = (rects[hcount][0].bottom-rects[hcount][0].top)+1;
64                for(int wcount=0; wcount<w; wcount++)
65                {
66                        pts_w = (rects[hcount][wcount].right-rects[hcount][wcount].left)+1;
67                        doLoad(&rects[hcount][wcount],map,index_w);
68                        index_w += pts_w;
69                }
70                index_h += pts_h*samplewidth; 
71        }
72}
73
74
75
76void SRTMDataGrid::doLoad(LATLON_TILE *rect,Map& map,int index)
77{
78        // Get the .hgt file for the current rectangle
79        char hgtfile[1024];
80        SRTMDataGrid::getHgtFilename(hgtfile,rect->origin);
81        FILE *fp=fopen(hgtfile,"rb");   
82
83        if(fp)
84        {
85
86                int width = (rect->right-rect->left)+1;
87                unsigned char *data = new unsigned char[width*2];
88
89                int datacount; 
90                double h;
91                EarthPoint curLatLon;
92
93                double frac = 0.0008333333333*f;
94                curLatLon.x=rect->origin.x+(((double)rect->left)/1200);
95                curLatLon.y=(rect->origin.y+1)-(((double)rect->top)/1200);
96                double origlong=curLatLon.x;
97                int i;
98
99                for(int row=rect->top; row<=rect->bottom; row++)
100                {
101                        curLatLon.x = origlong; 
102
103                        // 20/02/05 Only do every 'f' rows
104                        if((row-rect->top)%f == 0)
105                        {
106                                fseek(fp,(row*1201+rect->left)*2, SEEK_SET);
107                                fread (data,1,width*2,fp);
108                                datacount=0;
109                                i=0;
110                                for(int pt=row*1201+rect->left;pt<=row*1201+rect->right; pt++)
111                                {
112                                        // 20/02/05 Only do every 'f' columns
113                                        if( (pt-(row*1201+rect->left) )%f == 0) 
114                                        {
115                                                h=
116                                                ( ((double)data[datacount])*256+
117                                                  ((double)data[datacount+1]) ) * 3.28084;
118                                                points[index+i].hgt = 
119                                                        (h>=1 && h<4500) ? h: 1;
120                                                points[index+i].screenPos=
121                                                        (map.isGridRef()) ?
122                                                                map.getScreenPos(ll_to_gr(curLatLon)):
123                                                                map.getScreenPos (curLatLon);
124                                                i++;
125                                        }
126                                        datacount+=2;
127                                        curLatLon.x +=  frac;
128                                }
129                                index+=samplewidth;
130                        }
131                        curLatLon.y -= frac;
132                }
133                delete[] data;
134                fclose(fp);
135        }
136}
137       
138void SRTMDataGrid::getHgtFilename(char *hgtfile,EarthPoint& latlon)
139{
140        sprintf ( hgtfile,"data/N%02d%s%03d.hgt",
141                                                int(latlon.y),(latlon.x<0 ? "W":"E"),
142                                                abs(int(latlon.x)));
143}
144
145
146SRTMDataGrid::~SRTMDataGrid()
147{
148        delete[] points;
149}
150
151
152void SRTMDataGrid::setPoint (int row,int col)
153{
154        pt = row*samplewidth+col;
155        edges[0][0] = pt;
156        edges[0][1] = pt+f;
157        edges[1][0] = pt;
158        edges[1][1] = pt+samplewidth*f;
159        edges[2][0] = pt+f;
160        edges[2][1] = pt+samplewidth*f+f;
161        edges[3][0] = pt+samplewidth*f;
162        edges[3][1] = pt+samplewidth*f+f;
163}
164
165double SRTMDataGrid::startHeight(int interval)
166{
167        double ht = min (
168                                        min(points[pt].hgt, points[pt+f].hgt),
169                                        min(points[pt+samplewidth*f].hgt, 
170                                                 points[pt+samplewidth*f+f].hgt )
171                                        );
172        return ceil(ht/interval) * interval;
173}
174
175double SRTMDataGrid::endHeight(int interval)
176{
177
178        double ht = max (
179                                        max(points[pt].hgt, points[pt+f].hgt),
180                                        max(points[pt+samplewidth*f].hgt, 
181                                                 points[pt+samplewidth*f+f].hgt )
182                                        );
183        return floor(ht/interval) * interval;
184}
185
186// checked : the correct edges are being considered each time.
187void SRTMDataGrid::getLine(LINE *lines,int *n_lines, int ht)
188{
189        int go = 0;
190        int prevedges[2];
191        LINE line;
192        double eAh0, eAh1, eBh0, eBh1; 
193        int eAp0x, eAp0y, eAp1x, eAp1y, eBp0x, eBp0y, eBp1x, eBp1y;
194
195        // See addContour()
196        bool two_contours = (
197                        (points[pt].hgt<ht && points[pt+f].hgt>ht && 
198                         points[pt+samplewidth*f+f].hgt<ht && 
199                         points[pt+samplewidth*f].hgt > ht) 
200                         ||
201                        (points[pt].hgt>ht && points[pt+f].hgt<ht && 
202                         points[pt+samplewidth*f+f].hgt>ht && 
203                         points[pt+samplewidth*f].hgt < ht)
204                                ) ; 
205
206        for(int edge=0; edge<3; edge++)
207        {
208                if(between(ht,points[edges[edge][0]].hgt,points[edges[edge][1]].hgt))
209                {
210                        for(int edge2=edge+1; edge2<4; edge2++)
211                        {
212                                if(between
213                                        (ht,points[edges[edge2][0]].hgt,
214                                                 points[edges[edge2][1]].hgt))
215                                {
216                                        eAh0 = points[edges[edge][0]].hgt;
217                                        eAh1 = points[edges[edge][1]].hgt;
218                                        eBh0 = points[edges[edge2][0]].hgt;
219                                        eBh1 = points[edges[edge2][1]].hgt;
220
221                                        eAp0x = points[edges[edge][0]].screenPos.x;
222                                        eAp0y = points[edges[edge][0]].screenPos.y;
223                                        eAp1x = points[edges[edge][1]].screenPos.x;
224                                        eAp1y = points[edges[edge][1]].screenPos.y;
225
226                                        eBp0x = points[edges[edge2][0]].screenPos.x;
227                                        eBp0y = points[edges[edge2][0]].screenPos.y;
228                                        eBp1x = points[edges[edge2][1]].screenPos.x;
229                                        eBp1y = points[edges[edge2][1]].screenPos.y;
230
231                                        // We draw a line.
232                                        line.p[0].x =
233                                                        eAp0x + ( ((ht-eAh0) / (eAh1-eAh0))     
234                                                         * (eAp1x-eAp0x) );
235
236                                        line.p[0].y =
237                                                        eAp0y + ( ((ht-eAh0) / (eAh1-eAh0))     
238                                                         * (eAp1y-eAp0y) );
239
240                                        line.p[1].x =
241                                                        eBp0x + ( ((ht-eBh0) / (eBh1-eBh0))     
242                                                         * (eBp1x-eBp0x) );
243
244                                        line.p[1].y =
245                                                        eBp0y + ( ((ht-eBh0) / (eBh1-eBh0))     
246                                                         * (eBp1y-eBp0y) );
247
248                                        if(two_contours==false) 
249                                        {
250                                                lines[0] = line;
251                                                *n_lines=1;
252                                        }
253                                        else
254                                        {
255                                                SRTMDataGrid::addContour
256                                                                (lines,line,edge,edge2,prevedges,&go,
257                                                                 n_lines);
258                                        }
259                                }
260                        }
261                }
262        }
263}
264
265// This function will be called when the special case of two opposite corners
266// having a height above a contour, and the other two below, occurs. In this
267// case - and this case only - two contours of a given height will be drawn
268// through a quadrangle. This special case confuses the standard algorithm
269// no end :-) Thus, we need to ensure that the two contours are drawn
270// on the opposite side of the quadrangle (any other combination wouldn't
271// make sense), and this function does that. Which two opposite sides
272// doesn't actually matter; in fact, it's impossible to tell which two would
273// be correct.
274void SRTMDataGrid::addContour
275        (LINE *lines,LINE line,int edge1,int edge2,int *prevedges,int *go,
276        int *n_lines)
277{
278        /* Edge numbering assumed to be:
279
280             0
281           +---+
282          1|   |2
283           +---+
284             3
285
286        Don't draw contours through opposite edges. These don't apply in this
287        special case.
288
289        */     
290
291        if(!((edge1==0 && edge2==3) || (edge1==1 && edge2==2)))
292        {
293                if(*go==0)
294                {
295                        prevedges[0] = edge1;
296                        prevedges[1] = edge2;
297                        *go = 1;
298                        *n_lines = 1;
299                        lines[0] = line;
300                }
301                else if( (edge1==2 && edge2==3 && prevedges[0]==0 && prevedges[1]==1)
302                         || (edge1==1 && edge2==3 && prevedges[0]==0 && prevedges[1]==2)
303                          && *go==1)
304                {
305                        lines[1] = line;
306                        *n_lines = 2;
307                }
308        }
309}       
310
311// finds the minimum distance between a point and the list of previous points
312int SRTMDataGrid::hgtptDistance(vector<int>& prevs)
313{
314        int dist=sqrt(samplewidth*sampleheight), result; 
315        int dx, dy;
316        for(int count=0; count<prevs.size(); count++)
317        {
318                dx = prevs[count]%samplewidth - pt%samplewidth;
319                dy = prevs[count]/sampleheight - pt/sampleheight;
320                result= sqrt(dx*dx + dy*dy);
321                dist = (result<dist) ? result:dist;
322        }
323        return dist;
324}
325
326Colour SRTMDataGrid::getHeightShading (double shadingres)
327{
328        double est_ht = floor (((points[pt].hgt+points[pt+samplewidth*f+f].hgt)/2+
329                         (points[pt+f].hgt+points[pt+samplewidth*f].hgt)/2 
330                   ) / 2);
331        return doGetHeightShading(est_ht,shadingres);
332}
333
334Colour SRTMDataGrid::doGetHeightShading(double ht,double shadingres)
335{
336        Colour colour;
337        ht=round(ht/shadingres);
338        double a=1250/shadingres,
339        b=250/shadingres,
340        c=2000/shadingres,
341        d=500/shadingres,
342        e=1750/shadingres;
343       
344        colour.r = (ht<a) ? 191+(ht*(64/a)) : 
345        ((ht<c) ? 255-((ht-(1250/shadingres))*(16/b)) : 
346                                  255-((ht-d)*(16/d)) );
347
348        colour.g=(ht<c) ? 255-(ht*(16/b)) : 255-((ht+c)*(16/d));
349
350        colour.b = (ht<e) ? 191-(ht*(16/b)) :
351        ((ht<c) ? 95+((ht-e)*(16/b)) :
352                                  95+((ht-(1500/shadingres))*(16/d)) );
353        return colour;
354}
355
356}
Note: See TracBrowser for help on using the repository browser.