1 /******************************************************************************
2 *
3 * Project: OpenCPN
4 * Purpose: MBTiles chart type support
5 * Author: David Register
6 *
7 ***************************************************************************
8 * Copyright (C) 2018 by David S. Register *
9 * *
10 * This program is free software; you can redistribute it and/or modify *
11 * it under the terms of the GNU General Public License as published by *
12 * the Free Software Foundation; either version 2 of the License, or *
13 * (at your option) any later version. *
14 * *
15 * This program is distributed in the hope that it will be useful, *
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
18 * GNU General Public License for more details. *
19 * *
20 * You should have received a copy of the GNU General Public License *
21 * along with this program; if not, write to the *
22 * Free Software Foundation, Inc., *
23 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. *
24 ***************************************************************************
25 *
26 */
27
28
29 // ============================================================================
30 // declarations
31 // ============================================================================
32
33
34 // ----------------------------------------------------------------------------
35 // headers
36 // ----------------------------------------------------------------------------
37
38 // For compilers that support precompilation, includes "wx.h".
39 #include "wx/wxprec.h"
40
41 #ifndef WX_PRECOMP
42 #include "wx/wx.h"
43 #endif //precompiled headers
44
45
46 // Why are these not in wx/prec.h?
47 #include "wx/dir.h"
48 #include "wx/stream.h"
49 #include "wx/wfstream.h"
50 #include "wx/tokenzr.h"
51 #include "wx/filename.h"
52 #include <wx/image.h>
53 #include <wx/fileconf.h>
54 #include <wx/mstream.h>
55 #include <sys/stat.h>
56 #include <sstream>
57 #include <map>
58 #include <unordered_map>
59
60 #include <sqlite3.h> //We need some defines
61 #include <SQLiteCpp/SQLiteCpp.h>
62
63 #include "mbtiles.h"
64 #include "chart1.h"
65 #include "chcanv.h"
66 #include "glChartCanvas.h"
67
68 // Missing from MSW include files
69 #ifdef _MSC_VER
70 typedef __int32 int32_t;
71 typedef unsigned __int32 uint32_t;
72 typedef __int64 int64_t;
73 typedef unsigned __int64 uint64_t;
74 #endif
75
76 // ----------------------------------------------------------------------------
77 // Random Prototypes
78 // ----------------------------------------------------------------------------
79
80 #if !defined(NAN)
81 static const long long lNaN = 0xfff8000000000000;
82 #define NAN (*(double*)&lNaN)
83 #endif
84
85 #ifdef OCPN_USE_CONFIG
86 class MyConfig;
87 extern MyConfig *pConfig;
88 #endif
89
90 #define LON_UNDEF NAN
91 #define LAT_UNDEF NAN
92
93 // The OpenStreetMaps zommlevel translation tables
94 // https://wiki.openstreetmap.org/wiki/Zoom_levels
95
96 /*Level Degree Area m / pixel ~Scale # Tiles
97 0 360 whole world 156,412 1:500 million 1
98 1 180 78,206 1:250 million 4
99 2 90 39,103 1:150 million 16
100 3 45 19,551 1:70 million 64
101 4 22.5 9,776 1:35 million 256
102 5 11.25 4,888 1:15 million 1,024
103 6 5.625 2,444 1:10 million 4,096
104 7 2.813 1,222 1:4 million 16,384
105 8 1.406 610.984 1:2 million 65,536
106 9 0.703 wide area 305.492 1:1 million 262,144
107 10 0.352 152.746 1:500,000 1,048,576
108 11 0.176 area 76.373 1:250,000 4,194,304
109 12 0.088 38.187 1:150,000 16,777,216
110 13 0.044 village or town 19.093 1:70,000 67,108,864
111 14 0.022 9.547 1:35,000 268,435,456
112 15 0.011 4.773 1:15,000 1,073,741,824
113 16 0.005 small road 2.387 1:8,000 4,294,967,296
114 17 0.003 1.193 1:4,000 17,179,869,184
115 18 0.001 0.596 1:2,000 68,719,476,736
116 19 0.0005 0.298 1:1,000 274,877,906,944
117 */
118
119 // A "nominal" scale value, by zoom factor. Estimated at equator, with monitor pixel size of 0.3mm
120 static const double OSM_zoomScale[] = { 5e8,
121 2.5e8,
122 1.5e8,
123 7.0e7,
124 3.5e7,
125 1.5e7,
126 1.0e7,
127 4.0e6,
128 2.0e6,
129 1.0e6,
130 5.0e5,
131 2.5e5,
132 1.5e5,
133 7.0e4,
134 3.5e4,
135 1.5e4,
136 8.0e3,
137 4.0e3,
138 2.0e3,
139 1.0e3};
140
141 // Meters per pixel, by zoom factor
142 static const double OSM_zoomMPP[] = { 156412,
143 78206,
144 39103,
145 19551,
146 9776,
147 4888,
148 2444,
149 1222,
150 610,984,
151 305.492,
152 152.746,
153 76.373,
154 38.187,
155 19.093,
156 9.547,
157 4.773,
158 2.387,
159 1.193,
160 0.596,
161 0.298 };
162
163
164 static const double eps = 6e-6; // about 1cm on earth's surface at equator
165 extern MyFrame *gFrame;
166
167 #if defined( __UNIX__ ) && !defined(__WXOSX__) // high resolution stopwatch for profiling
168 class OCPNStopWatch
169 {
170 public:
OCPNStopWatch()171 OCPNStopWatch() { Reset(); }
Reset()172 void Reset() { clock_gettime(CLOCK_REALTIME, &tp); }
173
GetTime()174 double GetTime() {
175 timespec tp_end;
176 clock_gettime(CLOCK_REALTIME, &tp_end);
177 return (tp_end.tv_sec - tp.tv_sec) * 1.e3 + (tp_end.tv_nsec - tp.tv_nsec) / 1.e6;
178 }
179
180 private:
181 timespec tp;
182 };
183 #endif
184
185 // *********************************************
186 // Utility Functions
187 // *********************************************
188
189 // https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#C.2FC.2B.2B
long2tilex(double lon,int z)190 static int long2tilex(double lon, int z)
191 {
192 if(lon < -180)
193 lon += 360;
194
195 return (int)(floor((lon + 180.0) / 360.0 * pow(2.0, z)));
196 }
197
lat2tiley(double lat,int z)198 static int lat2tiley(double lat, int z)
199 {
200 int y = (int)(floor((1.0 - log( tan(lat * M_PI/180.0) + 1.0 / cos(lat * M_PI/180.0)) / M_PI) / 2.0 * pow(2.0, z)));
201 int ymax = 1 << z;
202 y = ymax - y - 1;
203 return y;
204 }
205
tilex2long(int x,int z)206 static double tilex2long(int x, int z)
207 {
208 return x / pow(2.0, z) * 360.0 - 180;
209 }
210
tiley2lat(int y,int z)211 static double tiley2lat(int y, int z)
212 {
213 double n = pow(2.0,z);
214 int ymax = 1 << z;
215 y = ymax - y - 1;
216 double latRad = atan(sinh(M_PI*(1-(2*y/n))));
217 return 180.0 / M_PI * latRad;
218 }
219
220
221
222 // ----------------------------------------------------------------------------
223 // private classes
224 // ----------------------------------------------------------------------------
225
226
227 // Per tile descriptor
228 class mbTileDescriptor
229 {
230 public:
mbTileDescriptor()231 mbTileDescriptor() { glTextureName = 0; m_bAvailable = false; m_bgeomSet = false;}
232
~mbTileDescriptor()233 virtual ~mbTileDescriptor() { }
234
235 int tile_x, tile_y;
236 int m_zoomLevel;
237 float latmin, lonmin, latmax, lonmax;
238 LLBBox box;
239
240 GLuint glTextureName;
241 bool m_bAvailable;
242 bool m_bgeomSet;
243
244 };
245
246 // Per zoomlevel descriptor of tile array for that zoomlevel
247 class mbTileZoomDescriptor
248 {
249 public:
mbTileZoomDescriptor()250 mbTileZoomDescriptor(){}
~mbTileZoomDescriptor()251 virtual ~mbTileZoomDescriptor(){}
252
253 int tile_x_min, tile_x_max;
254 int tile_y_min, tile_y_max;
255
256 int nx_tile, ny_tile;
257
258 //std::map<unsigned int, mbTileDescriptor *> tileMap;
259 std::unordered_map<unsigned int, mbTileDescriptor *> tileMap;
260 };
261
262
263
264
265
266
267 // ============================================================================
268 // ChartMBTiles implementation
269 // ============================================================================
270
271
ChartMBTiles()272 ChartMBTiles::ChartMBTiles()
273 {
274 // Init some private data
275 m_ChartFamily = CHART_FAMILY_RASTER;
276 m_ChartType = CHART_TYPE_MBTILES;
277
278 m_Chart_Skew = 0.0;
279
280 m_datum_str = _T("WGS84"); // assume until proven otherwise
281 m_projection = PROJECTION_WEB_MERCATOR;
282 m_imageType = wxBITMAP_TYPE_ANY;
283
284 m_b_cdebug = 0;
285
286 m_minZoom = 0;
287 m_maxZoom = 19;
288
289 m_nNoCOVREntries = 0;
290 m_nCOVREntries = 0;
291 m_pCOVRTablePoints = NULL;
292 m_pCOVRTable = NULL;
293 m_pNoCOVRTablePoints = NULL;
294 m_pNoCOVRTable = NULL;
295 m_tileArray = NULL;
296
297 m_LonMin = LON_UNDEF;
298 m_LonMax = LON_UNDEF;
299 m_LatMin = LAT_UNDEF;
300 m_LatMax = LAT_UNDEF;
301
302 #ifdef OCPN_USE_CONFIG
303 wxFileConfig *pfc = (wxFileConfig *)pConfig;
304 pfc->SetPath ( _T ( "/Settings" ) );
305 pfc->Read ( _T ( "DebugMBTiles" ), &m_b_cdebug, 0 );
306 #endif
307 m_pDB = NULL;
308
309 }
310
~ChartMBTiles()311 ChartMBTiles::~ChartMBTiles()
312 {
313 FlushTiles();
314 if(m_pDB){
315 delete m_pDB;
316 }
317 }
318
319 //-------------------------------------------------------------------------------------------------
320 // Get the Chart thumbnail data structure
321 // Creating the thumbnail bitmap as required
322 //-------------------------------------------------------------------------------------------------
323
GetThumbData()324 ThumbData *ChartMBTiles::GetThumbData()
325 {
326 return NULL;
327 }
328
GetThumbData(int tnx,int tny,float lat,float lon)329 ThumbData *ChartMBTiles::GetThumbData(int tnx, int tny, float lat, float lon)
330 {
331 return NULL;
332 }
333
UpdateThumbData(double lat,double lon)334 bool ChartMBTiles::UpdateThumbData(double lat, double lon)
335 {
336 return true;
337 }
338
AdjustVP(ViewPort & vp_last,ViewPort & vp_proposed)339 bool ChartMBTiles::AdjustVP(ViewPort &vp_last, ViewPort &vp_proposed)
340 {
341 return true;
342 }
343
344
345
346 // Report recommended minimum and maximum scale values for which use of this chart is valid
347
GetNormalScaleMin(double canvas_scale_factor,bool b_allow_overzoom)348 double ChartMBTiles::GetNormalScaleMin(double canvas_scale_factor, bool b_allow_overzoom)
349 {
350 // if(b_allow_overzoom)
351 return (canvas_scale_factor / m_ppm_avg) / 132; // allow wide range overzoom overscale
352 // else
353 // return (canvas_scale_factor / m_ppm_avg) / 2; // don't suggest too much overscale
354
355 }
356
GetNormalScaleMax(double canvas_scale_factor,int canvas_width)357 double ChartMBTiles::GetNormalScaleMax(double canvas_scale_factor, int canvas_width)
358 {
359 return (canvas_scale_factor / m_ppm_avg) * 40.0; // excessive underscale is slow, and unreadable
360 }
361
362
GetNearestPreferredScalePPM(double target_scale_ppm)363 double ChartMBTiles::GetNearestPreferredScalePPM(double target_scale_ppm)
364 {
365 return target_scale_ppm;
366 }
367
368 //Checks/corrects/completes the initialization based on real data from the tiles table
InitFromTiles(const wxString & name)369 void ChartMBTiles::InitFromTiles( const wxString& name )
370 {
371 try
372 {
373 // Open the MBTiles database file
374 const char *name_UTF8 = "";
375 wxCharBuffer utf8CB = name.ToUTF8(); // the UTF-8 buffer
376 if ( utf8CB.data() )
377 name_UTF8 = utf8CB.data();
378
379 SQLite::Database db(name_UTF8);
380
381 // Check if tiles with advertised min and max zoom level really exist, or correct the defaults
382 // We can't blindly use what we find though - the DB often contains empty cells at very low zoom levels, so if we have some info from metadata, we will use that if more conservative...
383 SQLite::Statement query(db, "SELECT min(zoom_level) AS min_zoom, max(zoom_level) AS max_zoom FROM tiles");
384 while (query.executeStep())
385 {
386 const char* colMinZoom = query.getColumn(0);
387 const char* colMaxZoom = query.getColumn(1);
388
389 int min_zoom, max_zoom;
390 sscanf( colMinZoom, "%i", &min_zoom );
391 m_minZoom = wxMax(m_minZoom, min_zoom);
392 sscanf( colMaxZoom, "%i", &max_zoom );
393 m_maxZoom = wxMin(m_maxZoom, max_zoom);
394 if(m_minZoom > m_maxZoom) {
395 //We are looking at total nonsense with wrong metatadata and actual tile coverage out of it, better use what's really in the data to be able to show at least something
396 m_minZoom = min_zoom;
397 m_maxZoom = max_zoom;
398 }
399 }
400
401 // std::cout << name.c_str() << " zoom_min: " << m_minZoom << " zoom_max: " << m_maxZoom << std::endl;
402
403 // Traversing the entire tile table can be expensive....
404 // Use declared bounds if present.
405
406 if(!std::isnan(m_LatMin) && !std::isnan(m_LatMax) && !std::isnan(m_LonMin) && !std::isnan(m_LonMax) )
407 return;
408
409 // Try to guess the coverage extents from the tiles. This will be hard to get right -
410 //the finest resolution likely does not cover the whole area, while the lowest resolution tiles probably contain a lot of theoretical space
411 //which actually is not covered. And some resolutions may be actually missing... What do we use?
412 // If we have the metadata and it is not completely off, we should probably prefer it.
413 SQLite::Statement query1(db, wxString::Format("SELECT min(tile_row) AS min_row, max(tile_row) as max_row, min(tile_column) as min_column, max(tile_column) as max_column, count(*) as cnt, zoom_level FROM tiles WHERE zoom_level >= %d AND zoom_level <= %d GROUP BY zoom_level ORDER BY zoom_level ASC", m_minZoom, m_maxZoom).c_str());
414 float minLat = 999., maxLat = -999.0, minLon = 999., maxLon = -999.0;
415 while (query1.executeStep())
416 {
417 const char* colMinRow = query1.getColumn(0);
418 const char* colMaxRow = query1.getColumn(1);
419 const char* colMinCol = query1.getColumn(2);
420 const char* colMaxCol = query1.getColumn(3);
421 const char* colCnt = query1.getColumn(4);
422 const char* colZoom = query1.getColumn(5);
423
424 int minRow, maxRow, minCol, maxCol, cnt, zoom;
425 sscanf( colMinRow, "%i", &minRow );
426 sscanf( colMaxRow, "%i", &maxRow );
427 sscanf( colMinCol, "%i", &minCol );
428 sscanf( colMaxCol, "%i", &maxCol );
429 sscanf( colMinRow, "%i", &minRow );
430 sscanf( colMaxRow, "%i", &maxRow );
431 sscanf( colCnt, "%i", &cnt );
432 sscanf( colZoom, "%i", &zoom );
433
434 // Let's try to use the simplest possible algo and just look for the zoom level with largest extent (Which probably be the one with lowest resolution?)...
435 minLat = wxMin(minLat, tiley2lat(minRow, zoom));
436 maxLat = wxMax(maxLat, tiley2lat(maxRow - 1, zoom));
437 minLon = wxMin(minLon, tilex2long(minCol, zoom));
438 maxLon = wxMax(maxLon, tilex2long(maxCol + 1, zoom));
439 //std::cout << "Zoom: " << zoom << " minlat: " << tiley2lat(minRow, zoom) << " maxlat: " << tiley2lat(maxRow - 1, zoom) << " minlon: " << tilex2long(minCol, zoom) << " maxlon: " << tilex2long(maxCol + 1, zoom) << std::endl;
440 }
441
442 // ... and use what we found only in case we miss some of the values from metadata...
443 if(std::isnan(m_LatMin))
444 m_LatMin = minLat;
445 if(std::isnan(m_LatMax))
446 m_LatMax = maxLat;
447 if(std::isnan(m_LonMin))
448 m_LonMin = minLon;
449 if(std::isnan(m_LonMax))
450 m_LonMax = maxLon;
451 }
452 catch (std::exception& e)
453 {
454 const char *t = e.what();
455 wxLogMessage("mbtiles exception: %s", e.what());
456 }
457 }
458
Init(const wxString & name,ChartInitFlag init_flags)459 InitReturn ChartMBTiles::Init( const wxString& name, ChartInitFlag init_flags )
460 {
461 m_global_color_scheme = GLOBAL_COLOR_SCHEME_RGB;
462
463 m_FullPath = name;
464 m_Description = m_FullPath;
465
466 try
467 {
468 // Open the MBTiles database file
469 const char *name_UTF8 = "";
470 wxCharBuffer utf8CB = name.ToUTF8(); // the UTF-8 buffer
471 if ( utf8CB.data() )
472 name_UTF8 = utf8CB.data();
473
474 SQLite::Database db(name_UTF8);
475
476 // Compile a SQL query, getting everything from the "metadata" table
477 SQLite::Statement query(db, "SELECT * FROM metadata ");
478
479 // Loop to execute the query step by step, to get rows of result
480 while (query.executeStep())
481 {
482 const char* colName = query.getColumn(0);
483 const char* colValue = query.getColumn(1);
484
485 //Get the geometric extent of the data
486 if(!strncmp(colName, "bounds", 6)){
487 float lon1, lat1, lon2, lat2;
488 sscanf( colValue, "%g,%g,%g,%g", &lon1, &lat1, &lon2, &lat2 );
489
490 // There is some confusion over the layout of this field...
491 m_LatMax = wxMax(lat1, lat2);
492 m_LatMin = wxMin(lat1, lat2);
493 m_LonMax = wxMax(lon1, lon2);
494 m_LonMin = wxMin(lon1, lon2);
495
496 }
497
498 // Not very interesting as it may be wrong, we better find out from first loaded tile and adjust m_imageType accordingly
499 //else if(!strncmp(colName, "format", 6) ){
500 // m_bPNG = !strncmp(colValue, "png", 3);
501 //}
502
503 //Get the min and max zoom values present in the db
504 else if(!strncmp(colName, "minzoom", 7)){
505 sscanf( colValue, "%i", &m_minZoom );
506 }
507 else if(!strncmp(colName, "maxzoom", 7)){
508 sscanf( colValue, "%i", &m_maxZoom );
509 }
510
511 else if(!strncmp(colName, "description", 11)){
512 m_Description = wxString(colValue, wxConvUTF8);
513 }
514 else if(!strncmp(colName, "name", 11)){
515 m_Name = wxString(colValue, wxConvUTF8);
516 }
517 else if(!strncmp(colName, "type", 11)){
518 m_Type = wxString(colValue, wxConvUTF8).Upper().IsSameAs("OVERLAY") ? MBTilesType::OVERLAY : MBTilesType::BASE;
519 }
520 else if(!strncmp(colName, "scheme", 11)){
521 m_Scheme = wxString(colValue, wxConvUTF8).Upper().IsSameAs("XYZ") ? MBTilesScheme::XYZ : MBTilesScheme::TMS;
522 }
523 }
524 }
525 catch (std::exception& e)
526 {
527 const char *t = e.what();
528 wxLogMessage("mbtiles exception: %s", e.what());
529 return INIT_FAIL_REMOVE;
530 }
531
532 // Fix the missing/wrong metadata values
533 InitFromTiles(name);
534
535
536 // set the chart scale parameters based on the max zoom factor
537 m_ppm_avg = 1.0 / OSM_zoomMPP[m_minZoom];
538 m_Chart_Scale = OSM_zoomScale[m_maxZoom];
539
540
541 PrepareTiles(); // Initialize the tile data structures
542
543
544 LLRegion covrRegion;
545
546 LLBBox extentBox;
547 extentBox.Set(m_LatMin, m_LonMin, m_LatMax, m_LonMax);
548
549 const char *name_UTF8 = "";
550 wxCharBuffer utf8CB = name.ToUTF8(); // the UTF-8 buffer
551 if ( utf8CB.data() )
552 name_UTF8 = utf8CB.data();
553
554 SQLite::Database db(name_UTF8);
555
556 int zoomFactor = m_minZoom;
557 int minRegionZoom = -1;
558 bool covr_populated = false;
559
560 m_nTiles = 0;
561 while( (zoomFactor <= m_maxZoom) && (minRegionZoom < 0) )
562 {
563 LLRegion covrRegionZoom;
564 wxRegion regionZoom;
565 char qrs[100];
566
567 // Protect against trying to create the exact coverage for the brutal large scale layers contianing tens of thousand tiles.
568 sprintf(qrs, "select count(*) from tiles where zoom_level = %d ", zoomFactor);
569 SQLite::Statement query_size(db, qrs);
570
571 if(query_size.executeStep())
572 {
573 const char* colValue = query_size.getColumn(0);
574 int tile_at_zoom = atoi(colValue);
575 m_nTiles += tile_at_zoom;
576
577 if (tile_at_zoom > 1000) {
578 zoomFactor++;
579 if(!covr_populated) {
580 covr_populated = true;
581 covrRegion = extentBox;
582 }
583 continue;
584 }
585 }
586
587 // query the database
588 sprintf(qrs, "select tile_column, tile_row from tiles where zoom_level = %d ", zoomFactor);
589
590
591 // Compile a SQL query, getting the specific data
592 SQLite::Statement query(db, qrs);
593 covr_populated = true;
594
595 while (query.executeStep())
596 {
597 const char* colValue = query.getColumn(0);
598 const char* c2 = query.getColumn(1);
599 int tile_x_found = atoi(colValue); // tile_x
600 int tile_y_found = atoi(c2); // tile_y
601
602 regionZoom.Union(tile_x_found, tile_y_found-1, 1, 1);
603
604 } // inner while
605
606 wxRegionIterator upd( regionZoom ); // get the rect list
607 double eps_factor = eps * 100; // roughly 1 m
608
609 while( upd ) {
610 wxRect rect = upd.GetRect();
611
612 double lonmin = round(tilex2long(rect.x, zoomFactor)/eps_factor)*eps_factor;
613 double lonmax = round(tilex2long(rect.x + rect.width, zoomFactor)/eps_factor)*eps_factor;
614 double latmin = round(tiley2lat(rect.y, zoomFactor)/eps_factor)*eps_factor;
615 double latmax = round(tiley2lat(rect.y + rect.height, zoomFactor)/eps_factor)*eps_factor;
616
617 LLBBox box;
618 box.Set(latmin, lonmin, latmax, lonmax);
619
620 LLRegion tileRegion(box);
621 //if(i <= 1)
622 covrRegionZoom.Union(tileRegion);
623
624 upd++;
625 minRegionZoom = zoomFactor; // We take the first populated (lowest) zoom level
626 // region as the final chart region
627 }
628
629 covrRegion.Union(covrRegionZoom);
630
631 zoomFactor++;
632
633 } // while
634
635
636 // The coverage region must be reduced if necessary to include only the db specified bounds.
637 covrRegion.Intersect(extentBox);
638
639 m_minZoomRegion = covrRegion;
640
641 // Populate M_COVR entries for the OCPN chart database
642 if(covrRegion.contours.size()){ // Check for no intersection caused by ??
643 m_nCOVREntries = covrRegion.contours.size();
644 m_pCOVRTablePoints = (int *)malloc(m_nCOVREntries * sizeof(int));
645 m_pCOVRTable = (float **)malloc(m_nCOVREntries * sizeof(float *));
646 std::list<poly_contour>::iterator it = covrRegion.contours.begin();
647 for(int i=0; i<m_nCOVREntries; i++) {
648 m_pCOVRTablePoints[i] = it->size();
649 m_pCOVRTable[i] = (float *)malloc(m_pCOVRTablePoints[i] * 2 * sizeof(float));
650 std::list<contour_pt>::iterator jt = it->begin();
651 for(int j=0; j<m_pCOVRTablePoints[i]; j++) {
652 m_pCOVRTable[i][2*j+0] = jt->y;
653 m_pCOVRTable[i][2*j+1] = jt->x;
654 jt++;
655 }
656 it++;
657 }
658 }
659
660 if(init_flags == HEADER_ONLY)
661 return INIT_OK;
662
663 InitReturn pi_ret = PostInit();
664 if( pi_ret != INIT_OK)
665 return pi_ret;
666 else
667 return INIT_OK;
668
669 }
670
PreInit(const wxString & name,ChartInitFlag init_flags,ColorScheme cs)671 InitReturn ChartMBTiles::PreInit( const wxString& name, ChartInitFlag init_flags, ColorScheme cs )
672 {
673 m_global_color_scheme = cs;
674 return INIT_OK;
675 }
676
677
PostInit(void)678 InitReturn ChartMBTiles::PostInit(void)
679 {
680 // Create the persistent MBTiles database file
681 const char *name_UTF8 = "";
682 wxCharBuffer utf8CB = m_FullPath.ToUTF8(); // the UTF-8 buffer
683 if ( utf8CB.data() )
684 name_UTF8 = utf8CB.data();
685
686 m_pDB = new SQLite::Database(name_UTF8);
687 m_pDB->exec("PRAGMA locking_mode=EXCLUSIVE");
688 m_pDB->exec("PRAGMA cache_size=-50000");
689
690 bReadyToRender = true;
691 return INIT_OK;
692 }
693
PrepareTiles()694 void ChartMBTiles::PrepareTiles()
695 {
696 //OCPNStopWatch sw;
697 m_tileArray = new mbTileZoomDescriptor* [(m_maxZoom - m_minZoom) + 1];
698
699 for(int i=0 ; i < (m_maxZoom - m_minZoom) + 1 ; i++){
700 PrepareTilesForZoom(m_minZoom + i, (i==0)); // Preset the geometry only on the minZoom tiles
701 }
702 //printf("PrepareTiles time: %f\n", sw.GetTime());
703
704 }
705
FlushTiles()706 void ChartMBTiles::FlushTiles()
707 {
708 if(!bReadyToRender || m_tileArray == nullptr)
709 return;
710 for(int iz=0 ; iz < (m_maxZoom - m_minZoom) + 1 ; iz++){
711 mbTileZoomDescriptor *tzd = m_tileArray[iz];
712
713 for (auto const &it : tzd->tileMap)
714 {
715 mbTileDescriptor *tile = it.second;
716 if( tile ){
717 if (tile->glTextureName > 0)
718 glDeleteTextures(1, &tile->glTextureName);
719 delete tile;
720 }
721 }
722 delete tzd;
723 }
724 }
725
FlushTextures()726 void ChartMBTiles::FlushTextures()
727 {
728 if (m_tileArray == nullptr) {
729 return;
730 }
731 for(int iz=0 ; iz < (m_maxZoom - m_minZoom) + 1 ; iz++){
732 mbTileZoomDescriptor *tzd = m_tileArray[iz];
733
734 for (auto const &it : tzd->tileMap)
735 {
736 mbTileDescriptor *tile = it.second;
737 if( tile && tile->glTextureName > 0){
738 glDeleteTextures(1, &tile->glTextureName);
739 tile->glTextureName = 0;
740 }
741 }
742 }
743 }
744
745
PrepareTilesForZoom(int zoomFactor,bool bset_geom)746 void ChartMBTiles::PrepareTilesForZoom(int zoomFactor, bool bset_geom)
747 {
748 mbTileZoomDescriptor *tzd = new mbTileZoomDescriptor;
749
750 m_tileArray[zoomFactor - m_minZoom] = tzd;
751
752 // Calculate the tile counts in x and y, based on zoomfactor and chart extents
753 tzd->tile_x_min = long2tilex(m_LonMin + eps, zoomFactor);
754 tzd->tile_x_max = long2tilex(m_LonMax- eps, zoomFactor);
755 tzd->tile_y_min = lat2tiley(m_LatMin + eps, zoomFactor);
756 tzd->tile_y_max = lat2tiley(m_LatMax - eps, zoomFactor);
757
758 tzd->nx_tile = abs(tzd->tile_x_max - tzd->tile_x_min) + 1;
759 tzd->ny_tile = tzd->tile_y_max - tzd->tile_y_min + 1;
760
761 return;
762 }
763
764
765
GetChartExtent(Extent * pext)766 bool ChartMBTiles::GetChartExtent(Extent *pext)
767 {
768 pext->NLAT = m_LatMax;
769 pext->SLAT = m_LatMin;
770 pext->ELON = m_LonMax;
771 pext->WLON = m_LonMin;
772
773 return true;
774 }
775
776
777
SetColorScheme(ColorScheme cs,bool bApplyImmediate)778 void ChartMBTiles::SetColorScheme(ColorScheme cs, bool bApplyImmediate)
779 {
780 if(m_global_color_scheme != cs){
781 m_global_color_scheme = cs;
782 FlushTextures();
783 }
784 }
785
786
787
788
789
GetValidCanvasRegion(const ViewPort & VPoint,OCPNRegion * pValidRegion)790 void ChartMBTiles::GetValidCanvasRegion(const ViewPort& VPoint, OCPNRegion *pValidRegion)
791 {
792
793 pValidRegion->Clear();
794 pValidRegion->Union(0, 0, VPoint.pix_width, VPoint.pix_height);
795 return;
796 }
797
GetValidRegion()798 LLRegion ChartMBTiles::GetValidRegion()
799 {
800 return m_minZoomRegion;
801 }
802
803
RenderViewOnDC(wxMemoryDC & dc,const ViewPort & VPoint)804 bool ChartMBTiles::RenderViewOnDC(wxMemoryDC& dc, const ViewPort& VPoint)
805 {
806 return true;
807 }
808
getTileTexture(mbTileDescriptor * tile)809 bool ChartMBTiles::getTileTexture( mbTileDescriptor *tile)
810 {
811 if(!m_pDB)
812 return false;
813
814 // Is the texture ready?
815 if(tile->glTextureName > 0){
816 glBindTexture( GL_TEXTURE_2D, tile->glTextureName );
817
818 return true;
819 }
820 else{
821 if(!tile->m_bAvailable)
822 return false;
823 // fetch the tile data from the mbtile database
824 try
825 {
826
827 char qrs[2100];
828 sprintf(qrs, "select tile_data, length(tile_data) from tiles where zoom_level = %d AND tile_column=%d AND tile_row=%d", tile->m_zoomLevel, tile->tile_x, tile->tile_y);
829
830 // Compile a SQL query, getting the specific blob
831 SQLite::Statement query(*m_pDB, qrs);
832
833 int queryResult = query.tryExecuteStep();
834 if(SQLITE_DONE == queryResult){
835 tile->m_bAvailable = false;
836 return false; // requested ROW not found, should never happen
837 }
838 else{
839 SQLite::Column blobColumn = query.getColumn(0); // Get the blob
840 const void* blob = blobColumn.getBlob();
841
842 int length = query.getColumn(1); // Get the length
843
844 wxMemoryInputStream blobStream(blob, length);
845 wxImage blobImage;
846
847 blobImage = wxImage(blobStream, m_imageType);
848
849 int blobWidth = blobImage.GetWidth();
850 int blobHeight = blobImage.GetHeight();
851 unsigned char *imgdata = blobImage.GetData();
852
853 if( (m_global_color_scheme != GLOBAL_COLOR_SCHEME_RGB) && (m_global_color_scheme != GLOBAL_COLOR_SCHEME_DAY) ){
854 double dimLevel;
855 switch( m_global_color_scheme ){
856 case GLOBAL_COLOR_SCHEME_DUSK: {
857 dimLevel = 0.8;
858 break;
859 }
860 case GLOBAL_COLOR_SCHEME_NIGHT: {
861 dimLevel = 0.3;
862 break;
863 }
864 default: {
865 dimLevel = 1.0;
866 break;
867 }
868 }
869
870 // for( int iy = 0; iy < blobHeight; iy++ ) {
871 // for( int ix = 0; ix < blobWidth; ix++ ) {
872 // wxImage::RGBValue rgb( blobImage.GetRed( ix, iy ), blobImage.GetGreen( ix, iy ), blobImage.GetBlue( ix, iy ) );
873 // wxImage::HSVValue hsv = wxImage::RGBtoHSV( rgb );
874 // hsv.value = hsv.value * dimLevel;
875 // wxImage::RGBValue nrgb = wxImage::HSVtoRGB( hsv );
876 // blobImage.SetRGB( ix, iy, nrgb.red, nrgb.green, nrgb.blue );
877 // }
878 // }
879
880 for( int j = 0; j < blobHeight*blobWidth; j++ ){
881 unsigned char *d = &imgdata[3*j];
882 wxImage::RGBValue rgb( *d, *(d+1), *(d+2) );
883 wxImage::HSVValue hsv = wxImage::RGBtoHSV( rgb );
884 hsv.value = hsv.value * dimLevel;
885 wxImage::RGBValue nrgb = wxImage::HSVtoRGB( hsv );
886 *d = nrgb.red; *(d+1) = nrgb.green; *(d+2) = nrgb.blue;
887 }
888
889 }
890
891
892 int stride = 4;
893 int tex_w = 256;
894 int tex_h = 256;
895 if( !imgdata )
896 return false;
897 m_imageType = blobImage.GetType();
898
899 unsigned char *teximage = (unsigned char *) malloc( stride * tex_w * tex_h );
900 bool transparent = blobImage.HasAlpha();
901
902 for( int j = 0; j < tex_w*tex_h; j++ ){
903 for( int k = 0; k < 3; k++ )
904 teximage[j * stride + k] = imgdata[3*j + k];
905
906 // Some NOAA Tilesets do not give transparent tiles, so we detect NOAA's idea of blank
907 // as RGB(1,0,0) and force alpha = 0;
908 if( imgdata[3*j] == 1 && imgdata[3*j + 1] == 0 && imgdata[3*j+2] == 0) {
909 teximage[j * stride + 3] = 0;
910 } else {
911 if( transparent ) {
912 teximage[j * stride + 3] = blobImage.GetAlpha(j % tex_w, j / tex_w);
913 } else {
914 teximage[j * stride + 3] = 255;
915 }
916 }
917 }
918
919
920 glGenTextures( 1, &tile->glTextureName );
921 glBindTexture( GL_TEXTURE_2D, tile->glTextureName );
922
923 glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE );
924 glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE );
925 glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR );
926 glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR );
927
928 glTexImage2D( GL_TEXTURE_2D, 0, GL_RGBA, tex_w, tex_h, 0, GL_RGBA, GL_UNSIGNED_BYTE, teximage );
929
930 free(teximage);
931
932 return true;
933 }
934
935 }
936 catch (std::exception& e)
937 {
938 const char *t = e.what();
939 wxLogMessage("mbtiles exception: %s", e.what());
940 }
941 }
942
943
944
945 return false;
946 }
947
948 class wxPoint2DDouble;
949
GetDoublePixFromLL(ViewPort & vp,double lat,double lon)950 wxPoint2DDouble ChartMBTiles::GetDoublePixFromLL( ViewPort& vp, double lat, double lon )
951 {
952 double easting = 0;
953 double northing = 0;
954 double xlon = lon - eps;
955
956 switch( vp.m_projection_type ) {
957 case PROJECTION_MERCATOR:
958 case PROJECTION_WEB_MERCATOR:
959 default:
960 const double z = WGS84_semimajor_axis_meters * mercator_k0;
961
962 easting = (xlon - vp.clon) * DEGREE * z;
963
964 // y =.5 ln( (1 + sin t) / (1 - sin t) )
965 const double s = sin(lat * DEGREE);
966 const double y3 = (.5 * log((1 + s) / (1 - s))) * z;
967
968 const double s0 = sin(vp.clat * DEGREE);
969 const double y30 = (.5 * log((1 + s0) / (1 - s0))) * z;
970 northing = y3 - y30;
971
972 break;
973 }
974
975 double epix = easting * vp.view_scale_ppm;
976 double npix = northing * vp.view_scale_ppm;
977 double dxr = epix;
978 double dyr = npix;
979
980 // Apply VP Rotation
981 double angle = vp.rotation;
982
983 if( angle ) {
984 dxr = epix * cos( angle ) + npix * sin( angle );
985 dyr = npix * cos( angle ) - epix * sin( angle );
986 }
987
988 //printf(" gdpll: %g %g %g\n", vp.clon, (vp.pix_width / 2.0 ) + dxr, ( vp.pix_height / 2.0 ) - dyr);
989
990 return wxPoint2DDouble(( vp.pix_width / 2.0 ) + dxr, ( vp.pix_height / 2.0 ) - dyr);
991 }
992
RenderTile(mbTileDescriptor * tile,int zoomLevel,const ViewPort & VPoint)993 bool ChartMBTiles::RenderTile( mbTileDescriptor *tile, int zoomLevel, const ViewPort& VPoint)
994 {
995 ViewPort vp = VPoint;
996
997 bool btexture = getTileTexture(tile);
998 if(!btexture) { // failed to load, draw NODTA on the minimum zoom
999 glDisable(GL_TEXTURE_2D);
1000 return false;
1001 }
1002 else{
1003 glEnable(GL_TEXTURE_2D);
1004 glEnable(GL_BLEND);
1005 glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
1006
1007 }
1008
1009 float coords[8];
1010 float texcoords[] = { 0., 1., 0., 0., 1., 0., 1., 1. };
1011
1012 ViewPort mvp = vp;
1013
1014 wxPoint2DDouble p;
1015
1016 p = GetDoublePixFromLL(mvp, tile->latmin, tile->lonmin); coords[0] = p.m_x; coords[1] = p.m_y;
1017 p = GetDoublePixFromLL(mvp, tile->latmax, tile->lonmin); coords[2] = p.m_x; coords[3] = p.m_y;
1018 p = GetDoublePixFromLL(mvp, tile->latmax, tile->lonmax); coords[4] = p.m_x; coords[5] = p.m_y;
1019 p = GetDoublePixFromLL(mvp, tile->latmin, tile->lonmax); coords[6] = p.m_x; coords[7] = p.m_y;
1020
1021 glChartCanvas::RenderSingleTexture(coords, texcoords, &vp, 0, 0, 0);
1022
1023 glDisable(GL_BLEND);
1024
1025 return true;
1026 }
1027
1028
RenderRegionViewOnGL(const wxGLContext & glc,const ViewPort & VPoint,const OCPNRegion & RectRegion,const LLRegion & Region)1029 bool ChartMBTiles::RenderRegionViewOnGL(const wxGLContext &glc, const ViewPort& VPoint, const OCPNRegion &RectRegion, const LLRegion &Region)
1030 {
1031 // Do not render if significantly underzoomed
1032 if( VPoint.chart_scale > (20 * OSM_zoomScale[m_minZoom])){
1033 if(m_nTiles > 500){
1034 return true;
1035 }
1036 }
1037
1038 ViewPort vp = VPoint;
1039
1040 OCPNRegion screen_region(wxRect(0, 0, vp.pix_width, vp.pix_height));
1041 LLRegion screenLLRegion = vp.GetLLRegion( screen_region );
1042 LLBBox screenBox = screenLLRegion.GetBox();
1043
1044 if((m_LonMax - m_LonMin) > 180){ // big chart
1045 LLRegion validRegion = m_minZoomRegion;
1046 validRegion.Intersect(screenLLRegion);
1047 glChartCanvas::SetClipRegion(vp, validRegion);
1048 }
1049 else
1050 glChartCanvas::SetClipRegion(vp, m_minZoomRegion);
1051
1052
1053 /* setup opengl parameters */
1054 glEnable( GL_TEXTURE_2D );
1055
1056 int viewZoom = m_maxZoom;
1057 double zoomMod = 2.0; // decrease to get more detail, nominal 4?, 2 works OK for NOAA.
1058
1059 for(int kz=m_minZoom ; kz <= 19 ; kz++){
1060 double db_mpp = OSM_zoomMPP[kz];
1061 double vp_mpp = 1. / vp.view_scale_ppm;
1062
1063 if(db_mpp < vp_mpp * zoomMod){
1064 viewZoom = kz;
1065 break;
1066 }
1067 }
1068
1069 viewZoom = wxMin(viewZoom, m_maxZoom);
1070 //printf("viewZoomCalc: %d %g %g\n", viewZoom, VPoint.view_scale_ppm, 1. / VPoint.view_scale_ppm);
1071
1072 int zoomFactor = m_minZoom;
1073
1074 // DEBUG TODO Show single zoom
1075 //zoomFactor = 5; //m_minZoom;
1076 //viewZoom = zoomFactor;
1077
1078 int maxrenZoom = m_minZoom;
1079
1080 LLBBox box = Region.GetBox();
1081
1082 // if the full screen box spans IDL,
1083 // we need to render the entire screen in two passes.
1084 bool btwoPass = false;
1085 if( ((screenBox.GetMinLon() < -180) && (screenBox.GetMaxLon() > -180)) ||
1086 ((screenBox.GetMinLon() < 180) && (screenBox.GetMaxLon() > 180)) ){
1087 //printf("\nTwoPass\n");
1088 btwoPass = true;
1089 box = screenBox;
1090 }
1091
1092
1093 while(zoomFactor <= viewZoom){
1094 //printf("zoomFactor: %d viewZoom: %d\n", zoomFactor, viewZoom);
1095 mbTileZoomDescriptor *tzd = m_tileArray[zoomFactor - m_minZoom];
1096
1097
1098 // Get the tile numbers of the box corners of this render region, at this zoom level
1099 int topTile = wxMin(tzd->tile_y_max, lat2tiley(box.GetMaxLat(), zoomFactor));
1100 int botTile = wxMax(tzd->tile_y_min, lat2tiley(box.GetMinLat(), zoomFactor));
1101 int leftTile = long2tilex(box.GetMinLon(), zoomFactor);
1102 int rightTile = long2tilex(box.GetMaxLon(), zoomFactor);
1103
1104 if(btwoPass){
1105 leftTile = long2tilex(-180+eps, zoomFactor);
1106 rightTile = long2tilex(box.GetMaxLon(), zoomFactor);
1107 vp = VPoint;
1108 if(vp.clon > 0)
1109 vp.clon -= 360;
1110
1111 }
1112 else
1113 vp = VPoint;
1114
1115 //botTile -= 1;
1116 topTile += 1;
1117
1118 //printf("limits: {%d %d} {%d %d}\n", botTile, topTile, leftTile, rightTile);
1119
1120 for(int i=botTile ; i < topTile ; i++){
1121 if( (i > tzd->tile_y_max) || (i < tzd->tile_y_min) )
1122 continue;
1123
1124 for(int j = leftTile ; j < rightTile+1 ; j++){
1125 if( (tzd->tile_x_max >= tzd->tile_x_min) && ((j > tzd->tile_x_max) || (j < tzd->tile_x_min)) )
1126 continue;
1127
1128 unsigned int index = ((i- tzd->tile_y_min) * (tzd->nx_tile + 1)) + j;
1129 //printf("pass 1: %d %d %d\n", zoomFactor, i, j);
1130 mbTileDescriptor *tile = NULL;
1131
1132 if(tzd->tileMap.find(index) != tzd->tileMap.end())
1133 tile = tzd->tileMap[index];
1134 if(NULL == tile){
1135 tile = new mbTileDescriptor;
1136 tile->tile_x = j;
1137 tile->tile_y = i;
1138 tile->m_zoomLevel = zoomFactor;
1139 tile->m_bAvailable = true;
1140
1141 tzd->tileMap[index] = tile;
1142 }
1143
1144 if(!tile->m_bgeomSet){
1145
1146 tile->lonmin = round(tilex2long(tile->tile_x, zoomFactor)/eps)*eps;
1147 tile->lonmax = round(tilex2long(tile->tile_x + 1, zoomFactor)/eps)*eps;
1148 tile->latmin = round(tiley2lat(tile->tile_y - 1, zoomFactor)/eps)*eps;
1149 tile->latmax = round(tiley2lat(tile->tile_y, zoomFactor)/eps)*eps;
1150
1151 tile->box.Set(tile->latmin, tile->lonmin, tile->latmax, tile->lonmax);
1152 tile->m_bgeomSet = true;
1153 }
1154
1155 if(!Region.IntersectOut(tile->box))
1156 {
1157 if(RenderTile(tile, zoomFactor, vp))
1158 maxrenZoom = zoomFactor;
1159 }
1160 }
1161 } // for
1162
1163 // second pass
1164 if(btwoPass){
1165 vp = VPoint;
1166 if(vp.clon < 0)
1167 vp.clon += 360;
1168
1169 // Get the tile numbers of the box corners of this render region, at this zoom level
1170 int topTile = wxMin(tzd->tile_y_max, lat2tiley(box.GetMaxLat(), zoomFactor));
1171 int botTile = wxMax(tzd->tile_y_min, lat2tiley(box.GetMinLat(), zoomFactor));
1172 int leftTile = long2tilex(box.GetMinLon(), zoomFactor);
1173 int rightTile = long2tilex(-180-eps/*box.GetMaxLon()*/, zoomFactor);
1174
1175 if(rightTile < leftTile)
1176 rightTile = leftTile;
1177 topTile += 1;
1178
1179
1180 for(int i=botTile ; i < topTile ; i++){
1181 for(int j = leftTile ; j < rightTile+1 ; j++){
1182 unsigned int index = ((i- tzd->tile_y_min) * (tzd->nx_tile + 1)) + j;
1183
1184 mbTileDescriptor *tile = NULL;
1185
1186 //printf("pass 2: %d %d %d\n", zoomFactor, i, j);
1187
1188 if(tzd->tileMap.find(index) != tzd->tileMap.end())
1189 tile = tzd->tileMap[index];
1190 if(NULL == tile){
1191 tile = new mbTileDescriptor;
1192 tile->tile_x = j;
1193 tile->tile_y = i;
1194 tile->m_zoomLevel = zoomFactor;
1195 tile->m_bAvailable = true;
1196
1197 tzd->tileMap[index] = tile;
1198 }
1199
1200 if(!tile->m_bgeomSet){
1201
1202 tile->lonmin = round(tilex2long(tile->tile_x, zoomFactor)/eps)*eps;
1203 tile->lonmax = round(tilex2long(tile->tile_x + 1, zoomFactor)/eps)*eps;
1204 tile->latmin = round(tiley2lat(tile->tile_y - 1, zoomFactor)/eps)*eps;
1205 tile->latmax = round(tiley2lat(tile->tile_y, zoomFactor)/eps)*eps;
1206
1207 tile->box.Set(tile->latmin, tile->lonmin, tile->latmax, tile->lonmax);
1208 tile->m_bgeomSet = true;
1209 }
1210
1211
1212 if(!Region.IntersectOut(tile->box))
1213 RenderTile(tile, zoomFactor, vp);
1214 }
1215 } // for
1216 }
1217
1218 zoomFactor++;
1219 }
1220
1221 glDisable(GL_TEXTURE_2D);
1222
1223 m_zoomScaleFactor = 2.0 * OSM_zoomMPP[maxrenZoom] * VPoint.view_scale_ppm;
1224
1225 glChartCanvas::DisableClipRegion();
1226
1227 return true;
1228 }
1229
1230
RenderRegionViewOnDC(wxMemoryDC & dc,const ViewPort & VPoint,const OCPNRegion & Region)1231 bool ChartMBTiles::RenderRegionViewOnDC(wxMemoryDC& dc, const ViewPort& VPoint, const OCPNRegion &Region)
1232 {
1233 gFrame->GetPrimaryCanvas()->SetAlertString( _("MBTile requires OpenGL to be enabled"));
1234
1235 return true;
1236 }
1237